3 ;;; Copyright (c) 2007, by A.J. Rossini <blindglobe@gmail.com>
4 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
6 ;;;; matrices for statistics. Extends CLEM
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.
17 (defpackage :lisp-stat-matrix
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
))
25 (in-package :lisp-stat-matrix
)
29 CL-USER
> (defparameter tr1
(make-instance 'number-matrix
))
32 #<NUMBER-MATRIX of dimensions
(1)>
33 CL-USER
> (defparameter tr2
(make-instance 'number-matrix
:rows
4 :cols
3))
36 #<NUMBER-MATRIX
[.000000000 .000000000 .000000000;
37 .000000000 .000000000 .000000000;
38 .000000000 .000000000 .000000000;
39 .000000000 .000000000 .000000000]>
40 CL-USER
> (transpose tr2
)
48 (defparameter tr3
(make-instance 'clem
::base-vector
:length
4))
52 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
54 ;;;; Array to Row-Major Data Vector Conversion Functions
56 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
58 (defun array-data-vector (a)
60 Displaces array A to a vector"
61 (make-array (array-total-size a
)
63 :element-type
(array-element-type a
)))
65 (defun vector-to-array (v dims
)
67 Displaces vector V to array with dimensions DIMS"
68 ;;; Yes, but using row or column first approach?
71 :element-type
(array-element-type v
)))
75 (defun check-matrix (a)
76 (if (not (and (arrayp a
) (= (array-rank a
) 2)))
77 (error "not a matrix - ~s" a
)
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
))
88 Returns T if X is a matrix, NIL otherwise."
89 (and (arrayp x
) (= (array-rank x
) 2)))
93 Returns number of rows in X."
94 (array-dimension x
0))
98 Returns number of columns in X."
99 (array-dimension x
1))
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
)))))
111 (setf a
(vector-to-array (coerce a
'vector
) (list 1 (length a
)))))
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"))
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
)))
123 (declare (fixnum n m p
))
132 (* (aref a i k
) (aref b k j
)))))
133 (setf (aref c i j
) x
)))
136 (number (aref c
0 0))
137 (t (coerce (compound-data-seq c
) rtype
)))))))
139 (defun identity-matrix (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
)
145 (setf (aref result i i
) 1))))
147 ;; this thing is not very efficient at this point - too much coercing
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
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
)))))
159 (let* ((x (coerce x
'vector
))
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"))))
168 Returns a list of the rows of M as vectors"
170 (let ((m (num-rows x
))
173 (declare (fixnum m n
))
176 (let ((row (make-array n
)))
179 (setf (aref row i
) (aref x k i
))))))
180 (dotimes (i m result
)
182 (setf result
(cons (get-row (- m i
1)) result
))))))
184 (defun column-list (x)
186 Returns a list of the columns of M as vectors"
188 (let ((m (num-rows x
))
191 (declare (fixnum m n
))
194 (let ((col (make-array m
)))
197 (setf (aref col i
) (aref x i k
))))))
198 (dotimes (i n result
)
200 (setf result
(cons (get-col (- n i
1)) result
))))))
202 (defun inner-product (x y
)
204 Returns inner product of sequences X and Y."
208 (cx (make-next-element x
))
209 (cy (make-next-element y
))
212 (if (/= n
(length y
)) (error "sequence lengths do not match"))
213 (dotimes (i n result
)
216 (+ result
(* (get-next-element cx i
) (get-next-element cy i
)))))))
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
))
227 (a (make-array (list m n
))))
228 (declare (fixnum m n
))
233 (setf (aref a i j
) (funcall f
(aref x i
) (aref y j
)))))))
235 (defun cross-product (x)
237 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
238 (inner-product X X)."
240 (let* ((n (num-rows x
))
242 (c (make-array (list p p
))))
243 (declare (fixnum n p
))
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
))))))
255 (defun transpose-list (x)
256 (let ((m (length (first 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
)))
264 (setf (first next
) (mapcar #'first cx
))
265 (do ((next cx
(cdr next
)))
267 (setf (first next
) (rest (first next
)))))))
271 Returns the transpose of the matrix M."
273 ((consp x
) (transpose-list x
))
276 (let* ((m (num-rows x
))
278 (tx (make-array (list n m
))))
279 (declare (fixnum m n
))
284 (setf (aref tx j i
) (aref x i j
))))))))
286 (defun bind-columns (&rest args
)
288 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
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
))
306 (x (first args
) (first args
)))
310 (let ((cx (make-next-element x
)))
312 (setf (aref result i firstcol
) (get-next-element cx i
)))))
314 (let ((k (arg-cols x
)))
317 (setf (aref result i
(+ firstcol j
)) (aref x i j
)))))))
318 (incf firstcol
(arg-cols x
))))))
320 (defun bind-rows (&rest args
)
322 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
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
))
340 (x (first args
) (first args
)))
344 (let ((cx (make-next-element x
)))
346 (setf (aref result firstrow i
) (get-next-element cx i
)))))
348 (let ((k (arg-rows x
)))
351 (setf (aref result
(+ firstrow j
) i
) (aref x j i
)))))))
352 (incf firstrow
(arg-rows x
))))))