;;; -*- Mode: Lisp -*-
;;;
;;; $Header: /home/gene/library/website/docsrc/linear/RCS/matrix.lisp,v 395.1 2008/04/20 17:25:55 gene Exp $
;;;
;;; Copyright (c) 2007 Gene Michael Stover.  All rights reserved.
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU Lesser 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 St, Fifth Floor, Boston, MA  02110-1301
;;; USA
;;;

(defpackage "COM.CYBERTIGGYR.GENE.MATRIX"
  (:use "COMMON-LISP")
  (:import-from "CYBERTIGGYR-TEST" "DEFTEST")
  (:export "IS-SAME-SIZE"
	   "MATRIXP"
	   "MATRIX"
	   "MREF"))
(in-package "COM.CYBERTIGGYR.GENE.MATRIX")

(deftest test0000 ()
  "Null test.  Always succeeds."
  'test0000)

(deftype matrix () `(simple-array real (* *)))

(deftest test0011 ()
  "Test that an array we create with the standard function MAKE-ARRAY
is a MATRIX."
  (typep (make-array '(2 3) :element-type 'real :initial-element 0)
	 'matrix))

(defun matrixp (x)
  "Return true if & only if X is of type MATRIX."
  (typep x 'matrix))

(deftest test0012 ()
  "Like TEST0011 but uses MATRIXP."
  (matrixp (make-array '(2 3) :element-type 'real :initial-element 0)))

(deftest test0013 ()
  "Test that a vector is not a MATRIX."
  (not (matrixp (make-array 3 :element-type 'real :initial-element 0))))

(defun mref (m i j)
  "Return value of element at row I, column J in matrix M.
Assumes I & J are within range; else, you get an error."
  (declare (type matrix m) (type fixnum i j))
  (assert (< -1 i (array-dimension m 0)))
  (assert (< -1 j (array-dimension m 1)))
  (aref m i j))

(deftest test0021 ()
  "Test MREF."
  (zerop (mref (make-array '(3 3) :element-type 'real :initial-element 0)
	       1 2)))

(deftest test0022 ()
  "Test MREF."
  (eql 2/3 (mref (make-array '(3 3) :element-type 'real :initial-element 2/3)
		 1 2)))

(defun mset (m i j value)
  "Assign VALUE to the element at row I, column J in matrix M.  Assume
I & J are within range.  Assume VALUE is of type REAL."
  (declare (type matrix m) (type fixnum i j) (type real value))
  (assert (< -1 i (array-dimension m 0)))
  (assert (< -1 j (array-dimension m 1)))
  (setf (aref m i j) value))

(deftest test0031 ()
  "Test that MSET returns its last argument."
  (let ((m (make-array '(3 3) :element-type 'real :initial-element 0)))
    (eql 2/3 (mset m 1 1 2/3))))

(deftest test0032 ()
  "Test that an element MSET returns the same value from MREF."
  (let ((m (make-array '(3 3) :element-type 'real :initial-element 0)))
    (mset m 1 1 2/3)
    (if (eql 2/3 (mref m 1 1))
	(progn
	  (mset m 1 1 -1)
	  (if (eql -1 (mref m 1 1))
	      'test0032                 ; good
	    (progn
	      (format t "~&~A: MREF did not return the expected value."
		      'test0032)
	      nil)))
      ;; else
      (progn 
	(format t "~&~A: MREF did not return the expected value." 'test0032)
	nil))))
  
(defsetf mref (m i j) (value)
  `(mset ,m ,i ,j ,value))

(deftest test0041 ()
  "Test that SETF MREF returns its final value when we set one element."
  (let ((a (make-array '(3 3) :element-type 'real :initial-element 0)))
    (eql 2/3 (setf (mref a 1 1) 2/3))))

(deftest test0042 ()
  "Test that SETF MREF returns its final value when we set three
elements."
  (let ((a (make-array '(3 3) :element-type 'real :initial-element 0)))
    (eql -2 (setf (mref a 0 0) 2/3
		  (mref a 0 1) 1
		  (mref a 1 0) -2))))

(defun is-same-size (a b)
  "Return true if & only if A & B have the same dimensions."
  (declare (type matrix a b))
  (equal (array-dimensions a) (array-dimensions b)))

(deftest test0051 ()
  "Test that IS-SAME-SIZE returns true for two matrices with the same
size."
  (let ((a (make-array '(2 3) :element-type 'real :initial-element 1))
	(b (make-array '(2 3) :element-type 'real :initial-element 2)))
    (is-same-size a b)))

(deftest test0052 ()
  "Test that IS-SAME-SIZE returns false for two matrices with the same
size."
  (let ((a (make-array '(3 2) :element-type 'real :initial-element 1))
	(b (make-array '(2 3) :element-type 'real :initial-element 2)))
    (not (is-same-size a b))))

(defun foreach (fn m)
  "Call FN for each element in M.  Call FN with the row & column
of the element.  Return the final value which FN returned."
  (declare (type function fn) (type matrix m))
  (let ((x nil)
	(rcount (array-dimension m 0))
	(ccount (array-dimension m 1)))
    (do ((r 0 (1+ r)))
	((>= r rcount) x)
	(declare (type fixnum r))
	(do ((c 0 (1+ c)))
	    ((>= c ccount))
	    (declare (type fixnum c))
	    (setq x (funcall fn m r c))))))

(deftest test0061 ()
  "Test FOREACH on a 3-by-3 matrix."
  (let ((x (foreach (constantly 123) (make-array '(3 3) :element-type 'real
						 :initial-element -1))))
    (eql x 123)))

(deftest test0062 ()
  "Test that FOREACH returns the final value returned by FN, even
if FN doesn't return the same value every time."
  ;; We'll call FOREACH on a 2-by-3 matrix.  It should call FN
  ;; six times.  FN increments a counter each time.
  (let ((count 0)
	(m (make-array '(2 3) :element-type 'real :initial-element 1)))
    (let ((x (foreach #'(lambda (m i j)
			  (declare (ignore m i j))
			  (incf count))
		      m)))
      (and (eql 6 x) (eql 6 count)))))

