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 (typep a 'sequence) (typep b 'sequence)) 'number)
 107                      ((typep a 'sequence) (if (consp a) 'list 'vector))
 108                      ((typep b 'sequence) (if (consp b) 'list 'vector)))))
 109
 110     (if (typep a 'sequence)
 111       (setf a (vector-to-array (coerce a 'vector) (list 1 (length a)))))
 112     (if (typep b 'sequence)
 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         ((typep x 'sequence)
 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 (typep x 'sequence) (matrixp x)))
 293                (error "bad argument type")))
 294          (arg-cols (x) (if (typep x 'sequence) 1 (num-cols x)))
 295          (arg-rows (x) (if (typep x 'sequence) (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 (typep x 'sequence) (matrixp x)))
 327                (error "bad argument type")))
 328          (arg-cols (x) (if (typep x 'sequence) (length x) (num-cols x)))
 329          (arg-rows (x) (if (typep x 'sequence) 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          ((typep x 'sequence)
 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))))))