matrix-matlisp.lisp

   1 i;;; -*- mode: lisp -*-
   2 ;;;
   3 ;;; Copyright (c) 2008--, 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 ;; matrix-matlisp -- extended functions for matrix and linear algebra
   9 ;; using MATLISP.
  10
  11 ;; Package Setup
  12
  13 (in-package :cl-user)
  14
  15 (defpackage :lisp-stat-matrix-matlisp
  16   (:use :common-lisp
  17         :matlisp)
  18   (:export sweep
  19            ))
  20
  21 (in-package :lisp-stat-matrix-matlisp)
  22
  23 ;; Might check with features that we've (pushnew :matlisp *features*)
  24 ;; and then we can do  #+matlisp    as needed.
  25
  26 ;;;
  27 ;;; SWEEP Operator, MATLISP implementation.
  28 ;;;
  29
  30 (defun make-sweep-front (x y w n p mode has_w x_mean result)
  31   "Compute the sweep matrix, a generalized computation based on the regression components Y and X.
  32 X is: numerical version of design/regt matrix
  33 y is: response
  34 w is: weights (might be simply 1
  35 n is: ?
  36 p is: ?
  37 mode is: real, complex, though complex is not yet supported.
  38 has_w just indicates if weights (w) are meaningful.
  39 x_mean:
  40 result: ? place to store?
  41 "
  42   (declare (fixnum n p mode has_w))
  43   (let ((x_data nil)
  44         (result_data nil)
  45         (val 0.0)
  46         (dxi 0.0)
  47         (dyi 0.0)
  48         (dv 0.0)
  49         (dw 0.0)
  50         (sum_w 0.0)
  51         (dxik 0.0)
  52         (dxjk 0.0)
  53         (dyj 0.0)
  54         (dx_meani 0.0)
  55         (dx_meanj 0.0)
  56         (dy_mean 0.0)
  57         (has-w (if (/= 0 has_w) t nil))
  58         (RE 1))
  59     (declare (long-float val dxi dyi dv dw sum_w dxik dxjk dyj
  60                     dx_meani dx_meanj dy_mean))
  61
  62     (setf x_data (compound-data-seq x))
  63     (setf result_data (compound-data-seq result))
  64
  65     ;; find the mean of y
  66     (dotimes (i n)
  67       (declare (fixnum i))
  68       (setf dyi (makedouble (aref y i)))
  69       (when has-w
  70         (setf dw (makedouble (aref w i)))
  71         (incf sum_w dw)
  72         (setf dyi (* dyi dw)))
  73       (incf val dyi))
  74     (if (not has-w) (setf sum_w (float n 0.0)))
  75     (if (<= sum_w 0.0) (error "non positive sum of weights"))
  76     (setf dy_mean (/ val sum_w))
  77
  78     ;; find the column means of X
  79     (dotimes (j p)
  80       (declare (fixnum j))
  81       (setf val 0.0)
  82       (dotimes (i n)
  83         (declare (fixnum i))
  84         (setf dxi (makedouble (aref x_data (+ (* p i) j))))
  85         (when has-w
  86           (setf dw (makedouble (aref w i)))
  87           (setf dxi (* dxi dw)))
  88         (incf val dxi))
  89       (setf (aref x_mean j) (/ val sum_w)))
  90
  91     ;; put 1/sum_w in topleft, means on left, minus means on top
  92     (setf (aref result_data 0) (/ 1.0 sum_w))
  93     (dotimes (i p)
  94       (declare (fixnum i))
  95       (setf dxi (makedouble (aref x_mean i)))
  96       (setf (aref result_data (+ i 1)) (- dxi))
  97       (setf (aref result_data (* (+ i 1) (+ p 2))) dxi))
  98     (setf (aref result_data (+ p 1)) (- dy_mean))
  99     (setf (aref result_data (* (+ p 1) (+ p 2))) dy_mean)
 100
 101     ;; put sums of adjusted cross products in body
 102     (dotimes (i p)
 103       (declare (fixnum i))
 104       (dotimes (j p)
 105         (declare (fixnum j))
 106         (setf val 0.0)
 107         (dotimes (k n)
 108           (declare (fixnum k))
 109           (setf dxik (makedouble (aref x_data (+ (* p k) i))))
 110           (setf dxjk (makedouble (aref x_data (+ (* p k) j))))
 111           (setf dx_meani (makedouble (aref x_mean i)))
 112           (setf dx_meanj (makedouble (aref x_mean j)))
 113           (setf dv (* (- dxik dx_meani) (- dxjk dx_meanj)))
 114           (when has-w
 115             (setf dw (makedouble (aref w k)))
 116             (setf dv (* dv dw)))
 117           (incf val dv))
 118         (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ j 1))) val)
 119         (setf (aref result_data (+ (* (+ j 1) (+ p 2)) (+ i 1))) val))
 120       (setf val 0.0)
 121       (dotimes (j n)
 122         (declare (fixnum j))
 123         (setf dxik (makedouble (aref x_data (+ (* p j) i))))
 124         (setf dyj (makedouble (aref y j)))
 125         (setf dx_meani (makedouble (aref x_mean i)))
 126         (setf dv (* (- dxik dx_meani) (- dyj dy_mean)))
 127         (when has-w
 128           (setf dw (makedouble (aref w j)))
 129           (setf dv (* dv dw)))
 130         (incf val dv))
 131       (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ p 1))) val)
 132       (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ i 1))) val))
 133     (setf val 0.0)
 134     (dotimes (j n)
 135       (declare (fixnum j))
 136       (setf dyj (makedouble (aref y j)))
 137       (setf dv (* (- dyj dy_mean) (- dyj dy_mean)))
 138       (when has-w
 139         (setf dw (makedouble (aref w j)))
 140         (setf dv (* dv dw)))
 141       (incf val dv))
 142     (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ p 1))) val)))
 143
 144 ;;; FIXME: use matlisp
 145 (defun sweep-in-place-front (a rows cols mode k tol)
 146   "Sweep algorithm for linear regression."
 147   (declare (long-float tol))
 148   (declare (fixnum rows cols mode k))
 149   (let ((data nil)
 150         (pivot 0.0)
 151         (aij 0.0)
 152         (aik 0.0)
 153         (akj 0.0)
 154         (akk 0.0)
 155         (RE 1))
 156     (declare (long-float pivot aij aik akj akk))
 157
 158     (if (> mode RE) (error "not supported for complex data yet"))
 159     (if (or (< k 0) (>= k rows) (>= k cols)) (error "index out of range"))
 160
 161     (setf tol (max tol machine-epsilon))
 162     (setf data (compound-data-seq a))
 163
 164     (setf pivot (makedouble (aref data (+ (* cols k) k))))
 165
 166     (cond
 167      ((or (> pivot tol) (< pivot (- tol)))
 168       (dotimes (i rows)
 169         (declare (fixnum i))
 170         (dotimes (j cols)
 171           (declare (fixnum j))
 172           (when (and (/= i k) (/= j k))
 173             (setf aij (makedouble (aref data (+ (* cols i) j))))
 174             (setf aik (makedouble (aref data (+ (* cols i) k))))
 175             (setf akj (makedouble (aref data (+ (* cols k) j))))
 176             (setf aij (- aij (/ (* aik akj) pivot)))
 177             (setf (aref data (+ (* cols i) j)) aij))))
 178
 179       (dotimes (i rows)
 180         (declare (fixnum i))
 181         (setf aik (makedouble (aref data (+ (* cols i) k))))
 182         (when (/= i k)
 183           (setf aik (/ aik pivot))
 184           (setf (aref data (+ (* cols i) k)) aik)))
 185
 186       (dotimes (j cols)
 187         (declare (fixnum j))
 188         (setf akj (makedouble (aref data (+ (* cols k) j))))
 189         (when (/= j k)
 190           (setf akj (- (/ akj pivot)))
 191           (setf (aref data (+ (* cols k) j)) akj)))
 192
 193       (setf akk (/ 1.0 pivot))
 194       (setf (aref data (+ (* cols k) k)) akk)
 195       1)
 196      (t 0))))
 197
 198 ;; FIXME: use matlisp
 199 (defun make-sweep-matrix (x y &optional w)
 200   "Args: (x y &optional weights)
 201  X is matrix, Y and WEIGHTS are sequences. Returns the sweep matrix of the
 202  (weighted) regression of Y on X"
 203    (if (not (typep x 'real-matrix)) (error "Make Sweep Matrix: X not real-matrix"))
 204    (if (not (typep y 'real-matrix)) (error "Make Sweep Matrix: Y not real-matrix"))
 205    (if w (if (not (typep w 'real-matrix)) (error "Make Sweep Matrix: W not real-matrix")))
 206    (let ((n (n x))
 207          (p (m x)))
 208      (if (/= n (length y)) (error "dimensions do not match"))
 209      (if (and w (/= n (length w))) (error "dimensions do not match"))
 210      (let ((mode (max (la-data-mode x)
 211                       (la-data-mode x)
 212                       (if w (la-data-mode w) 0)))
 213            (result (make-array (list (+ p 2) (+ p 2))))
 214            (x-mean (make-array p))
 215            (y (coerce y 'vector))
 216            (w (if w (coerce w 'vector)))
 217            (has-w (if w 1 0)))
 218        (make-sweep-front x y w n p mode has-w x-mean result)
 219        result)))
 220
 221  (defun sweep-in-place (a k tol)
 222   (check-matrix a)
 223   (check-one-fixnum k)
 224   (check-one-real tol)
 225   (let ((rows (num-rows a))
 226         (cols (num-cols a))
 227         (mode (la-data-mode a)))
 228     (let ((swept (sweep-in-place-front a rows cols mode k tol)))
 229       (if (/= 0 swept) t nil))))
 230
 231 (defun sweep-operator (a columns &optional tolerances)
 232   "Args: (a indices &optional tolerances)
 233
 234 A is a matrix, INDICES a sequence of the column indices to be
 235 swept. Returns a list of the swept result and the list of the columns
 236 actually swept. (See MULTREG documentation.) If supplied, TOLERANCES
 237 should be a list of real numbers the same length as INDICES. An index
 238 will only be swept if its pivot element is larger than the
 239 corresponding element of TOLERANCES."
 240
 241   ;; Why we should use generics/methods....
 242   (if (not (typep a 'real-matrix))
 243       (error "a is not a matrix"))
 244   ;; following: we need integer?
 245   (if (and (not (typep columns 'sequence))
 246            (not (typep (setf columns (list columns)) 'sequence)))
 247       (error "columns not coerceable to sequence"))
 248   (if tolerances
 249       (if (and (not (typep tolerances 'sequence))
 250                (not (setf tolerances (list tolerances))))
 251           (error "tolerances not coercable to sequence.")))
 252
 253   ;; (check-real a) done by matlisp.
 254   (check-fixnum columns)
 255   (if tolerances (check-real tolerances))
 256
 257   (do ((tol .0000001)
 258        (result (copy-array a))
 259        (swept-columns nil)
 260        (columns (coerce columns 'list) (cdr columns))
 261        (tolerances (if (consp tolerances) (coerce tolerances 'list))
 262                    (if (consp tolerances) (cdr tolerances))))
 263       ((null columns) (list result swept-columns))
 264     ;; this should be done in the context of matlisp, I think...?
 265     (let ((col (first columns))
 266             (tol (if (consp tolerances) (first tolerances) tol)))
 267         (if (sweep-in-place result col tol)
 268             (setf swept-columns (cons col swept-columns))))))
 269