src/numerics/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 ;;; Issues:
  14 ;;; #1 - Need to extend to use lisp-matrix?
  15 ;;; #2 - do as a second alternative?
  16
  17 (in-package :lisp-stat-matrix)
  18
  19 (deftype matrix () 'array)  ;; temp fix
  20
  21 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  22 ;;;;
  23 ;;;;          Array to Row-Major Data Vector Conversion Functions
  24 ;;;;
  25 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  26
  27 (defun array-data-vector (a)
  28 "Args: (a)
  29 Displaces array A to a vector"
  30   (make-array (array-total-size a)
  31               :displaced-to a
  32               :element-type (array-element-type a)))
  33
  34 (defun vector-to-array (v dims)
  35 "Args: (v dims)
  36 Displaces vector V to array with dimensions DIMS"
  37   (make-array dims
  38               :displaced-to v
  39               :element-type (array-element-type v)))
  40
  41 ;;;;
  42
  43 (defun check-matrix (a)
  44   (if (not (and (typep a' array)
  45                 (= (array-rank a) 2)))
  46       (error "not a matrix - ~s" a)
  47       t))
  48
  49 (defun check-square-matrix (a)
  50   (if (and (check-matrix a)
  51            (/= (array-dimension a 0) (array-dimension a 1))
  52            (error "matrix not square - ~s" a))
  53       t))
  54
  55 (defun matrixp (x)
  56 "Args: (x)
  57
  58 Returns:
  59    T if X is a 2-d array (i.e. a matrix),
  60    NIL otherwise."
  61   (and (typep x 'array)
  62        (= (array-rank x) 2)))
  63
  64 (defun num-rows (x)
  65   "Args: (x)
  66
  67 Returns number of rows in X."
  68   (if (matrixp x)
  69       (array-dimension x 0)
  70       (error "only useful for matrices.")))
  71
  72 (defun num-cols (x)
  73   "Args: (x)
  74
  75 Returns number of columns in X."
  76   (if (matrixp x)
  77       (array-dimension x 1)
  78       (error "only useful for matrices.")))
  79
  80
  81 ;;; Look at this!  Prime target for generic function dispatch!
  82 (defun matmult (a b &rest args)
  83   "Args: (a b &rest args)
  84 Returns the matrix product of matrices a, b, etc. If a is a vector it is
  85 treated as a row vector; if b is a vector it is treated as a column vector."
  86   ;; fixme: why does SBCL claim this is unreachable?
  87   (let ((rtype (cond ((and (typep a 'matrix)
  88                            (typep b 'matrix)) 'matrix)
  89                      ((and (typep a 'matrix)
  90                            (typep b 'sequence)) 'vector)
  91                      ((and (typep a 'sequence)
  92                            (typep b 'matrix)) 'vector)
  93                      ((and (typep a 'sequence)
  94                            (typep b 'sequence)) 'number)
  95                      ((typep a 'sequence)
  96                       (if (consp a) 'list 'vector))
  97                      ((typep b 'sequence)
  98                       (if (consp b) 'list 'vector)))))
  99
 100     (if (typep a 'sequence)
 101       (setf a (vector-to-array (coerce a 'vector) (list 1 (length a)))))
 102     (if (typep b 'sequence)
 103       (setf b (vector-to-array (coerce b 'vector) (list (length b) 1))))
 104     (if (not (= (array-dimension a 1) (array-dimension b 0)))
 105       (error "dimensions do not match"))
 106     (if args
 107       (reduce #'matmult args :initial-value (matmult a b))
 108       (let* ((n (array-dimension a 0))
 109              (m (array-dimension b 1))
 110              (p (array-dimension a 1))
 111              (c (make-array (list n m)))
 112              x)
 113         (declare (fixnum n m p))
 114         (dotimes (i n)
 115           (declare (fixnum i))
 116           (dotimes (j m)
 117             (declare (fixnum j))
 118             (setq x 0)
 119             (dotimes (k p)
 120               (declare (fixnum k))
 121               (setq x (+ x
 122                          (* (aref a i k)
 123                             (aref b k j)))))
 124             (setf (aref c i j) x)))
 125         (case rtype
 126           (matrix c)
 127           (number (aref c 0 0))
 128           (t (coerce (compound-data-seq c) rtype)))))))
 129
 130 (defun identity-matrix (n)
 131 "Args: (n)
 132 Returns the identity matrix of rank N."
 133   (let ((result (make-array (list n n) :initial-element 0)))
 134     (dotimes (i n result)
 135              (declare (fixnum i))
 136              (setf (aref result i i) 1))))
 137
 138 ;; this thing is not very efficient at this point - too much coercing
 139 (defun diagonal (x)
 140 "Args: (x)
 141 If X is a matrix, returns the diagonal of X. If X is a sequence, returns a
 142 diagonal matrix of rank (length X) with diagonal elements eq to the elements
 143 of X."
 144   (cond ((typep x 'matrix)
 145          (let* ((n (min (num-rows x) (num-cols x)))
 146                 (result (make-array n)))
 147            (dotimes (i n (coerce result 'list))
 148                     (setf (aref result i) (aref x i i)))))
 149         ((typep x 'sequence)
 150          (let* ((x (coerce x 'vector))
 151                 (n (length x))
 152                 (result (make-array (list n n) :initial-element 0)))
 153            (dotimes (i n result)
 154                     (setf (aref result i i) (aref x i)))))
 155         (t (error "argument must be a matrix or a sequence"))))
 156
 157 (defun row-list (x)
 158   "Args: (m)
 159
 160 Returns a list of the rows of M as vectors"
 161   (check-matrix x)
 162   (let ((m (num-rows x))
 163         (n (num-cols x))
 164         (result nil))
 165     (declare (fixnum m n))
 166     (flet ((get-row (k)
 167              (declare (fixnum k))
 168              (let ((row (make-array n)))
 169                (dotimes (i n row)
 170                  (declare (fixnum i))
 171                  (setf (aref row i) (aref x k i))))))
 172       (dotimes (i m result)
 173         (declare (fixnum i))
 174         (setf result (cons (get-row (- m i 1)) result))))))
 175
 176 (defun column-list (x)
 177   "Args: (m)
 178
 179 Returns a list of the columns of M as vectors"
 180   (check-matrix x)
 181   (let ((m (num-rows x))
 182         (n (num-cols x))
 183         (result nil))
 184     (declare (fixnum m n))
 185     (flet ((get-col (k)
 186              (declare (fixnum k))
 187              (let ((col (make-array m)))
 188                (dotimes (i m col)
 189                  (declare (fixnum i))
 190                  (setf (aref col i) (aref x i k))))))
 191       (dotimes (i n result)
 192         (declare (fixnum i))
 193         (setf result (cons (get-col (- n i 1)) result))))))
 194
 195 (defun inner-product (x y)
 196   "Args: (x y)
 197
 198 Returns inner product of sequences X and Y."
 199   (check-sequence x)
 200   (check-sequence y)
 201   (let ((n (length x))
 202         (cx (make-next-element x))
 203         (cy (make-next-element y))
 204         (result 0))
 205     (declare (fixnum n))
 206     (if (/= n (length y)) (error "sequence lengths do not match"))
 207     (dotimes (i n result)
 208       (declare (fixnum i))
 209       (setf result
 210             (+ result (* (get-next-element cx i)
 211                          (get-next-element cy i)))))))
 212
 213 (defun outer-product (x y &optional (f #'*))
 214   "Args: (x y &optional (fcn #'*))
 215
 216 Returns the generalized outer product of x and y, using fcn. Tat is, the result
 217 is a matrix of dimension ((length x) (length y)) and the (i j) element of the
 218 result is computed as (apply fcn (aref x i) (aref y j))."
 219
 220   (let* ((x (coerce x 'vector))
 221          (y (coerce y 'vector))
 222          (m (length x))
 223          (n (length y))
 224          (a (make-array (list m n))))
 225     (declare (fixnum m n))
 226     (dotimes (i m a)
 227       (declare (fixnum i))
 228       (dotimes (j n)
 229         (declare (fixnum j))
 230         (setf (aref a i j) (funcall f (aref x i) (aref y j)))))))
 231
 232 (defun cross-product (x)
 233   "Args: (x)
 234
 235 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
 236 (inner-product X X)."
 237
 238   (check-matrix x)
 239   (let* ((n (num-rows x))
 240          (p (num-cols x))
 241          (c (make-array (list p p))))
 242     (declare (fixnum n p))
 243     (dotimes (i p c)
 244       (declare (fixnum i))
 245       (dotimes (j (+ i 1))
 246         (declare (fixnum j))
 247         (let ((val 0))
 248           (dotimes (k n)
 249             (declare (fixnum k))
 250             (incf val (* (aref x k i) (aref x k j))))
 251           (setf (aref c i j) val)
 252           (setf (aref c j i) val))))))
 253
 254 (defun transpose-list (x)
 255   (let ((m (length (first x))))
 256     (dolist (next x)
 257       (if (not (consp next)) (error "not a list - ~a" x))
 258       (if (/= m (length next)) (error "sublists not the same length")))
 259     (do* ((cx (copy-list x))
 260           (result (make-list m))
 261           (next result (cdr next)))
 262          ((null next) result)
 263       (setf (first next) (mapcar #'first cx))
 264       (do ((next cx (cdr next)))
 265           ((null next))
 266         (setf (first next) (rest (first next)))))))
 267
 268 (defun transpose (x)
 269 "Args: (m)
 270 Returns the transpose of the matrix M."
 271   (cond
 272    ((consp x) (transpose-list x))
 273    (t
 274     (check-matrix x)
 275     (let* ((m (num-rows x))
 276            (n (num-cols x))
 277            (tx (make-array (list n m))))
 278       (declare (fixnum m n))
 279       (dotimes (i m tx)
 280         (declare (fixnum i))
 281         (dotimes (j n)
 282           (declare (fixnum j))
 283           (setf (aref tx j i) (aref x i j))))))))
 284
 285 (defun bind-columns (&rest args)
 286   "Args (&rest args)
 287
 288 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
 289 along their columns.  Example:
 290     (bind-columns #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2 5)(3 4 6))"
 291
 292   (flet ((check-arg (x)
 293            (if (not (or (typep x 'sequence) (typep x 'matrix)))
 294                (error "bad argument type")))
 295          (arg-cols (x) (if (typep x 'sequence) 1 (num-cols x)))
 296          (arg-rows (x) (if (typep x 'sequence) (length x) (num-rows x))))
 297     (dolist (x args) (check-arg x)) ;; verify data structure conformance.
 298     (let ((m (arg-rows (first args)))
 299           (n (arg-cols (first args))))
 300       (declare (fixnum m n))
 301       (dolist (x (rest args))
 302         (if (/= m (arg-rows x)) (error "column lengths do not match"))
 303         (incf n (arg-cols x)))
 304       (do* ((result (make-array (list m n)))
 305             (args args (rest args))
 306             (firstcol 0)
 307             (x (first args) (first args)))
 308            ((null args) result)
 309         (cond
 310           ((typep x 'sequence)
 311            (let ((cx (make-next-element x)))
 312             (dotimes (i m)
 313               (setf (aref result i firstcol) (get-next-element cx i)))))
 314           (t
 315            (let ((k (arg-cols x)))
 316              (dotimes (i m)
 317                (dotimes (j k)
 318                  (setf (aref result i (+ firstcol j)) (aref x i j)))))))
 319         (incf firstcol (arg-cols x))))))
 320
 321 (defun bind-rows (&rest args)
 322   "Args (&rest args)
 323
 324 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
 325 along their rows. Example:
 326    (bind-rows #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2)(3 4)(5 6))"
 327
 328   (flet ((check-arg (x)
 329            (if (not (or (typep x 'sequence)
 330                         (typep x 'matrix)))
 331                (error "bad argument type")))
 332          (arg-cols (x) (if (typep x 'sequence) (length x) (num-cols x)))
 333          (arg-rows (x) (if (typep x 'sequence) 1 (num-rows x))))
 334     (dolist (x args) (check-arg x))
 335     (let ((m (arg-rows (first args)))
 336           (n (arg-cols (first args))))
 337       (declare (fixnum m n))
 338       (dolist (x (rest args))
 339         (if (/= n (arg-cols x)) (error "row lengths do not match"))
 340         (incf m (arg-rows x)))
 341       (do* ((result (make-array (list m n)))
 342             (args args (rest args))
 343             (firstrow 0)
 344             (x (first args) (first args)))
 345            ((null args) result)
 346         (cond
 347          ((typep x 'sequence)
 348           (let ((cx (make-next-element x)))
 349             (dotimes (i n)
 350               (setf (aref result firstrow i) (get-next-element cx i)))))
 351          (t
 352           (let ((k (arg-rows x)))
 353             (dotimes (i n)
 354               (dotimes (j k)
 355                 (setf (aref result (+ firstrow j) i) (aref x j i)))))))
 356         (incf firstrow (arg-rows x))))))
 357
 358
 359 ;;;
 360 ;;; Copying Functions
 361 ;;;
 362
 363 (defun copy-vector (x)
 364   "Args: (x)
 365
 366 Returns a copy of the vector X"
 367   (copy-seq x))
 368
 369 (defun copy-array (a)
 370   "Args: (a)
 371
 372 Returns a copy of the array A"
 373   (vector-to-array (copy-seq (array-data-vector a))
 374                    (array-dimensions a)))
 375
 376
 377 #|
 378 (defgeneric copy (x)
 379    (:documentation "methods for copying linaar algebra forms."))
 380
 381 (defmethod copy ((x vector)))
 382
 383 (defmethod copy ((x matrix)))
 384 |#