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