;;; Lisplab, level2-generic.lisp
;;; Level2, non-specialized methods.

;;; Copyright (C) 2009 Joern Inge Vestgaarden
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU General Public License as published by
;;; the Free Software Foundation; either version 2 of the License, or
;;; (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License along
;;; with this program; if not, write to the Free Software Foundation, Inc.,
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.


;;; Implementation principles:
;;; - all operators in this film should specialize for matrix-base and only
;;;   assume level0 and level1 generic function (mref, vref, size, dim, etc.)
;;; - The methods in this file should not assume anything about implementation of
;;;   the matrices.
;;; - The methods in this file should be as short and clean as possible.
;;; - Avoid optimizations (Exept: call other level2 functions, such as mmap, as much as possible.)
;;;


(in-package :lisplab)

;;; For general vector

(defmethod vdot ((a vector-base) (b vector-base))
  (msum (.* a b)))

(defmethod vcross :before ((a vector-base) (b vector-base))
  (assert (= (size a) (size b) 3)))

(defmethod vcross ((a vector-base) (b vector-base))
  (let ((out (mcreate a)))
    (setf (vref out 0) (.- (.* (vref a 1) (vref b 2))
			   (.* (vref a 2) (vref b 1)))
	  (vref out 1) (.- (.* (vref a 2) (vref b 0))
			   (.* (vref a 0) (vref b 2)))
	  (vref out 2) (.- (.* (vref a 0) (vref b 1))
			   (.* (vref a 1) (vref b 0))))
    out))

(defmethod vnorm ((a vector-base))
  (.sqrt (vdot (.conj a) a)))

;;; Vector operations (ignore structure)

(defmethod copy ((a vector-base))
  (let ((x (make-matrix-instance (class-of a) (dim a) 0)))
    (dotimes (i (size x))
      (setf (vref x i) (vref a i)))
    x))

(defmethod mmap ((type (eql t)) f (a vector-base) &rest args)
  "Maps with output type given by first matrix."
  (apply #'mmap (type-of a) f a args)) 

(defmethod mmap ((b (eql nil)) f (a vector-base) &rest args)  
  (cond ((not args)
	 (dotimes (i (size a))
	   (funcall f (vref a i))))
	((not (cdr args))
	 (let ((c (car args)))
	   (dotimes (i (size a))
	     (funcall f (vref a i) (vref c i)))))
	(t (dotimes (i (size a))
	     (apply f (vref a i) 
		    (mapcar (lambda (x) 
			      (vref x i)) 
			    args)))))
  nil)

(defmethod mmap ((type symbol) f (a vector-base) &rest args)
  (apply #'mmap-into (make-matrix-instance type (dim a) 0)  f a args)) 

(defmethod mmap ((type list) f (a vector-base) &rest args)
  ;; The type here is a spec
  (apply #'mmap-into (make-matrix-instance type (dim a) 0)  f a args)) 

;; TODO map of matrix desciptions
(defmethod mmap-into ((b vector-base) f (a vector-base) &rest args)  
  (cond ((not args)
	 (dotimes (i (size a))
	   (setf (vref b i) (funcall f (vref a i)))))
	((not (cdr args))
	 (let ((c (car args)))
	   (dotimes (i (size a))
	     (setf (vref b i) (funcall f (vref a i) (vref c i))))))
	(t (dotimes (i (size a))
	     (setf (vref b i) (apply f (vref a i) 
				     (mapcar (lambda (x) 
					       (vref x i)) 
					     args))))))
  b)

(defmethod msum ((m vector-base))
  (let ((sum 0))
    (dotimes (i (size m))
      (setf sum (.+ sum (vref m i))))
    sum))

(defmethod mmax ((m vector-base))
  (let ((max (vref m 0))
	(idx 0))
    (dotimes (i (size m))
      (when (.> (vref m i) max)
	(setf max (vref m i)
	      idx i)))
    (values max idx)))

(defmethod mmin ((m vector-base))
  (let ((min (vref m 0))
	(idx 0))
    (dotimes (i (size m))
      (when (.< (vref m i) min)
	(setf min (vref m i)
	      idx i)))
    (values min idx)))

(defmethod mabsmax ((m vector-base))
  (let ((max (vref m 0))
	(idx 0))
    (dotimes (i (size m))
      (when (.> (abs (vref m i)) (abs max))
	(setf max (vref m i)
	      idx i)))
    (values max idx)))

(defmethod mabsmin ((m vector-base))
  (let ((min (vref m 0))
	(idx 0))	  
    (dotimes (i (size m))
      (when (.< (abs (vref m i)) (abs min))
	(setf min (vref m i)
	      idx i)))
    (values min idx)))

(defmethod mminmax ((m vector-base))
  (let ((max (vref m 0))
	(min (vref m 0)))    
    (dotimes (i (size m))
      (when (.> (vref m i) max)
	(setf max (vref m i)))
      (when (.< (vref m i) min)
	(setf min (vref m i))))
    (list min max)))

(defmethod mfill ((a vector-base) val)
  (dotimes (i (size a))
    (setf (vref a i) val))
  val)