Update local CFFI to darcs from 1.6.08
[CommonLispStat.git] / matrices.lsp
blob9a09bee8e08db360df36f0b3661ce58108037638
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 ;; matrix -- conflicts!
27 num-rows num-cols matmult identity-matrix diagonal
28 row-list column-list inner-product outer-product
29 cross-product transpose bind-columns bind-rows
30 array-data-vector vector-to-array
32 check-matrix check-square-matrix
34 copy-array copy-vector
37 (in-package :lisp-stat-matrix)
39 (deftype matrix () 'array) ;; temp fix
41 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
42 ;;;;
43 ;;;; Array to Row-Major Data Vector Conversion Functions
44 ;;;;
45 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
47 (defun array-data-vector (a)
48 "Args: (a)
49 Displaces array A to a vector"
50 (make-array (array-total-size a)
51 :displaced-to a
52 :element-type (array-element-type a)))
54 (defun vector-to-array (v dims)
55 "Args: (v dims)
56 Displaces vector V to array with dimensions DIMS"
57 (make-array dims
58 :displaced-to v
59 :element-type (array-element-type v)))
61 ;;;;
63 (defun check-matrix (a)
64 (if (not (and (typep a' array)
65 (= (array-rank a) 2)))
66 (error "not a matrix - ~s" a)
67 t))
69 (defun check-square-matrix (a)
70 (if (and (check-matrix a)
71 (/= (array-dimension a 0) (array-dimension a 1))
72 (error "matrix not square - ~s" a))
73 t))
75 (defun matrixp (x)
76 "Args: (x)
78 Returns:
79 T if X is a 2-d array (i.e. a matrix),
80 NIL otherwise."
81 (and (typep x 'array)
82 (= (array-rank x) 2)))
84 (defun num-rows (x)
85 "Args: (x)
87 Returns number of rows in X."
88 (if (matrixp x)
89 (array-dimension x 0)
90 (error "only useful for matrices.")))
92 (defun num-cols (x)
93 "Args: (x)
95 Returns number of columns in X."
96 (if (matrixp x)
97 (array-dimension x 1)
98 (error "only useful for matrices.")))
101 ;;; Look at this! Prime target for generic function dispatch!
102 (defun matmult (a b &rest args)
103 "Args: (a b &rest args)
104 Returns the matrix product of matrices a, b, etc. If a is a vector it is
105 treated as a row vector; if b is a vector it is treated as a column vector."
106 ;; fixme: why does SBCL claim this is unreachable?
107 (let ((rtype (cond ((and (typep a 'matrix)
108 (typep b 'matrix)) 'matrix)
109 ((and (typep a 'matrix)
110 (typep b 'sequence)) 'vector)
111 ((and (typep a 'sequence)
112 (typep b 'matrix)) 'vector)
113 ((and (typep a 'sequence)
114 (typep b 'sequence)) 'number)
115 ((typep a 'sequence)
116 (if (consp a) 'list 'vector))
117 ((typep b 'sequence)
118 (if (consp b) 'list 'vector)))))
120 (if (typep a 'sequence)
121 (setf a (vector-to-array (coerce a 'vector) (list 1 (length a)))))
122 (if (typep b 'sequence)
123 (setf b (vector-to-array (coerce b 'vector) (list (length b) 1))))
124 (if (not (= (array-dimension a 1) (array-dimension b 0)))
125 (error "dimensions do not match"))
126 (if args
127 (reduce #'matmult args :initial-value (matmult a b))
128 (let* ((n (array-dimension a 0))
129 (m (array-dimension b 1))
130 (p (array-dimension a 1))
131 (c (make-array (list n m)))
133 (declare (fixnum n m p))
134 (dotimes (i n)
135 (declare (fixnum i))
136 (dotimes (j m)
137 (declare (fixnum j))
138 (setq x 0)
139 (dotimes (k p)
140 (declare (fixnum k))
141 (setq x (+ x
142 (* (aref a i k)
143 (aref b k j)))))
144 (setf (aref c i j) x)))
145 (case rtype
146 (matrix c)
147 (number (aref c 0 0))
148 (t (coerce (compound-data-seq c) rtype)))))))
150 (defun identity-matrix (n)
151 "Args: (n)
152 Returns the identity matrix of rank N."
153 (let ((result (make-array (list n n) :initial-element 0)))
154 (dotimes (i n result)
155 (declare (fixnum i))
156 (setf (aref result i i) 1))))
158 ;; this thing is not very efficient at this point - too much coercing
159 (defun diagonal (x)
160 "Args: (x)
161 If X is a matrix, returns the diagonal of X. If X is a sequence, returns a
162 diagonal matrix of rank (length X) with diagonal elements eq to the elements
163 of X."
164 (cond ((typep x 'matrix)
165 (let* ((n (min (num-rows x) (num-cols x)))
166 (result (make-array n)))
167 (dotimes (i n (coerce result 'list))
168 (setf (aref result i) (aref x i i)))))
169 ((typep x 'sequence)
170 (let* ((x (coerce x 'vector))
171 (n (length x))
172 (result (make-array (list n n) :initial-element 0)))
173 (dotimes (i n result)
174 (setf (aref result i i) (aref x i)))))
175 (t (error "argument must be a matrix or a sequence"))))
177 (defun row-list (x)
178 "Args: (m)
180 Returns a list of the rows of M as vectors"
181 (check-matrix x)
182 (let ((m (num-rows x))
183 (n (num-cols x))
184 (result nil))
185 (declare (fixnum m n))
186 (flet ((get-row (k)
187 (declare (fixnum k))
188 (let ((row (make-array n)))
189 (dotimes (i n row)
190 (declare (fixnum i))
191 (setf (aref row i) (aref x k i))))))
192 (dotimes (i m result)
193 (declare (fixnum i))
194 (setf result (cons (get-row (- m i 1)) result))))))
196 (defun column-list (x)
197 "Args: (m)
199 Returns a list of the columns of M as vectors"
200 (check-matrix x)
201 (let ((m (num-rows x))
202 (n (num-cols x))
203 (result nil))
204 (declare (fixnum m n))
205 (flet ((get-col (k)
206 (declare (fixnum k))
207 (let ((col (make-array m)))
208 (dotimes (i m col)
209 (declare (fixnum i))
210 (setf (aref col i) (aref x i k))))))
211 (dotimes (i n result)
212 (declare (fixnum i))
213 (setf result (cons (get-col (- n i 1)) result))))))
215 (defun inner-product (x y)
216 "Args: (x y)
218 Returns inner product of sequences X and Y."
219 (check-sequence x)
220 (check-sequence y)
221 (let ((n (length x))
222 (cx (make-next-element x))
223 (cy (make-next-element y))
224 (result 0))
225 (declare (fixnum n))
226 (if (/= n (length y)) (error "sequence lengths do not match"))
227 (dotimes (i n result)
228 (declare (fixnum i))
229 (setf result
230 (+ result (* (get-next-element cx i)
231 (get-next-element cy i)))))))
233 (defun outer-product (x y &optional (f #'*))
234 "Args: (x y &optional (fcn #'*))
236 Returns the generalized outer product of x and y, using fcn. Tat is, the result
237 is a matrix of dimension ((length x) (length y)) and the (i j) element of the
238 result is computed as (apply fcn (aref x i) (aref y j))."
240 (let* ((x (coerce x 'vector))
241 (y (coerce y 'vector))
242 (m (length x))
243 (n (length y))
244 (a (make-array (list m n))))
245 (declare (fixnum m n))
246 (dotimes (i m a)
247 (declare (fixnum i))
248 (dotimes (j n)
249 (declare (fixnum j))
250 (setf (aref a i j) (funcall f (aref x i) (aref y j)))))))
252 (defun cross-product (x)
253 "Args: (x)
255 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
256 (inner-product X X)."
258 (check-matrix x)
259 (let* ((n (num-rows x))
260 (p (num-cols x))
261 (c (make-array (list p p))))
262 (declare (fixnum n p))
263 (dotimes (i p c)
264 (declare (fixnum i))
265 (dotimes (j (+ i 1))
266 (declare (fixnum j))
267 (let ((val 0))
268 (dotimes (k n)
269 (declare (fixnum k))
270 (incf val (* (aref x k i) (aref x k j))))
271 (setf (aref c i j) val)
272 (setf (aref c j i) val))))))
274 (defun transpose-list (x)
275 (let ((m (length (first x))))
276 (dolist (next x)
277 (if (not (consp next)) (error "not a list - ~a" x))
278 (if (/= m (length next)) (error "sublists not the same length")))
279 (do* ((cx (copy-list x))
280 (result (make-list m))
281 (next result (cdr next)))
282 ((null next) result)
283 (setf (first next) (mapcar #'first cx))
284 (do ((next cx (cdr next)))
285 ((null next))
286 (setf (first next) (rest (first next)))))))
288 (defun transpose (x)
289 "Args: (m)
290 Returns the transpose of the matrix M."
291 (cond
292 ((consp x) (transpose-list x))
294 (check-matrix x)
295 (let* ((m (num-rows x))
296 (n (num-cols x))
297 (tx (make-array (list n m))))
298 (declare (fixnum m n))
299 (dotimes (i m tx)
300 (declare (fixnum i))
301 (dotimes (j n)
302 (declare (fixnum j))
303 (setf (aref tx j i) (aref x i j))))))))
305 (defun bind-columns (&rest args)
306 "Args (&rest args)
308 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
309 along their columns. Example:
310 (bind-columns #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2 5)(3 4 6))"
312 (flet ((check-arg (x)
313 (if (not (or (typep x 'sequence) (typep x 'matrix)))
314 (error "bad argument type")))
315 (arg-cols (x) (if (typep x 'sequence) 1 (num-cols x)))
316 (arg-rows (x) (if (typep x 'sequence) (length x) (num-rows x))))
317 (dolist (x args) (check-arg x)) ;; verify data structure conformance.
318 (let ((m (arg-rows (first args)))
319 (n (arg-cols (first args))))
320 (declare (fixnum m n))
321 (dolist (x (rest args))
322 (if (/= m (arg-rows x)) (error "column lengths do not match"))
323 (incf n (arg-cols x)))
324 (do* ((result (make-array (list m n)))
325 (args args (rest args))
326 (firstcol 0)
327 (x (first args) (first args)))
328 ((null args) result)
329 (cond
330 ((typep x 'sequence)
331 (let ((cx (make-next-element x)))
332 (dotimes (i m)
333 (setf (aref result i firstcol) (get-next-element cx i)))))
335 (let ((k (arg-cols x)))
336 (dotimes (i m)
337 (dotimes (j k)
338 (setf (aref result i (+ firstcol j)) (aref x i j)))))))
339 (incf firstcol (arg-cols x))))))
341 (defun bind-rows (&rest args)
342 "Args (&rest args)
344 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
345 along their rows. Example:
346 (bind-rows #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2)(3 4)(5 6))"
348 (flet ((check-arg (x)
349 (if (not (or (typep x 'sequence)
350 (typep x 'matrix)))
351 (error "bad argument type")))
352 (arg-cols (x) (if (typep x 'sequence) (length x) (num-cols x)))
353 (arg-rows (x) (if (typep x 'sequence) 1 (num-rows x))))
354 (dolist (x args) (check-arg x))
355 (let ((m (arg-rows (first args)))
356 (n (arg-cols (first args))))
357 (declare (fixnum m n))
358 (dolist (x (rest args))
359 (if (/= n (arg-cols x)) (error "row lengths do not match"))
360 (incf m (arg-rows x)))
361 (do* ((result (make-array (list m n)))
362 (args args (rest args))
363 (firstrow 0)
364 (x (first args) (first args)))
365 ((null args) result)
366 (cond
367 ((typep x 'sequence)
368 (let ((cx (make-next-element x)))
369 (dotimes (i n)
370 (setf (aref result firstrow i) (get-next-element cx i)))))
372 (let ((k (arg-rows x)))
373 (dotimes (i n)
374 (dotimes (j k)
375 (setf (aref result (+ firstrow j) i) (aref x j i)))))))
376 (incf firstrow (arg-rows x))))))
380 ;;; Copying Functions
383 (defun copy-vector (x)
384 "Args: (x)
386 Returns a copy of the vector X"
387 (copy-seq x))
389 (defun copy-array (a)
390 "Args: (a)
392 Returns a copy of the array A"
393 (vector-to-array (copy-seq (array-data-vector a))
394 (array-dimensions a)))
398 (defgeneric copy (x)
399 (:documentation "methods for copying linaar algebra forms."))
401 (defmethod copy ((x vector)))
403 (defmethod copy ((x matrix)))