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