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