merged from ansiClib
[CommonLispStat.git] / matrices.lsp
blob41f981d85c8867ea065bbd85d837ba1157a3ace6
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.
13 ;;(provide "matrices")
15 ;;;;
16 ;;;; Package Setup
17 ;;;;
19 ;;(in-package :lisp-stat-basics)
21 (defpackage :lisp-stat-matrix
22 (:use :common-lisp
23 :lisp-stat-sequence)
24 (:export
26 ;;(export '(
27 matrixp num-rows num-cols matmult identity-matrix diagonal
28 row-list column-list inner-product outer-product cross-product
29 transpose bind-columns bind-rows))
31 (in-package :lisp-stat-matrix)
33 (deftype matrix () 'array) ;; temp fix
35 (defun check-matrix (a)
36 (if (not (and (arrayp a) (= (array-rank a) 2)))
37 (error "not a matrix - ~s" a)
38 t))
40 (defun check-square-matrix (a)
41 (if (and (check-matrix a)
42 (/= (array-dimension a 0) (array-dimension a 1))
43 (error "matrix not square - ~s" a))
44 t))
46 (defun matrixp (x)
47 "Args: (x)
48 Returns T if X is a matrix, NIL otherwise."
49 (and (arrayp x) (= (array-rank x) 2)))
51 (defun num-rows (x)
52 "Args: (x)
53 Returns number of rows in X."
54 (array-dimension x 0))
56 (defun num-cols (x)
57 "Args: (x)
58 Returns number of columns in X."
59 (array-dimension x 1))
61 (defun matmult (a b &rest args)
62 "Args: (a b &rest args)
63 Returns the matrix product of matrices a, b, etc. If a is a vector it is
64 treated as a row vector; if b is a vector it is treated as a column vector."
65 (let ((rtype (cond ((and (matrixp a) (matrixp b)) 'matrix)
66 ((and (sequencep a) (sequencep b)) 'number)
67 ((sequencep a) (if (consp a) 'list 'vector))
68 ((sequencep b) (if (consp b) 'list 'vector)))))
70 (if (sequencep a)
71 (setf a (vector-to-array (coerce a 'vector) (list 1 (length a)))))
72 (if (sequencep b)
73 (setf b (vector-to-array (coerce b 'vector) (list (length b) 1))))
74 (if (not (= (array-dimension a 1) (array-dimension b 0)))
75 (error "dimensions do not match"))
76 (if args
77 (reduce #'matmult args :initial-value (matmult a b))
78 (let* ((n (array-dimension a 0))
79 (m (array-dimension b 1))
80 (p (array-dimension a 1))
81 (c (make-array (list n m)))
83 (declare (fixnum n m p))
84 (dotimes (i n)
85 (declare (fixnum i))
86 (dotimes (j m)
87 (declare (fixnum j))
88 (setq x 0)
89 (dotimes (k p)
90 (declare (fixnum k))
91 (setq x (+ x
92 (* (aref a i k) (aref b k j)))))
93 (setf (aref c i j) x)))
94 (case rtype
95 (matrix c)
96 (number (aref c 0 0))
97 (t (coerce (compound-data-seq c) rtype)))))))
99 (defun identity-matrix (n)
100 "Args: (n)
101 Returns the identity matrix of rank N."
102 (let ((result (make-array (list n n) :initial-element 0)))
103 (dotimes (i n result)
104 (declare (fixnum i))
105 (setf (aref result i i) 1))))
107 ;; this thing is not very efficient at this point - too much coercing
108 (defun diagonal (x)
109 "Args: (x)
110 If X is a matrix, returns the diagonal of X. If X is a sequence, returns a
111 diagonal matrix of rank (length X) with diagonal elements eq to the elements
112 of X."
113 (cond ((matrixp x)
114 (let* ((n (min (num-rows x) (num-cols x)))
115 (result (make-array n)))
116 (dotimes (i n (coerce result 'list))
117 (setf (aref result i) (aref x i i)))))
118 ((sequencep x)
119 (let* ((x (coerce x 'vector))
120 (n (length x))
121 (result (make-array (list n n) :initial-element 0)))
122 (dotimes (i n result)
123 (setf (aref result i i) (aref x i)))))
124 (t (error "argument must be a matrix or a sequence"))))
126 (defun row-list (x)
127 "Args: (m)
128 Returns a list of the rows of M as vectors"
129 (check-matrix x)
130 (let ((m (num-rows x))
131 (n (num-cols x))
132 (result nil))
133 (declare (fixnum m n))
134 (flet ((get-row (k)
135 (declare (fixnum k))
136 (let ((row (make-array n)))
137 (dotimes (i n row)
138 (declare (fixnum i))
139 (setf (aref row i) (aref x k i))))))
140 (dotimes (i m result)
141 (declare (fixnum i))
142 (setf result (cons (get-row (- m i 1)) result))))))
144 (defun column-list (x)
145 "Args: (m)
146 Returns a list of the columns of M as vectors"
147 (check-matrix x)
148 (let ((m (num-rows x))
149 (n (num-cols x))
150 (result nil))
151 (declare (fixnum m n))
152 (flet ((get-col (k)
153 (declare (fixnum k))
154 (let ((col (make-array m)))
155 (dotimes (i m col)
156 (declare (fixnum i))
157 (setf (aref col i) (aref x i k))))))
158 (dotimes (i n result)
159 (declare (fixnum i))
160 (setf result (cons (get-col (- n i 1)) result))))))
162 (defun inner-product (x y)
163 "Args: (x y)
164 Returns inner product of sequences X and Y."
165 (check-sequence x)
166 (check-sequence y)
167 (let ((n (length x))
168 (cx (make-next-element x))
169 (cy (make-next-element y))
170 (result 0))
171 (declare (fixnum n))
172 (if (/= n (length y)) (error "sequence lengths do not match"))
173 (dotimes (i n result)
174 (declare (fixnum i))
175 (setf result
176 (+ result (* (get-next-element cx i) (get-next-element cy i)))))))
178 (defun outer-product (x y &optional (f #'*))
179 "Args: (x y &optional (fcn #'*))
180 Returns the generalized outer product of x and y, using fcn. Tat is, the result
181 is a matrix of dimension ((length x) (length y)) and the (i j) element of the
182 result is computed as (apply fcn (aref x i) (aref y j))."
183 (let* ((x (coerce x 'vector))
184 (y (coerce y 'vector))
185 (m (length x))
186 (n (length y))
187 (a (make-array (list m n))))
188 (declare (fixnum m n))
189 (dotimes (i m a)
190 (declare (fixnum i))
191 (dotimes (j n)
192 (declare (fixnum j))
193 (setf (aref a i j) (funcall f (aref x i) (aref y j)))))))
195 (defun cross-product (x)
196 "Args: (x)
197 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
198 (inner-product X X)."
199 (check-matrix x)
200 (let* ((n (num-rows x))
201 (p (num-cols x))
202 (c (make-array (list p p))))
203 (declare (fixnum n p))
204 (dotimes (i p c)
205 (declare (fixnum i))
206 (dotimes (j (+ i 1))
207 (declare (fixnum j))
208 (let ((val 0))
209 (dotimes (k n)
210 (declare (fixnum k))
211 (incf val (* (aref x k i) (aref x k j))))
212 (setf (aref c i j) val)
213 (setf (aref c j i) val))))))
215 (defun transpose-list (x)
216 (let ((m (length (first x))))
217 (dolist (next x)
218 (if (not (consp next)) (error "not a list - ~a" x))
219 (if (/= m (length next)) (error "sublists not the same length")))
220 (do* ((cx (copy-list x))
221 (result (make-list m))
222 (next result (cdr next)))
223 ((null next) result)
224 (setf (first next) (mapcar #'first cx))
225 (do ((next cx (cdr next)))
226 ((null next))
227 (setf (first next) (rest (first next)))))))
229 (defun transpose (x)
230 "Args: (m)
231 Returns the transpose of the matrix M."
232 (cond
233 ((consp x) (transpose-list x))
235 (check-matrix x)
236 (let* ((m (num-rows x))
237 (n (num-cols x))
238 (tx (make-array (list n m))))
239 (declare (fixnum m n))
240 (dotimes (i m tx)
241 (declare (fixnum i))
242 (dotimes (j n)
243 (declare (fixnum j))
244 (setf (aref tx j i) (aref x i j))))))))
246 (defun bind-columns (&rest args)
247 "Args (&rest args)
248 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
249 along their columns.
250 Example: (bind-columns #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2 5)(3 4 6))"
251 (flet ((check-arg (x)
252 (if (not (or (sequencep x) (matrixp x)))
253 (error "bad argument type")))
254 (arg-cols (x) (if (sequencep x) 1 (num-cols x)))
255 (arg-rows (x) (if (sequencep x) (length x) (num-rows x))))
256 (dolist (x args) (check-arg x))
257 (let ((m (arg-rows (first args)))
258 (n (arg-cols (first args))))
259 (declare (fixnum m n))
260 (dolist (x (rest args))
261 (if (/= m (arg-rows x)) (error "column lengths do not match"))
262 (incf n (arg-cols x)))
263 (do* ((result (make-array (list m n)))
264 (args args (rest args))
265 (firstcol 0)
266 (x (first args) (first args)))
267 ((null args) result)
268 (cond
269 ((sequencep x)
270 (let ((cx (make-next-element x)))
271 (dotimes (i m)
272 (setf (aref result i firstcol) (get-next-element cx i)))))
274 (let ((k (arg-cols x)))
275 (dotimes (i m)
276 (dotimes (j k)
277 (setf (aref result i (+ firstcol j)) (aref x i j)))))))
278 (incf firstcol (arg-cols x))))))
280 (defun bind-rows (&rest args)
281 "Args (&rest args)
282 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
283 along their rows.
284 Example: (bind-rows #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2)(3 4)(5 6))"
285 (flet ((check-arg (x)
286 (if (not (or (sequencep x) (matrixp x)))
287 (error "bad argument type")))
288 (arg-cols (x) (if (sequencep x) (length x) (num-cols x)))
289 (arg-rows (x) (if (sequencep x) 1 (num-rows x))))
290 (dolist (x args) (check-arg x))
291 (let ((m (arg-rows (first args)))
292 (n (arg-cols (first args))))
293 (declare (fixnum m n))
294 (dolist (x (rest args))
295 (if (/= n (arg-cols x)) (error "row lengths do not match"))
296 (incf m (arg-rows x)))
297 (do* ((result (make-array (list m n)))
298 (args args (rest args))
299 (firstrow 0)
300 (x (first args) (first args)))
301 ((null args) result)
302 (cond
303 ((sequencep x)
304 (let ((cx (make-next-element x)))
305 (dotimes (i n)
306 (setf (aref result firstrow i) (get-next-element cx i)))))
308 (let ((k (arg-rows x)))
309 (dotimes (i n)
310 (dotimes (j k)
311 (setf (aref result (+ firstrow j) i) (aref x j i)))))))
312 (incf firstrow (arg-rows x))))))