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