src/data/matrix.lisp

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