;;; Lisplab, matrix2-generic.lisp
;;; Level2, non-specialized matrix 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)

;;; This is OK, but could be optimzied!
(defmacro w/mat (a args &body body)
  (let ((a2 (gensym))
	(x (first args))
	(i (second args))
	(j (third args)))
  `(let ((,a2 ,a))
     (dotimes (,i (rows ,a2))
       (dotimes (,j (cols ,a2))
	 (let ((,x (mref ,a2 ,i ,j)))
	   (setf (mref ,a2 ,i ,j)
		 ,@body))))
     ,a2)))

(defmethod copy-contents ((a matrix-base) (b matrix-base) 
			  &optional (converter #'identity))
  (dotimes (i (rows a))
    (dotimes (j (cols a))
      (setf (mref b i j) (funcall converter (mref a i j))))
    b))

(defmethod sub-matrix ((m matrix-base) rr cc)
  (unless (cddr rr)
    (setf rr (cons (car rr) (cons 1 (cdr rr)))))
  (unless (cddr cc)
    (setf cc (cons (car cc) (cons 1 (cdr cc)))))
  (destructuring-bind (r0 r-step r1) rr
    (destructuring-bind (c0 c-step c1) cc
      (when (>= r1 (rows m))
        (setf r1 (1- (rows m))))
      (when (>= c1 (cols m))
        (setf c1 (1- (cols m))))
      (let* ((rows (1+ (floor (- r1 r0) r-step)))
             (cols (1+ (floor (- c1 c0) c-step)))
             (m1 (mcreate m 0 (list rows cols))))
        (dotimes (i rows)
          (dotimes (j cols)
            (setf (mref m1 i j) 
                  (mref m 
			(+ r0 (* r-step i)) 
			(+ c0 (* c-step j))))))
        m1))))

(defmethod get-row ((m matrix-base) row)
  (sub-matrix m (list row row) (list 0 (1- (cols m)))))

(defmethod get-col ((m matrix-base) col)
  (sub-matrix m (list 0 (1- (rows m))) (list col col)))

(defmethod circ-shift ((A matrix-base) shift)
  ;; TODO move to level3
  (let ((B (mcreate A))	 
	(rows (rows A))
	(cols (cols A))
	(dr (first shift))
	(dc (second shift)))
    (dotimes (i rows)
      (dotimes (j cols)
	(setf (mref B (mod (+ i dr) rows) (mod (+ j dc) cols)) 
	      (mref A i j))))
      B))

(defmethod pad-shift ((A matrix-base) shift &optional (value 0))
  ;; TODO move to level3
  (let ((B (mcreate A value))
	(rows (rows A))
	(cols (cols A))
	(dr (first shift))
	(dc (second shift)))
    (loop for i from (max 0 dr) below  (min rows (+ rows dr)) do
	 (loop for j from (max 0 dc) below (min cols (+ cols dc)) do
	      (setf (mref B i j) 
		    (mref A (- i dr) (- j dc)))))
      B))

(defmethod mreverse ((A matrix-base))
  (let ((B (mcreate A))
	(len (size A)))
    (dotimes (i len)
      (setf (vref B (- len i 1))
	    (vref A i)))
    B))

(defmethod export-list ((m matrix-base))
  (let ((res nil))
    (dotimes (i (size m))
      (push (vref m i) res))
    (nreverse res)))

(defmethod import-list ((m matrix-base) list)
  (let ((tmp list))
    (dotimes (i (size m))
      (unless tmp
	(return-from import-list m))
      (setf (vref m i) (car tmp)
	    tmp (cdr tmp)))
    m))

(defmethod reshape ((a matrix-base) shape)
  (let ((B (mcreate a 0 shape)))
    (dotimes (i (size B))
      (setf (vref B i) (vref A i)))
    B))

(defmethod to-vector ((a matrix-base))
  (reshape a (list (size a) 1)))

(defmethod to-matrix ((a matrix-base) rows)
  (reshape a (list rows (/ (size a) rows) 1)))


(defmethod row-swap! ((A matrix-base) i j)
  (dotimes (c (cols A))
    (psetf (mref A i c) (mref A j c)
	   (mref A j c) (mref A i c)))
  A)

(defmethod row-mul! ((A matrix-base) i num)
  (dotimes (c (cols A))
    (setf (mref A i c) (.* num (mref A i c))))
  A)

(defmethod row-add! ((A matrix-base) i j num)
  (dotimes (c (cols A))
    (setf (mref A i c) (.+ (mref A i c) (.* num (mref A j c)))))
  A)

;;; The column operations 

(defmethod col-swap! ((A matrix-base) i j)
  (dotimes (r (rows A))
    (psetf (mref A r i) (mref A r j)
	   (mref A r j) (mref A r i)))
  A)

(defmethod col-smul! ((A matrix-base) i num)
  (dotimes (r (rows A))
    (setf (mref A r i) (.* num (mref A r i))))
  A)
  
(defmethod col-sum ((A matrix-base) i)
  (let ((sum 0))
    (dotimes (r (rows A))
      (setf sum (.+ sum (mref A r i))))
    sum))

(defmethod col-col-mul-sum ((A matrix-base) i (B matrix-base) j)
  (let ((sum 0))
    (dotimes (r (rows A))
      (setf sum (.+ sum (.* (mref A r i)
			    (mref B r j)))))
    sum))