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