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