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