src/data/matrix.lisp

   1 ;;; -*- mode: lisp -*-
   2 ;;;
   3 ;;; Copyright (c) 2007--, by A.J. Rossini <blindglobe@gmail.com>
   4 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
   5
   6 ;;; What is this talk of 'release'? Klingons do not make software
   7 ;;; 'releases'.  Our software 'escapes', leaving a bloody trail of
   8 ;;; designers and quality assurance people in its wake.
   9
  10 (in-package :cls-matrix)
  11
  12 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  13 ;;;;
  14 ;;;;          Array to Row-Major Data Vector Conversion Functions
  15 ;;;;
  16 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  17
  18 (defun array-data-vector (a)
  19 "Args: (a)
  20 Displaces array A to a vector"
  21   (make-array (array-total-size a)
  22               :displaced-to a
  23               :element-type (array-element-type a)))
  24
  25 (defun vector-to-array (v dims)
  26 "Args: (v dims)
  27 Displaces vector V to array with dimensions DIMS"
  28 ;;; Yes, but using row or column first approach?
  29   (make-array dims
  30               :displaced-to v
  31               :element-type (array-element-type v)))
  32
  33 (defun lisp-array-matrix-p (a)
  34   "Verifies that a lisp array could be considered a matrix.  FIXME: Do
  35 we care whether the entries are numeric or even commonly-numeric?"
  36   (check-type a array)
  37   (assert (= (array-rank a) 2) (a) "lisp array is not 2-dim: ~A"  a)
  38   t)
  39
  40 (defun check-square-matrix (a)
  41   (if (and (lisp-array-matrix-p a)
  42            (/= (array-dimension a 0) (array-dimension a 1))
  43            (error "matrix not square - ~s" a))
  44       t))
  45
  46 (defun matrixp (x)
  47 "Args: (x)
  48 Returns T if X is a matrix, NIL otherwise."
  49   (and (arrayp x) (= (array-rank x) 2)))
  50
  51 (defun num-rows (x)
  52 "Args: (x)
  53 Returns number of rows in X."
  54   (array-dimension x 0))
  55
  56 (defun num-cols (x)
  57 "Args: (x)
  58 Returns number of columns in X."
  59   (array-dimension x 1))
  60
  61 (defun matmult (a b &rest args)
  62 "Args: (a b &rest args)
  63 Returns the matrix product of matrices a, b, etc. If a is a vector it is
  64 treated as a row vector; if b is a vector it is treated as a column vector."
  65   (let ((rtype (cond ((and (matrixp a) (matrixp b)) 'matrix)
  66                      ((and (typep a 'sequence) (typep b 'sequence)) 'number)
  67                      ((typep a 'sequence) (if (consp a) 'list 'vector))
  68                      ((typep b 'sequence) (if (consp b) 'list 'vector)))))
  69
  70     (if (typep a 'sequence)
  71       (setf a (vector-to-array (coerce a 'vector) (list 1 (length a)))))
  72     (if (typep b 'sequence)
  73       (setf b (vector-to-array (coerce b 'vector) (list (length b) 1))))
  74     (if (not (= (array-dimension a 1) (array-dimension b 0)))
  75       (error "dimensions do not match"))
  76     (if args
  77       (reduce #'matmult args :initial-value (matmult a b))
  78       (let* ((n (array-dimension a 0))
  79              (m (array-dimension b 1))
  80              (p (array-dimension a 1))
  81              (c (make-array (list n m)))
  82              x)
  83         (declare (fixnum n m p))
  84         (dotimes (i n)
  85                  (declare (fixnum i))
  86                  (dotimes (j m)
  87                           (declare (fixnum j))
  88                           (setq x 0)
  89                           (dotimes (k p)
  90                                    (declare (fixnum k))
  91                                    (setq x (+ x
  92                                               (* (aref a i k) (aref b k j)))))
  93                           (setf (aref c i j) x)))
  94         (case rtype
  95           (matrix c)
  96           (number (aref c 0 0))
  97           (t (coerce (compound-data-seq c) rtype)))))))
  98
  99 (defun identity-matrix (n)
 100   "Return the identity matrix of rank N as a fixnum-valued array."
 101   (let ((result (make-array (list n n) :initial-element 0)))
 102     (dotimes (i n result)
 103              (declare (fixnum i))
 104              (setf (aref result i i) 1))))
 105
 106 (defun diagonal (x)
 107   "If X is a square matrix, returns the diagonal of X. If MxN,
 108 M>N, return the diagonal of the NxN upper-left sub-matrix.  If X is a
 109 sequence, returns a diagonal matrix of rank (length X) with diagonal
 110 elements eq to the elements of X.  Return type is a vector.  FIXME:
 111 could be made more efficient?"
 112   (cond ((matrixp x)
 113          (let* ((n (min (num-rows x) (num-cols x)))
 114                 (result (make-array n)))
 115            (dotimes (i n (coerce result 'list))
 116                     (setf (aref result i) (aref x i i)))))
 117         ((typep x 'sequence)
 118          (let* ((x (coerce x 'vector))
 119                 (n (length x))
 120                 (result (make-array (list n n) :initial-element 0)))
 121            (dotimes (i n result)
 122                     (setf (aref result i i) (aref x i)))))
 123         (t (error "argument must be a matrix or a sequence"))))
 124
 125 (defun row-list (x)
 126   "Returns a list of the rows of X as vectors"
 127   (lisp-array-matrix-p x)
 128   (let ((m (num-rows x))
 129         (n (num-cols x))
 130         (result nil))
 131     (declare (fixnum m n))
 132     (flet ((get-row (k)
 133              (declare (fixnum k))
 134              (let ((row (make-array n)))
 135                (dotimes (i n row)
 136                  (declare (fixnum i))
 137                  (setf (aref row i) (aref x k i))))))
 138       (dotimes (i m result)
 139         (declare (fixnum i))
 140         (setf result (cons (get-row (- m i 1)) result))))))
 141
 142 (defun column-list (x)
 143 "Args: (m)
 144 Returns a list of vectors, each vector represents a column of M."
 145   (lisp-array-matrix-p x)
 146   (let ((m (num-rows x))
 147         (n (num-cols x))
 148         (result nil))
 149     (declare (fixnum m n))
 150     (flet ((get-col (k)
 151              (declare (fixnum k))
 152              (let ((col (make-array m)))
 153                (dotimes (i m col)
 154                  (declare (fixnum i))
 155                  (setf (aref col i) (aref x i k))))))
 156       (dotimes (i n result)
 157         (declare (fixnum i))
 158         (setf result (cons (get-col (- n i 1)) result))))))
 159
 160 (defun inner-product (x y)
 161 "Args: (x y)
 162 Returns inner product of sequences X and Y."
 163   (check-type x sequence)
 164   (check-type y sequence)
 165   (assert (= (length x) (length y))
 166           (x y)
 167           "INNER-PRODUCT: sequences of different length")
 168   (let ((n (length x))
 169         (cx (make-next-element x))
 170         (cy (make-next-element y))
 171         (result 0))
 172     (declare (fixnum n))
 173     (dotimes (i n result)
 174       (declare (fixnum i))
 175       (setf result
 176             (+ result (* (get-next-element cx i) (get-next-element cy i)))))))
 177
 178 (defun outer-product (x y &optional (f #'*))
 179 "Args: (x y &optional (fcn #'*))
 180 Returns the generalized outer product of x and y, using fcn. Tat is, the result
 181 is a matrix of dimension ((length x) (length y)) and the (i j) element of the
 182 result is computed as (apply fcn (aref x i) (aref y j))."
 183   (let* ((x (coerce x 'vector))
 184          (y (coerce y 'vector))
 185          (m (length x))
 186          (n (length y))
 187          (a (make-array (list m n))))
 188     (declare (fixnum m n))
 189     (dotimes (i m a)
 190       (declare (fixnum i))
 191       (dotimes (j n)
 192         (declare (fixnum j))
 193         (setf (aref a i j) (funcall f (aref x i) (aref y j)))))))
 194
 195 (defun cross-product (x)
 196   "Args: (x)
 197 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
 198  (inner-product X X).  This function works with raw LISP ARRAYS."
 199   (lisp-array-matrix-p x)
 200   (let* ((n (num-rows x))
 201          (p (num-cols x))
 202          (c (make-array (list p p))))
 203     (declare (fixnum n p))
 204     (dotimes (i p c)
 205       (declare (fixnum i))
 206       (dotimes (j (+ i 1))
 207         (declare (fixnum j))
 208         (let ((val 0))
 209           (dotimes (k n)
 210             (declare (fixnum k))
 211             (incf val (* (aref x k i) (aref x k j))))
 212           (setf (aref c i j) val)
 213           (setf (aref c j i) val))))))
 214
 215 (defun transpose-list (x)
 216   (let ((m (length (first x))))
 217     (dolist (next x)
 218       (if (not (consp next)) (error "not a list - ~a" x))
 219       (if (/= m (length next)) (error "sublists not the same length")))
 220     (do* ((cx (copy-list x))
 221           (result (make-list m))
 222           (next result (cdr next)))
 223          ((null next) result)
 224       (setf (first next) (mapcar #'first cx))
 225       (do ((next cx (cdr next)))
 226           ((null next))
 227         (setf (first next) (rest (first next)))))))
 228
 229 (defun transpose (x)
 230 "Args: (m)
 231 Returns the transpose of the matrix M."
 232   (cond
 233    ((consp x) (transpose-list x))
 234    (t
 235     (lisp-array-matrix-p x)
 236     (let* ((m (num-rows x))
 237            (n (num-cols x))
 238            (tx (make-array (list n m))))
 239       (declare (fixnum m n))
 240       (dotimes (i m tx)
 241         (declare (fixnum i))
 242         (dotimes (j n)
 243           (declare (fixnum j))
 244           (setf (aref tx j i) (aref x i j))))))))
 245
 246 (defun bind-columns (&rest args)
 247 "Args (&rest args)
 248 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
 249 along their columns.
 250 Example: (bind-columns #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2 5)(3 4 6))"
 251   (flet ((check-arg (x)
 252            (if (not (or (typep x 'sequence) (matrixp x)))
 253                (error "bad argument type")))
 254          (arg-cols (x) (if (typep x 'sequence) 1 (num-cols x)))
 255          (arg-rows (x) (if (typep x 'sequence) (length x) (num-rows x))))
 256     (dolist (x args) (check-arg x))
 257     (let ((m (arg-rows (first args)))
 258           (n (arg-cols (first args))))
 259       (declare (fixnum m n))
 260       (dolist (x (rest args))
 261         (if (/= m (arg-rows x)) (error "column lengths do not match"))
 262         (incf n (arg-cols x)))
 263       (do* ((result (make-array (list m n)))
 264             (args args (rest args))
 265             (firstcol 0)
 266             (x (first args) (first args)))
 267            ((null args) result)
 268         (cond
 269          ((typep x 'sequence)
 270           (let ((cx (make-next-element x)))
 271             (dotimes (i m)
 272               (setf (aref result i firstcol) (get-next-element cx i)))))
 273          (t
 274           (let ((k (arg-cols x)))
 275             (dotimes (i m)
 276               (dotimes (j k)
 277                 (setf (aref result i (+ firstcol j)) (aref x i j)))))))
 278         (incf firstcol (arg-cols x))))))
 279
 280 (defun bind-rows (&rest args)
 281 "Args (&rest args)
 282 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
 283 along their rows.
 284 Example: (bind-rows #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2)(3 4)(5 6))"
 285   (flet ((check-arg (x)
 286            (if (not (or (typep x 'sequence) (matrixp x)))
 287                (error "bad argument type")))
 288          (arg-cols (x) (if (typep x 'sequence) (length x) (num-cols x)))
 289          (arg-rows (x) (if (typep x 'sequence) 1 (num-rows x))))
 290     (dolist (x args) (check-arg x))
 291     (let ((m (arg-rows (first args)))
 292           (n (arg-cols (first args))))
 293       (declare (fixnum m n))
 294       (dolist (x (rest args))
 295         (if (/= n (arg-cols x)) (error "row lengths do not match"))
 296         (incf m (arg-rows x)))
 297       (do* ((result (make-array (list m n)))
 298             (args args (rest args))
 299             (firstrow 0)
 300             (x (first args) (first args)))
 301            ((null args) result)
 302         (cond
 303          ((typep x 'sequence)
 304           (let ((cx (make-next-element x)))
 305             (dotimes (i n)
 306               (setf (aref result firstrow i) (get-next-element cx i)))))
 307          (t
 308           (let ((k (arg-rows x)))
 309             (dotimes (i n)
 310               (dotimes (j k)
 311                 (setf (aref result (+ firstrow j) i) (aref x j i)))))))
 312         (incf firstrow (arg-rows x))))))