defvar to fix var unknowns, declared away some obvious unused vars, but still need...
[CommonLispStat.git] / matrices.lsp
blob3e38490efcd50bc96ded64cd2cb2c4f639ed0dd0
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 (arrayp a) (= (array-rank a) 2)))
64 (error "not a matrix - ~s" a)
65 t))
67 (defun check-square-matrix (a)
68 (if (and (check-matrix a)
69 (/= (array-dimension a 0) (array-dimension a 1))
70 (error "matrix not square - ~s" a))
71 t))
73 (defun matrixp (x)
74 "Args: (x)
75 Returns T if X is a matrix, NIL otherwise."
76 (and (arrayp x) (= (array-rank x) 2)))
78 (defun num-rows (x)
79 "Args: (x)
80 Returns number of rows in X."
81 (array-dimension x 0))
83 (defun num-cols (x)
84 "Args: (x)
85 Returns number of columns in X."
86 (array-dimension x 1))
88 (defun matmult (a b &rest args)
89 "Args: (a b &rest args)
90 Returns the matrix product of matrices a, b, etc. If a is a vector it is
91 treated as a row vector; if b is a vector it is treated as a column vector."
92 ;; fixme: why does SBCL claim this is unreachable?
93 (let ((rtype (cond ((and (matrixp a) (matrixp b)) 'matrix)
94 ((and (sequencep a) (sequencep b)) 'number)
95 ((sequencep a) (if (consp a) 'list 'vector))
96 ((sequencep b) (if (consp b) 'list 'vector)))))
98 (if (sequencep a)
99 (setf a (vector-to-array (coerce a 'vector) (list 1 (length a)))))
100 (if (sequencep b)
101 (setf b (vector-to-array (coerce b 'vector) (list (length b) 1))))
102 (if (not (= (array-dimension a 1) (array-dimension b 0)))
103 (error "dimensions do not match"))
104 (if args
105 (reduce #'matmult args :initial-value (matmult a b))
106 (let* ((n (array-dimension a 0))
107 (m (array-dimension b 1))
108 (p (array-dimension a 1))
109 (c (make-array (list n m)))
111 (declare (fixnum n m p))
112 (dotimes (i n)
113 (declare (fixnum i))
114 (dotimes (j m)
115 (declare (fixnum j))
116 (setq x 0)
117 (dotimes (k p)
118 (declare (fixnum k))
119 (setq x (+ x
120 (* (aref a i k) (aref b k j)))))
121 (setf (aref c i j) x)))
122 (case rtype
123 (matrix c)
124 (number (aref c 0 0))
125 (t (coerce (compound-data-seq c) rtype)))))))
127 (defun identity-matrix (n)
128 "Args: (n)
129 Returns the identity matrix of rank N."
130 (let ((result (make-array (list n n) :initial-element 0)))
131 (dotimes (i n result)
132 (declare (fixnum i))
133 (setf (aref result i i) 1))))
135 ;; this thing is not very efficient at this point - too much coercing
136 (defun diagonal (x)
137 "Args: (x)
138 If X is a matrix, returns the diagonal of X. If X is a sequence, returns a
139 diagonal matrix of rank (length X) with diagonal elements eq to the elements
140 of X."
141 (cond ((matrixp x)
142 (let* ((n (min (num-rows x) (num-cols x)))
143 (result (make-array n)))
144 (dotimes (i n (coerce result 'list))
145 (setf (aref result i) (aref x i i)))))
146 ((sequencep x)
147 (let* ((x (coerce x 'vector))
148 (n (length x))
149 (result (make-array (list n n) :initial-element 0)))
150 (dotimes (i n result)
151 (setf (aref result i i) (aref x i)))))
152 (t (error "argument must be a matrix or a sequence"))))
154 (defun row-list (x)
155 "Args: (m)
156 Returns a list of the rows of M as vectors"
157 (check-matrix x)
158 (let ((m (num-rows x))
159 (n (num-cols x))
160 (result nil))
161 (declare (fixnum m n))
162 (flet ((get-row (k)
163 (declare (fixnum k))
164 (let ((row (make-array n)))
165 (dotimes (i n row)
166 (declare (fixnum i))
167 (setf (aref row i) (aref x k i))))))
168 (dotimes (i m result)
169 (declare (fixnum i))
170 (setf result (cons (get-row (- m i 1)) result))))))
172 (defun column-list (x)
173 "Args: (m)
174 Returns a list of the columns of M as vectors"
175 (check-matrix x)
176 (let ((m (num-rows x))
177 (n (num-cols x))
178 (result nil))
179 (declare (fixnum m n))
180 (flet ((get-col (k)
181 (declare (fixnum k))
182 (let ((col (make-array m)))
183 (dotimes (i m col)
184 (declare (fixnum i))
185 (setf (aref col i) (aref x i k))))))
186 (dotimes (i n result)
187 (declare (fixnum i))
188 (setf result (cons (get-col (- n i 1)) result))))))
190 (defun inner-product (x y)
191 "Args: (x y)
192 Returns inner product of sequences X and Y."
193 (check-sequence x)
194 (check-sequence y)
195 (let ((n (length x))
196 (cx (make-next-element x))
197 (cy (make-next-element y))
198 (result 0))
199 (declare (fixnum n))
200 (if (/= n (length y)) (error "sequence lengths do not match"))
201 (dotimes (i n result)
202 (declare (fixnum i))
203 (setf result
204 (+ result (* (get-next-element cx i) (get-next-element cy i)))))))
206 (defun outer-product (x y &optional (f #'*))
207 "Args: (x y &optional (fcn #'*))
208 Returns the generalized outer product of x and y, using fcn. Tat is, the result
209 is a matrix of dimension ((length x) (length y)) and the (i j) element of the
210 result is computed as (apply fcn (aref x i) (aref y j))."
211 (let* ((x (coerce x 'vector))
212 (y (coerce y 'vector))
213 (m (length x))
214 (n (length y))
215 (a (make-array (list m n))))
216 (declare (fixnum m n))
217 (dotimes (i m a)
218 (declare (fixnum i))
219 (dotimes (j n)
220 (declare (fixnum j))
221 (setf (aref a i j) (funcall f (aref x i) (aref y j)))))))
223 (defun cross-product (x)
224 "Args: (x)
225 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
226 (inner-product X X)."
227 (check-matrix x)
228 (let* ((n (num-rows x))
229 (p (num-cols x))
230 (c (make-array (list p p))))
231 (declare (fixnum n p))
232 (dotimes (i p c)
233 (declare (fixnum i))
234 (dotimes (j (+ i 1))
235 (declare (fixnum j))
236 (let ((val 0))
237 (dotimes (k n)
238 (declare (fixnum k))
239 (incf val (* (aref x k i) (aref x k j))))
240 (setf (aref c i j) val)
241 (setf (aref c j i) val))))))
243 (defun transpose-list (x)
244 (let ((m (length (first x))))
245 (dolist (next x)
246 (if (not (consp next)) (error "not a list - ~a" x))
247 (if (/= m (length next)) (error "sublists not the same length")))
248 (do* ((cx (copy-list x))
249 (result (make-list m))
250 (next result (cdr next)))
251 ((null next) result)
252 (setf (first next) (mapcar #'first cx))
253 (do ((next cx (cdr next)))
254 ((null next))
255 (setf (first next) (rest (first next)))))))
257 (defun transpose (x)
258 "Args: (m)
259 Returns the transpose of the matrix M."
260 (cond
261 ((consp x) (transpose-list x))
263 (check-matrix x)
264 (let* ((m (num-rows x))
265 (n (num-cols x))
266 (tx (make-array (list n m))))
267 (declare (fixnum m n))
268 (dotimes (i m tx)
269 (declare (fixnum i))
270 (dotimes (j n)
271 (declare (fixnum j))
272 (setf (aref tx j i) (aref x i j))))))))
274 (defun bind-columns (&rest args)
275 "Args (&rest args)
276 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
277 along their columns.
278 Example: (bind-columns #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2 5)(3 4 6))"
279 (flet ((check-arg (x)
280 (if (not (or (sequencep x) (matrixp x)))
281 (error "bad argument type")))
282 (arg-cols (x) (if (sequencep x) 1 (num-cols x)))
283 (arg-rows (x) (if (sequencep x) (length x) (num-rows x))))
284 (dolist (x args) (check-arg x)) ;; verify data structure conformance.
285 (let ((m (arg-rows (first args)))
286 (n (arg-cols (first args))))
287 (declare (fixnum m n))
288 (dolist (x (rest args))
289 (if (/= m (arg-rows x)) (error "column lengths do not match"))
290 (incf n (arg-cols x)))
291 (do* ((result (make-array (list m n)))
292 (args args (rest args))
293 (firstcol 0)
294 (x (first args) (first args)))
295 ((null args) result)
296 (cond
297 ((sequencep x)
298 (let ((cx (make-next-element x)))
299 (dotimes (i m)
300 (setf (aref result i firstcol) (get-next-element cx i)))))
302 (let ((k (arg-cols x)))
303 (dotimes (i m)
304 (dotimes (j k)
305 (setf (aref result i (+ firstcol j)) (aref x i j)))))))
306 (incf firstcol (arg-cols x))))))
308 (defun bind-rows (&rest args)
309 "Args (&rest args)
310 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
311 along their rows.
312 Example: (bind-rows #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2)(3 4)(5 6))"
313 (flet ((check-arg (x)
314 (if (not (or (sequencep x) (matrixp x)))
315 (error "bad argument type")))
316 (arg-cols (x) (if (sequencep x) (length x) (num-cols x)))
317 (arg-rows (x) (if (sequencep x) 1 (num-rows x))))
318 (dolist (x args) (check-arg x))
319 (let ((m (arg-rows (first args)))
320 (n (arg-cols (first args))))
321 (declare (fixnum m n))
322 (dolist (x (rest args))
323 (if (/= n (arg-cols x)) (error "row lengths do not match"))
324 (incf m (arg-rows x)))
325 (do* ((result (make-array (list m n)))
326 (args args (rest args))
327 (firstrow 0)
328 (x (first args) (first args)))
329 ((null args) result)
330 (cond
331 ((sequencep x)
332 (let ((cx (make-next-element x)))
333 (dotimes (i n)
334 (setf (aref result firstrow i) (get-next-element cx i)))))
336 (let ((k (arg-rows x)))
337 (dotimes (i n)
338 (dotimes (j k)
339 (setf (aref result (+ firstrow j) i) (aref x j i)))))))
340 (incf firstrow (arg-rows x))))))
341 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
342 ;;;;
343 ;;;; Copying Functions
344 ;;;;
345 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
348 ;;; COPY-VECTOR function
351 (defun copy-vector (x)
352 "Args: (x)
353 Returns a copy of the vector X"
354 (copy-seq x))
357 ;;; COPY-ARRAY function
360 (defun copy-array (a)
361 "Args: (a)
362 Returns a copy of the array A"
363 (vector-to-array (copy-seq (array-data-vector a))
364 (array-dimensions a)))