matrices.lsp

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