(defun map-matrix (fn &rest mlst)
  "Apply function FN to all corresponding elements of M & the other
matrices in LST.  Return final value which was returned by FN."
  (declare (type function fn) (type list mlst))
  (assert (every #'matrixp mlst))
  (assert (every #'(lambda (m) (is-same-size m (first mlst))) (rest mlst)))
  (let ((x nil)
	(rcount (array-dimension (first mlst) 0))
	(ccount (array-dimension (first mlst) 1)))
    (do ((r 0 (1+ r)))
	((>= r rcount) x)
	(declare (type fixnum r))
	(do ((c 0 (1+ c)))
	    ((>= c ccount))
	    (declare (type fixnum c))
	    (setq x (apply fn r c (mapcar #'(lambda (y)
					      (declare (type matrix y))
					      (mref y r c))
					  mlst)))))))

(deftest test0071 ()
  "Test that MAP-MATRIX iterates over each element & returns the
last value returned by FN."
  ;; An easy way to ensure that MAP-MATRIX visits every element is
  ;; to sum the elements.  We also count the number of visits.
  (let ((count 0)
	(sum 0)
	(a (make-array '(1 2) :element-type 'real :initial-element 0))
	(b (make-array '(1 2) :element-type 'real :initial-element 0)))
    (setf (mref a 0 0) 1 (mref a 0 1) 2
	  (mref b 0 0) 3 (mref b 0 1) 4)
    (let ((x (map-matrix #'(lambda (r c vala valb)
			     (declare (ignore r c) (type real vala valb))
			     (incf sum (+ vala valb))
			     (incf count))
			 a b)))
      (and (eql 2 x)
	   (eql 10 sum)))))

(defun make-matrix (rows columns &key (initial-element 0))
  "Return a new matrix with the indicated number of rows & columns.
Every element will be initialized to INITIAL-ELEMENT, which defaults
to zero.  INITIAL-ELEMENT must be of type REAL."
  (declare (type fixnum rows columns) (type real initial-element))
  (make-array (list rows columns) :element-type 'real
	      :initial-element initial-element))

(deftest test0081 ()
  "Test MAKE-MATRIX for a matrix with 2 rows, 1 column & a default
initial element."
  (let ((m (make-matrix 2 1)))
    (and (matrixp m)
	 (equal (array-dimensions m) '(2 1))
	 (eql (mref m 0 0) 0)
	 (eql (mref m 1 0) 0))))

(deftest test0082 ()
  "Test MAKE-MATRIX for a matrix with 2 rows, 3 columns & an initial
element of -3."
  (let* ((m (make-matrix 2 3 :initial-element -3))
	 (is-good (and (matrixp m)
		       (equal (array-dimensions m) '(2 3))
		       (eql (mref m 0 0) -3)
		       (eql (mref m 0 1) -3)
		       (eql (mref m 0 2) -3)
		       (eql (mref m 1 0) -3)
		       (eql (mref m 1 1) -3)
		       (eql (mref m 1 2) -3))))
    (unless is-good (print m))
    is-good))

(defun equal-matrix (a b)
  "Return true if & only if matrices A & B are of the same size & their
corresponding elements are EQL."
  (declare (type matrix a b))
  (and (equal (array-dimensions a) (array-dimensions b))
       (let ((is-equal t))
	 (map-matrix #'(lambda (r c vala valb)
			 (declare (ignore r c) (type real vala valb))
			 (unless (eql vala valb) (setq is-equal nil))
			 is-equal)
		     a b))))

(deftest test0091 ()
  "Test EQUAL-MATRIX on two matrices which are equal."
  (equal-matrix #2A((1 2 3) (4 5 6))
		#2A((1 2 3) (4 5 6))))

(deftest test0092 ()
  "Test EQUAL-MATRIX on two matrices which have the same size, but some
of their elements differ."
  (not (equal-matrix #2A((-1 2 3) (4 5 6))
		     #2A(( 1 2 3) (4 5 6)))))

(deftest test0093 ()
  "Test EQUAL-MATRIX on two matrices which have different sizes."
  (not (equal-matrix #2A((1 2 3) (4 5 6))
		     #2A((1 2 3)))))
	     
(defun copy-matrix (src)
  "Return a copy of matrix M."
  (declare (type matrix src))
  (let ((dst (apply #'make-matrix (array-dimensions src))))
    ;; FOREACH returns the last value returned by lambda.  Our
    ;; lambda always returns DST.  So FOREACH returns DST, & we
    ;; (COPY-MATRIX) return DST.
    (foreach #'(lambda (m r c)
		 (declare (ignore m) (type fixnum r c))
		 (setf (mref dst r c) (mref src r c))
		 dst)
	     src)))

(defun test0101 ()
  "Test COPY-MATRIX.  It should return a matrix which is not EQ to the
original matrix but whose elements are equivalent to those of the
original."
  (let ((a (make-matrix 3 4)))
    ;; Set a couple of A's elements so we know COPY-MATRIX copied
    ;; everything.
    (setf (mref a 0 3) 1
	  (mref a 1 2) 2
	  (mref a 2 1) 3)
    (let ((b (copy-matrix a)))
      (and (not (eq a b))
	   (equal-matrix a b)))))

(defun test0102 ()
  "Test COPY-MATRIX on a 1-by-1 matrix."
  (let ((a (make-matrix 1 1 :initial-element 3)))
    (let ((b (copy-matrix a)))
      (and (not (eq a b))
	   (equal-matrix a b)))))

(defun madd (a &rest mlst)
  "Return a new matrix which is the sum of other matrix A & other
matrices.  All of the matrices must have the same dimensions."
  (declare (type matrix a) (type list mlst))
  (assert (every #'matrixp mlst))
  (assert (every #'(lambda (b) (is-same-size a b)) mlst))
  (let ((c (copy-matrix a)))
    ;; Because our lambda always returns C, MAP-MATRIX will always
    ;; return C.  And since MAP-MATRIX is the last thing we do, we'll
    ;; return C.
    (apply
     #'map-matrix
     #'(lambda (i j &rest vlst)
	 (declare (type fixnum i j) (type list vlst))
	 (assert (every #'realp vlst))
	 ;; If VLST was long, it'd be safer to use REDUCE
	 ;; instead of APPLY, but VLST will almost always
	 ;; have only 1 element, & I doubt that it'll ever
	 ;; have more than two or three.  So APPLY is safe
	 ;; in practice, & I suspect it is faster than
	 ;; REDUCE.
	 (incf (mref c i j) (apply #'+ vlst))
	 c)
     mlst)))

(deftest test0111 ()
  "Test MADD on a couple of 2-by-3 matrices."
  (let* ((a #2A((1 2 3) ( 4  5  6)))
	 (b #2A((7 8 9) (10 11 12)))
	 (c (madd a b))
	 (expect #2A(( 8 10 12) (14 16 18)))
	 (is-good (equal-matrix c expect)))
    (unless is-good (print c))
    is-good))

(defun msub (a &rest mlst)
  "Return a new matrix which is the difference of matrix A & other
matrices.  All of the matrices must have the same dimensions."
  (declare (type matrix a) (type list mlst))
  (assert (every #'matrixp mlst))
  (assert (every #'(lambda (b) (is-same-size a b)) mlst))
  (let ((c (copy-matrix a)))
    ;; Because our lambda always returns C, MAP-MATRIX will always
    ;; return C.  And since MAP-MATRIX is the last thing we do, we'll
    ;; return C.
    (apply
     #'map-matrix
     #'(lambda (i j &rest vlst)
	 (declare (type fixnum i j) (type list vlst))
	 (assert (every #'realp vlst))
	 ;; If VLST was long, it'd be safer to use REDUCE
	 ;; instead of APPLY, but VLST will almost always
	 ;; have only 1 element, & I doubt that it'll ever
	 ;; have more than two or three.  So APPLY is safe
	 ;; in practice, & I suspect it is faster than
	 ;; REDUCE.
	 (decf (mref c i j) (apply #'+ vlst))
	 c)
     mlst)))

(deftest test0113 ()
  "Test MADD on a couple of 2-by-3 matrices."
  (let* ((a #2A((12 11 10) ( 9  8  7)))
	 (b #2A(( 6  5  4) ( 3  2  1)))
	 (c (msub a b))
	 (expect #2A(( 6  6  6) ( 6  6  6)))
	 (is-good (equal-matrix c expect)))
    (unless is-good (print c))
    is-good))

(deftest test0114 ()
  "Test MADD on  three of 2-by-3 matrices.  It's important to ensure
that (MSUB A B C) is A - B - C."
  (let* ((a #2A((12 11 10) ( 9  8  7)))
	 (b #2A(( 6  5  4) ( 3  2  1)))
	 (c #2A(( 5  4  3) ( 2  1  0)))
	 (d (msub a b c))
	 (expect #2A(( 1  2  3) ( 4  5  6)))
	 (is-good (equal-matrix d expect)))
    (unless is-good (print c))
    is-good))

(defun mmul-sum (i j a b)
  "Return sum of A[i][k] * B[k][j] for k = 0 ... m-1, where m is the
number of columns in A (& also the number of rows in B)."
  (declare (type fixnum i j) (type matrix a b))
  ;; These variable names are from Definition 6.6 on page 336 of The
  ;; Book (Burden & Faires).
  (assert (eql (array-dimension a 1) (array-dimension b 0)))
  (symbol-macrolet ((m (array-dimension a 1)))
    (do ((k 0 (1+ k))
	 (sum 0 (+ (* (mref a i k) (mref b k j)) sum)))
	((>= k m) sum)
	(declare (type fixnum k) (type real sum)))))

(defun mmul (a b)
  "Return new matrix which is product of matrices A & B."
  (declare (type matrix a b))
  (assert (eql (array-dimension a 1) (array-dimension b 0)))
  ;; These variable names are from Definition 6.6 on page 336 of The
  ;; Book (Burden & Faires).
  (symbol-macrolet ((n (array-dimension a 0))
		    (m (array-dimension a 1))
		    (p (array-dimension b 1)))
    (let ((c (make-matrix n p)))
      (declare (type matrix c))
      (do ((i 0 (1+ i)))
	  ((>= i n))
	  (declare (type fixnum i))
	  (do ((j 0 (1+ j)))
	      ((>= j p))
	      (declare (type fixnum j))
	      (setf (mref c i j) (mmul-sum i j a b))))
      c)))

(deftest test0121 ()
  "Test MMUL on a hard-coded input.  Input A is a 3-by-3 matrix.  Input
B is a 3-by-3 matrix.  The output should be a 3-by-3 matrix.  These
values are from the A & D matrices of Example 3 on page 337 of The Book."
  (let ((a #2A((2 1 -1) (3 1 2) (0 -2 -3)))
	(b #2A((1 -1 1) (2 -1 2) (3 0 3)))
	(expect #2A((1 -3 1) (11 -4 11) (-13 2 -13))))
    (declare (type matrix a b expect))
    (let ((x (mmul a b)))
      (declare (type matrix x))
      (unless (equal-matrix x expect) (print x))
      (equal-matrix x expect))))

;;; --- end of file ---