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