3 ;;; Copyright (c) 2007--, by A.J. Rossini <blindglobe@gmail.com>
4 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
6 ;;; What is this talk of 'release'? Klingons do not make software
7 ;;; 'releases'. Our software 'escapes', leaving a bloody trail of
8 ;;; designers and quality assurance people in its wake.
10 (in-package :cls-matrix
)
12 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
14 ;;;; Array to Row-Major Data Vector Conversion Functions
16 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
18 (defun array-data-vector (a)
20 Displaces array A to a vector"
21 (make-array (array-total-size a
)
23 :element-type
(array-element-type a
)))
25 (defun vector-to-array (v dims
)
27 Displaces vector V to array with dimensions DIMS"
28 ;;; Yes, but using row or column first approach?
31 :element-type
(array-element-type v
)))
33 (defun lisp-array-matrix-p (a)
34 "Verifies that a lisp array could be considered a matrix. FIXME: Do
35 we care whether the entries are numeric or even commonly-numeric?"
37 (assert (= (array-rank a
) 2) (a) "lisp array is not 2-dim: ~A" a
)
40 (defun check-square-matrix (a)
41 (if (and (lisp-array-matrix-p a
)
42 (/= (array-dimension a
0) (array-dimension a
1))
43 (error "matrix not square - ~s" a
))
48 Returns T if X is a matrix, NIL otherwise."
49 (and (arrayp x
) (= (array-rank x
) 2)))
53 Returns number of rows in X."
54 (array-dimension x
0))
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 (typep a
'sequence
) (typep b
'sequence
)) 'number
)
67 ((typep a
'sequence
) (if (consp a
) 'list
'vector
))
68 ((typep b
'sequence
) (if (consp b
) 'list
'vector
)))))
70 (if (typep a
'sequence
)
71 (setf a
(vector-to-array (coerce a
'vector
) (list 1 (length a
)))))
72 (if (typep b
'sequence
)
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"))
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
))
92 (* (aref a i k
) (aref b k j
)))))
93 (setf (aref c i j
) x
)))
97 (t (coerce (compound-data-seq c
) rtype
)))))))
99 (defun identity-matrix (n)
100 "Return the identity matrix of rank N as a fixnum-valued array."
101 (let ((result (make-array (list n n
) :initial-element
0)))
102 (dotimes (i n result
)
104 (setf (aref result i i
) 1))))
107 "If X is a square matrix, returns the diagonal of X. If MxN,
108 M>N, return the diagonal of the NxN upper-left sub-matrix. If X is a
109 sequence, returns a diagonal matrix of rank (length X) with diagonal
110 elements eq to the elements of X. Return type is a vector. FIXME:
111 could be made more efficient?"
113 (let* ((n (min (num-rows x
) (num-cols x
)))
114 (result (make-array n
)))
115 (dotimes (i n
(coerce result
'list
))
116 (setf (aref result i
) (aref x i i
)))))
118 (let* ((x (coerce x
'vector
))
120 (result (make-array (list n n
) :initial-element
0)))
121 (dotimes (i n result
)
122 (setf (aref result i i
) (aref x i
)))))
123 (t (error "argument must be a matrix or a sequence"))))
126 "Returns a list of the rows of X as vectors"
127 (lisp-array-matrix-p x
)
128 (let ((m (num-rows x
))
131 (declare (fixnum m n
))
134 (let ((row (make-array n
)))
137 (setf (aref row i
) (aref x k i
))))))
138 (dotimes (i m result
)
140 (setf result
(cons (get-row (- m i
1)) result
))))))
142 (defun column-list (x)
144 Returns a list of vectors, each vector represents a column of M."
145 (lisp-array-matrix-p x
)
146 (let ((m (num-rows x
))
149 (declare (fixnum m n
))
152 (let ((col (make-array m
)))
155 (setf (aref col i
) (aref x i k
))))))
156 (dotimes (i n result
)
158 (setf result
(cons (get-col (- n i
1)) result
))))))
160 (defun inner-product (x y
)
162 Returns inner product of sequences X and Y."
163 (check-type x sequence
)
164 (check-type y sequence
)
165 (assert (= (length x
) (length y
))
167 "INNER-PRODUCT: sequences of different length")
169 (cx (make-next-element x
))
170 (cy (make-next-element y
))
173 (dotimes (i n 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
))
187 (a (make-array (list m n
))))
188 (declare (fixnum m n
))
193 (setf (aref a i j
) (funcall f
(aref x i
) (aref y j
)))))))
195 (defun cross-product (x)
197 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
198 (inner-product X X). This function works with raw LISP ARRAYS."
199 (lisp-array-matrix-p x
)
200 (let* ((n (num-rows x
))
202 (c (make-array (list p p
))))
203 (declare (fixnum n p
))
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
))))
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
)))
224 (setf (first next
) (mapcar #'first cx
))
225 (do ((next cx
(cdr next
)))
227 (setf (first next
) (rest (first next
)))))))
231 Returns the transpose of the matrix M."
233 ((consp x
) (transpose-list x
))
235 (lisp-array-matrix-p x
)
236 (let* ((m (num-rows x
))
238 (tx (make-array (list n m
))))
239 (declare (fixnum m n
))
244 (setf (aref tx j i
) (aref x i j
))))))))
246 (defun bind-columns (&rest args
)
248 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
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 (typep x
'sequence
) (matrixp x
)))
253 (error "bad argument type")))
254 (arg-cols (x) (if (typep x
'sequence
) 1 (num-cols x
)))
255 (arg-rows (x) (if (typep x
'sequence
) (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
))
266 (x (first args
) (first args
)))
270 (let ((cx (make-next-element x
)))
272 (setf (aref result i firstcol
) (get-next-element cx i
)))))
274 (let ((k (arg-cols x
)))
277 (setf (aref result i
(+ firstcol j
)) (aref x i j
)))))))
278 (incf firstcol
(arg-cols x
))))))
280 (defun bind-rows (&rest args
)
282 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
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 (typep x
'sequence
) (matrixp x
)))
287 (error "bad argument type")))
288 (arg-cols (x) (if (typep x
'sequence
) (length x
) (num-cols x
)))
289 (arg-rows (x) (if (typep x
'sequence
) 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
))
300 (x (first args
) (first args
)))
304 (let ((cx (make-next-element x
)))
306 (setf (aref result firstrow i
) (get-next-element cx i
)))))
308 (let ((k (arg-rows x
)))
311 (setf (aref result
(+ firstrow j
) i
) (aref x j i
)))))))
312 (incf firstrow
(arg-rows x
))))))