;;; Lisplab, level2-constructors.lisp
;;; Possible and impossible ways to create matrices.

;;; 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.

;;; TODO: needs constructors for diagonal matrices.

;;; TODO: specialize convert for standard-class.


(in-package :lisplab)

;;; Creates matrices with general structure, e.g., #md((1 2) (3 4))

;;; TODO: error check is important here!

(set-dispatch-macro-character #\# #\M
  (lambda (stream c1 c2)
    (let* ((s1 (make-string 1)))
      (setf (aref s1 0) (read-char stream))
      (setf s1 (string-capitalize s1)) 
      (let ((type (cond ((string= s1 "D") :d)
			((string= s1 "Z") :z)
			(t :any)))
	    (cont (read stream t nil t)))
	(list 'mmat 
	      (list 'list type :ge :any)
	      (cons 'list (mapcar (lambda (x) (cons 'list x)) cont)))))))
			  
(defmethod mcreate ((m number) &optional (val 0) dim)
  (declare (ignore dim))
  ;; This is not about matrices at all, but is usefull 
  ;; when you use the dotted algebra and is not sure is input is numbers or matrices.
  ;; TODO what the dim, should I use it or ignore it
  val)

(defmethod mcreate* ((m number) &key (value 0) dim element-type structure implementation)
  ;; This is not about matrices at all, but is usefull 
  ;; when you use the dotted algebra and is not sure is input is numbers or matrices.
  ;; TODO what the dim, should I use it or ignore it
  (declare (ignore dim element-type structure implementation))
  value)

(defmethod mcreate ((a matrix-base) &optional (value 0) dim)
  (unless dim
    (setf dim (dim a)))
  (make-matrix-instance (class-of a) dim value)) 

(defmethod mcreate* ((a matrix-base) &key (value 0) dim element-type structure implementation)
  (unless dim
    (setf dim (dim a)))
  (unless element-type 
    (setf element-type (element-type-spec a)))
  (unless structure 
    (setf structure (structure-spec a)))  
  (unless implementation
    (setf implementation (implementation-spec a)))
  (make-matrix-instance (list element-type structure implementation) dim value)) 
(defmethod convert ((x matrix-base) type)
  (let ((y (make-matrix-instance type (dim x) 0)))
    (copy-contents x y)
    y))

(defun fill-matrix-with-list (m x)  
  "Helper function for convert."
  (let* ((rows (rows m))
	 (cols (cols m)))
    (do ((xx x (cdr xx))
	 (i 0 (1+ i)))
	((= i rows))
      (do ((yy (car xx) (cdr yy))
	   (j 0 (1+ j)))
	  ((= j cols))
	(setf (mref m i j) (car yy))))
    m))

(defmethod mnew (type value rows &optional cols) 
  (make-matrix-instance type 
			(if cols 
			    (list rows cols)
			    rows) 
			value))

(defun mmat (type x)
  "Creates a matrix from the supplied contents."
  (convert x type))

(defun mcol (type &rest args)
  "Creates a column matrix."
  (mmat type (mapcar #'list args)))

(defun mrow (type &rest args)
  "Creates a row matrix."
  (mmat type (list args)))

;;; Constructors for matrix-dge

(defun drandom (rows cols)
  "Creates a double matrix with random element between 0 and 1."
  (mmap t #'random (dnew 1d0 rows cols)))

(defun dmat (args)
  "Creates a matrix-dge from supplied contents."
  (convert args 'matrix-dge))

(defun dcol (&rest args)
  "Creates a matrix-dge column matrix."
  (apply #'mcol 'matrix-dge args))

(defun drow (&rest args)
  "Creates a matrix-dge row matrix."
  (apply #'mrow 'matrix-dge args))

(defun dnew (value rows &optional (cols 1))
  "Creates a matrix-dge matrix"
  (mnew 'matrix-dge value rows cols))

(defun drange (n from to &optional (shift 0))
  "Creates a column vector of length n, with elements of equal spacing 
between from and to. The shift is the a number between 0 and 1 and 
shifts the start position.

For example: (drange 4 0 1) -> 0 1 2 3, while 
             (drange 4 0 1 0.5) -> 0.5 1.5 2.5 3.5."
  (let* ((x (dnew 0 n 1))
	 (store (vector-store x))
	 (dx (/ (- (to-df to) (to-df from)) (to-df N)))
	 (shift (* dx (to-df shift))))
    (declare (type type-blas-store store)
	     (type type-blas-idx n)
	     (type double-float dx shift))
    (do ((i 0 (1+ i))
	 (v (to-df from) (+ v dx)))
	((>= i n) x)
      (declare (type double-float v))
      (setf (aref store i) (+ shift v)))))

(defun dgrid (xv yv)
  ;; TODO: change name to mgrid, since it can take type from xv and yv.
  "Creates grid matrices from input vectors. Input are the x and y vectors 
and outputs are a list of x and y matrices. The input vectors are  
typically created with drange."
   (let* ((r (size xv))
	  (c (size yv))
	  (x (dnew 0 r c))
	  (y (dnew 0 r c))	  
	  (xv* (vector-store xv))
	  (yv* (vector-store yv))
	  (x* (vector-store x))
	  (y* (vector-store y)))
    (declare (type type-blas-store xv* yv* x* y*)
             (type type-blas-idx r c))
    (dotimes (i r)
      (dotimes (j c)    
        (let ((k (column-major-idx i j r)))
          (declare (type type-blas-idx k))
          (setf (aref x* k) (aref xv* i)
                (aref y* k) (aref yv* j)))))
    (list x y)))

;;; Constructors for matrix-zge

(defun zmat (args)
  "Creates a matrix-zge from the supplied contents."
  (convert args 'matrix-zge))

(defun zcol (&rest args)
  "Creates a matrix-zge column matrix."
  (apply #'mcol 'matrix-zge args))

(defun zrow (&rest args)
  "Creates a matrix-zge row matrix."
  (apply #'mrow 'matrix-zge args))

(defun znew (value rows &optional (cols 1))
  "Creates a matrix-zge matrix"
  (mnew 'matrix-zge value rows cols))

(defun zgrid (xv yv)
  (let* ((r (size xv))
         (c (size yv))
         (x (znew 0 r c))
         (y (znew 0 r c)))
    (dotimes (i r)
      (dotimes (j c)
        (setf (mref x i j) (vref xv i)
              (mref y i j) (vref yv j))))
    (list x y)))


;;; Function matrix

(defun funmat (dim fun)
  "Creates a read only function matrix with column major order."
  (let ((rows (first dim))
	(cols (second dim)))
    (make-instance 'function-matrix 
		   :rows rows
		   :cols cols
		   :mref (lambda (self i j) 	
			   (declare (ignore self))
			   (funcall fun i j))
		   :vref (lambda (self i)
			   (declare (ignore self))
			   (multiple-value-bind (c r)
			       (floor i rows)
			     (funcall fun r c))))))

(defun fmat (type dim fun)
  "Creates a matrix of of type type, dim dim from the function definition. 
Row major order"
  (convert (funmat dim fun)
	    type))

(defmethod mcreate ((a function-matrix)
		    &optional (value 0) (dim (dim a)))
  (make-matrix-instance
   (list (element-type-spec a)
	 (structure-spec a)
	 :any)
   dim
   value))