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