more tests to track down multiple regression failures
[CommonLispStat.git] / matrices.lsp
blobdd318bb44bad3557de7baa869135fd146219ebb0
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.
8 ;;;; matrices -- Basic matrix operations
9 ;;;;
10 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
11 ;;;; unrestricted use.
14 ;;; Need to extend to use CLEM
16 ;;;
17 ;;; Package Setup
18 ;;;
20 (in-package :cl-user)
22 (defpackage :lisp-stat-matrix
23 (:use :common-lisp
24 :cffi
25 :lisp-stat-compound-data)
26 (:export matrixp num-rows num-cols matmult identity-matrix diagonal
27 row-list column-list inner-product outer-product
28 cross-product transpose bind-columns bind-rows
29 array-data-vector vector-to-array
31 check-matrix check-square-matrix
33 copy-array copy-vector
36 (in-package :lisp-stat-matrix)
38 (deftype matrix () 'array) ;; temp fix
40 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
41 ;;;;
42 ;;;; Array to Row-Major Data Vector Conversion Functions
43 ;;;;
44 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
46 (defun array-data-vector (a)
47 "Args: (a)
48 Displaces array A to a vector"
49 (make-array (array-total-size a)
50 :displaced-to a
51 :element-type (array-element-type a)))
53 (defun vector-to-array (v dims)
54 "Args: (v dims)
55 Displaces vector V to array with dimensions DIMS"
56 (make-array dims
57 :displaced-to v
58 :element-type (array-element-type v)))
60 ;;;;
62 (defun check-matrix (a)
63 (if (not (and (typep a' array)
64 (= (array-rank a) 2)))
65 (error "not a matrix - ~s" a)
66 t))
68 (defun check-square-matrix (a)
69 (if (and (check-matrix a)
70 (/= (array-dimension a 0) (array-dimension a 1))
71 (error "matrix not square - ~s" a))
72 t))
74 (defun matrixp (x)
75 "Args: (x)
76 Returns T if X is a matrix, NIL otherwise."
77 (and (typep x 'array)
78 (= (array-rank x) 2)))
80 (defun num-rows (x)
81 "Args: (x)
82 Returns number of rows in X."
83 (array-dimension x 0))
85 (defun num-cols (x)
86 "Args: (x)
87 Returns number of columns in X."
88 (array-dimension x 1))
91 ;;; Look at this! Prime target for generic function dispatch!
92 (defun matmult (a b &rest args)
93 "Args: (a b &rest args)
94 Returns the matrix product of matrices a, b, etc. If a is a vector it is
95 treated as a row vector; if b is a vector it is treated as a column vector."
96 ;; fixme: why does SBCL claim this is unreachable?
97 (let ((rtype (cond ((and (typep a 'matrix)
98 (typep b 'matrix)) 'matrix)
99 ((and (typep a 'matrix)
100 (typep b 'sequence)) 'vector)
101 ((and (typep a 'sequence)
102 (typep b 'matrix)) 'vector)
103 ((and (typep a 'sequence)
104 (typep b 'sequence)) 'number)
105 ((typep a 'sequence)
106 (if (consp a) 'list 'vector))
107 ((typep b 'sequence)
108 (if (consp b) 'list 'vector)))))
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)))
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)
133 (aref b k j)))))
134 (setf (aref c i j) x)))
135 (case rtype
136 (matrix c)
137 (number (aref c 0 0))
138 (t (coerce (compound-data-seq c) rtype)))))))
140 (defun identity-matrix (n)
141 "Args: (n)
142 Returns the identity matrix of rank N."
143 (let ((result (make-array (list n n) :initial-element 0)))
144 (dotimes (i n result)
145 (declare (fixnum i))
146 (setf (aref result i i) 1))))
148 ;; this thing is not very efficient at this point - too much coercing
149 (defun diagonal (x)
150 "Args: (x)
151 If X is a matrix, returns the diagonal of X. If X is a sequence, returns a
152 diagonal matrix of rank (length X) with diagonal elements eq to the elements
153 of X."
154 (cond ((typep x 'matrix)
155 (let* ((n (min (num-rows x) (num-cols x)))
156 (result (make-array n)))
157 (dotimes (i n (coerce result 'list))
158 (setf (aref result i) (aref x i i)))))
159 ((typep x 'sequence)
160 (let* ((x (coerce x 'vector))
161 (n (length x))
162 (result (make-array (list n n) :initial-element 0)))
163 (dotimes (i n result)
164 (setf (aref result i i) (aref x i)))))
165 (t (error "argument must be a matrix or a sequence"))))
167 (defun row-list (x)
168 "Args: (m)
169 Returns a list of the rows of M as vectors"
170 (check-matrix x)
171 (let ((m (num-rows x))
172 (n (num-cols x))
173 (result nil))
174 (declare (fixnum m n))
175 (flet ((get-row (k)
176 (declare (fixnum k))
177 (let ((row (make-array n)))
178 (dotimes (i n row)
179 (declare (fixnum i))
180 (setf (aref row i) (aref x k i))))))
181 (dotimes (i m result)
182 (declare (fixnum i))
183 (setf result (cons (get-row (- m i 1)) result))))))
185 (defun column-list (x)
186 "Args: (m)
187 Returns a list of the columns of M as vectors"
188 (check-matrix x)
189 (let ((m (num-rows x))
190 (n (num-cols x))
191 (result nil))
192 (declare (fixnum m n))
193 (flet ((get-col (k)
194 (declare (fixnum k))
195 (let ((col (make-array m)))
196 (dotimes (i m col)
197 (declare (fixnum i))
198 (setf (aref col i) (aref x i k))))))
199 (dotimes (i n result)
200 (declare (fixnum i))
201 (setf result (cons (get-col (- n i 1)) result))))))
203 (defun inner-product (x y)
204 "Args: (x y)
205 Returns inner product of sequences X and Y."
206 (check-sequence x)
207 (check-sequence y)
208 (let ((n (length x))
209 (cx (make-next-element x))
210 (cy (make-next-element y))
211 (result 0))
212 (declare (fixnum n))
213 (if (/= n (length y)) (error "sequence lengths do not match"))
214 (dotimes (i n result)
215 (declare (fixnum i))
216 (setf result
217 (+ result (* (get-next-element cx i) (get-next-element cy i)))))))
219 (defun outer-product (x y &optional (f #'*))
220 "Args: (x y &optional (fcn #'*))
221 Returns the generalized outer product of x and y, using fcn. Tat is, the result
222 is a matrix of dimension ((length x) (length y)) and the (i j) element of the
223 result is computed as (apply fcn (aref x i) (aref y j))."
224 (let* ((x (coerce x 'vector))
225 (y (coerce y 'vector))
226 (m (length x))
227 (n (length y))
228 (a (make-array (list m n))))
229 (declare (fixnum m n))
230 (dotimes (i m a)
231 (declare (fixnum i))
232 (dotimes (j n)
233 (declare (fixnum j))
234 (setf (aref a i j) (funcall f (aref x i) (aref y j)))))))
236 (defun cross-product (x)
237 "Args: (x)
238 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
239 (inner-product X X)."
240 (check-matrix x)
241 (let* ((n (num-rows x))
242 (p (num-cols x))
243 (c (make-array (list p p))))
244 (declare (fixnum n p))
245 (dotimes (i p c)
246 (declare (fixnum i))
247 (dotimes (j (+ i 1))
248 (declare (fixnum j))
249 (let ((val 0))
250 (dotimes (k n)
251 (declare (fixnum k))
252 (incf val (* (aref x k i) (aref x k j))))
253 (setf (aref c i j) val)
254 (setf (aref c j i) val))))))
256 (defun transpose-list (x)
257 (let ((m (length (first x))))
258 (dolist (next x)
259 (if (not (consp next)) (error "not a list - ~a" x))
260 (if (/= m (length next)) (error "sublists not the same length")))
261 (do* ((cx (copy-list x))
262 (result (make-list m))
263 (next result (cdr next)))
264 ((null next) result)
265 (setf (first next) (mapcar #'first cx))
266 (do ((next cx (cdr next)))
267 ((null next))
268 (setf (first next) (rest (first next)))))))
270 (defun transpose (x)
271 "Args: (m)
272 Returns the transpose of the matrix M."
273 (cond
274 ((consp x) (transpose-list x))
276 (check-matrix x)
277 (let* ((m (num-rows x))
278 (n (num-cols x))
279 (tx (make-array (list n m))))
280 (declare (fixnum m n))
281 (dotimes (i m tx)
282 (declare (fixnum i))
283 (dotimes (j n)
284 (declare (fixnum j))
285 (setf (aref tx j i) (aref x i j))))))))
287 (defun bind-columns (&rest args)
288 "Args (&rest args)
289 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
290 along their columns.
291 Example: (bind-columns #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2 5)(3 4 6))"
292 (flet ((check-arg (x)
293 (if (not (or (typep x 'sequence) (typep x 'matrix)))
294 (error "bad argument type")))
295 (arg-cols (x) (if (typep x 'sequence) 1 (num-cols x)))
296 (arg-rows (x) (if (typep x 'sequence) (length x) (num-rows x))))
297 (dolist (x args) (check-arg x)) ;; verify data structure conformance.
298 (let ((m (arg-rows (first args)))
299 (n (arg-cols (first args))))
300 (declare (fixnum m n))
301 (dolist (x (rest args))
302 (if (/= m (arg-rows x)) (error "column lengths do not match"))
303 (incf n (arg-cols x)))
304 (do* ((result (make-array (list m n)))
305 (args args (rest args))
306 (firstcol 0)
307 (x (first args) (first args)))
308 ((null args) result)
309 (cond
310 ((typep x 'sequence)
311 (let ((cx (make-next-element x)))
312 (dotimes (i m)
313 (setf (aref result i firstcol) (get-next-element cx i)))))
315 (let ((k (arg-cols x)))
316 (dotimes (i m)
317 (dotimes (j k)
318 (setf (aref result i (+ firstcol j)) (aref x i j)))))))
319 (incf firstcol (arg-cols x))))))
321 (defun bind-rows (&rest args)
322 "Args (&rest args)
323 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
324 along their rows.
325 Example: (bind-rows #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2)(3 4)(5 6))"
326 (flet ((check-arg (x)
327 (if (not (or (typep x 'sequence)
328 (typep x 'matrix)))
329 (error "bad argument type")))
330 (arg-cols (x) (if (typep x 'sequence) (length x) (num-cols x)))
331 (arg-rows (x) (if (typep x 'sequence) 1 (num-rows x))))
332 (dolist (x args) (check-arg x))
333 (let ((m (arg-rows (first args)))
334 (n (arg-cols (first args))))
335 (declare (fixnum m n))
336 (dolist (x (rest args))
337 (if (/= n (arg-cols x)) (error "row lengths do not match"))
338 (incf m (arg-rows x)))
339 (do* ((result (make-array (list m n)))
340 (args args (rest args))
341 (firstrow 0)
342 (x (first args) (first args)))
343 ((null args) result)
344 (cond
345 ((typep x 'sequence)
346 (let ((cx (make-next-element x)))
347 (dotimes (i n)
348 (setf (aref result firstrow i) (get-next-element cx i)))))
350 (let ((k (arg-rows x)))
351 (dotimes (i n)
352 (dotimes (j k)
353 (setf (aref result (+ firstrow j) i) (aref x j i)))))))
354 (incf firstrow (arg-rows x))))))
355 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
356 ;;;;
357 ;;;; Copying Functions
358 ;;;;
359 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
362 ;;; COPY-VECTOR function
365 (defun copy-vector (x)
366 "Args: (x)
367 Returns a copy of the vector X"
368 (copy-seq x))
371 ;;; COPY-ARRAY function
374 (defun copy-array (a)
375 "Args: (a)
376 Returns a copy of the array A"
377 (vector-to-array (copy-seq (array-data-vector a))
378 (array-dimensions a)))