src/numerics/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 ;;; Issues:
  14 ;;; #1 - Need to extend to use lisp-matrix?
  15 ;;; #2 - do as a second alternative?
  16
  17 (in-package :lisp-stat-matrix)
  18
  19
  20 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  21 ;;;;
  22 ;;;;          Array to Row-Major Data Vector Conversion Functions
  23 ;;;;
  24 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  25
  26 (defun array-data-vector (a)
  27 "Args: (a)
  28 Displaces array A to a vector"
  29   (make-array (array-total-size a)
  30               :displaced-to a
  31               :element-type (array-element-type a)))
  32
  33 (defun vector-to-array (v dims)
  34 "Args: (v dims)
  35 Displaces vector V to array with dimensions DIMS"
  36   (make-array dims
  37               :displaced-to v
  38               :element-type (array-element-type v)))
  39
  40
  41 (defun check-square-matrix (a)
  42   (if (and (check-matrix a)
  43            (/= (array-dimension a 0) (array-dimension a 1))
  44            (error "matrix not square - ~s" a))
  45       t))
  46
  47 (defun inner-product (x y)
  48   "Args: (x y)
  49
  50 Returns inner product of sequences X and Y."
  51   (check-sequence x)
  52   (check-sequence y)
  53   (let ((n (length x))
  54         (cx (make-next-element x))
  55         (cy (make-next-element y))
  56         (result 0))
  57     (declare (fixnum n))
  58     (if (/= n (length y)) (error "sequence lengths do not match"))
  59     (dotimes (i n result)
  60       (declare (fixnum i))
  61       (setf result
  62             (+ result (* (get-next-element cx i)
  63                          (get-next-element cy i)))))))
  64
  65 (defun outer-product (x y &optional (f #'*))
  66   "Args: (x y &optional (fcn #'*))
  67
  68 Returns the generalized outer product of x and y, using fcn. Tat is, the result
  69 is a matrix of dimension ((length x) (length y)) and the (i j) element of the
  70 result is computed as (apply fcn (aref x i) (aref y j))."
  71
  72   (let* ((x (coerce x 'vector))
  73          (y (coerce y 'vector))
  74          (m (length x))
  75          (n (length y))
  76          (a (make-array (list m n))))
  77     (declare (fixnum m n))
  78     (dotimes (i m a)
  79       (declare (fixnum i))
  80       (dotimes (j n)
  81         (declare (fixnum j))
  82         (setf (aref a i j) (funcall f (aref x i) (aref y j)))))))
  83
  84 (defun cross-product (x)
  85   "Args: (x)
  86
  87 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
  88 (inner-product X X)."
  89
  90   (check-matrix x)
  91   (let* ((n (num-rows x))
  92          (p (num-cols x))
  93          (c (make-array (list p p))))
  94     (declare (fixnum n p))
  95     (dotimes (i p c)
  96       (declare (fixnum i))
  97       (dotimes (j (+ i 1))
  98         (declare (fixnum j))
  99         (let ((val 0))
 100           (dotimes (k n)
 101             (declare (fixnum k))
 102             (incf val (* (aref x k i) (aref x k j))))
 103           (setf (aref c i j) val)
 104           (setf (aref c j i) val))))))
 105
 106 (defun transpose-list (x)
 107   (let ((m (length (first x))))
 108     (dolist (next x)
 109       (if (not (consp next)) (error "not a list - ~a" x))
 110       (if (/= m (length next)) (error "sublists not the same length")))
 111     (do* ((cx (copy-list x))
 112           (result (make-list m))
 113           (next result (cdr next)))
 114          ((null next) result)
 115       (setf (first next) (mapcar #'first cx))
 116       (do ((next cx (cdr next)))
 117           ((null next))
 118         (setf (first next) (rest (first next)))))))
 119
 120 (defun transpose (x)
 121 "Args: (m)
 122 Returns the transpose of the matrix M."
 123   (cond
 124    ((consp x) (transpose-list x))
 125    (t
 126     (check-matrix x)
 127     (let* ((m (num-rows x))
 128            (n (num-cols x))
 129            (tx (make-array (list n m))))
 130       (declare (fixnum m n))
 131       (dotimes (i m tx)
 132         (declare (fixnum i))
 133         (dotimes (j n)
 134           (declare (fixnum j))
 135           (setf (aref tx j i) (aref x i j))))))))
 136
 137