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
8 ;;;; matrices -- Basic matrix operations
10 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
11 ;;;; unrestricted use.
14 ;;; #1 - Need to extend to use lisp-matrix?
15 ;;; #2 - do as a second alternative?
17 (in-package :lisp-stat-matrix
)
19 (deftype matrix
() 'array
) ;; temp fix
21 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
23 ;;;; Array to Row-Major Data Vector Conversion Functions
25 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
27 (defun array-data-vector (a)
29 Displaces array A to a vector"
30 (make-array (array-total-size a
)
32 :element-type
(array-element-type a
)))
34 (defun vector-to-array (v dims
)
36 Displaces vector V to array with dimensions DIMS"
39 :element-type
(array-element-type v
)))
43 (defun check-matrix (a)
44 (if (not (and (typep a
' array
)
45 (= (array-rank a
) 2)))
46 (error "not a matrix - ~s" a
)
49 (defun check-square-matrix (a)
50 (if (and (check-matrix a
)
51 (/= (array-dimension a
0) (array-dimension a
1))
52 (error "matrix not square - ~s" a
))
59 T if X is a 2-d array (i.e. a matrix),
62 (= (array-rank x
) 2)))
67 Returns number of rows in X."
70 (error "only useful for matrices.")))
75 Returns number of columns in X."
78 (error "only useful for matrices.")))
81 ;;; Look at this! Prime target for generic function dispatch!
82 (defun matmult (a b
&rest args
)
83 "Args: (a b &rest args)
84 Returns the matrix product of matrices a, b, etc. If a is a vector it is
85 treated as a row vector; if b is a vector it is treated as a column vector."
86 ;; fixme: why does SBCL claim this is unreachable?
87 (let ((rtype (cond ((and (typep a
'matrix
)
88 (typep b
'matrix
)) 'matrix
)
89 ((and (typep a
'matrix
)
90 (typep b
'sequence
)) 'vector
)
91 ((and (typep a
'sequence
)
92 (typep b
'matrix
)) 'vector
)
93 ((and (typep a
'sequence
)
94 (typep b
'sequence
)) 'number
)
96 (if (consp a
) 'list
'vector
))
98 (if (consp b
) 'list
'vector
)))))
100 (if (typep a
'sequence
)
101 (setf a
(vector-to-array (coerce a
'vector
) (list 1 (length a
)))))
102 (if (typep b
'sequence
)
103 (setf b
(vector-to-array (coerce b
'vector
) (list (length b
) 1))))
104 (if (not (= (array-dimension a
1) (array-dimension b
0)))
105 (error "dimensions do not match"))
107 (reduce #'matmult args
:initial-value
(matmult a b
))
108 (let* ((n (array-dimension a
0))
109 (m (array-dimension b
1))
110 (p (array-dimension a
1))
111 (c (make-array (list n m
)))
113 (declare (fixnum n m p
))
124 (setf (aref c i j
) x
)))
127 (number (aref c
0 0))
128 (t (coerce (compound-data-seq c
) rtype
)))))))
130 (defun identity-matrix (n)
132 Returns the identity matrix of rank N."
133 (let ((result (make-array (list n n
) :initial-element
0)))
134 (dotimes (i n result
)
136 (setf (aref result i i
) 1))))
138 ;; this thing is not very efficient at this point - too much coercing
141 If X is a matrix, returns the diagonal of X. If X is a sequence, returns a
142 diagonal matrix of rank (length X) with diagonal elements eq to the elements
144 (cond ((typep x
'matrix
)
145 (let* ((n (min (num-rows x
) (num-cols x
)))
146 (result (make-array n
)))
147 (dotimes (i n
(coerce result
'list
))
148 (setf (aref result i
) (aref x i i
)))))
150 (let* ((x (coerce x
'vector
))
152 (result (make-array (list n n
) :initial-element
0)))
153 (dotimes (i n result
)
154 (setf (aref result i i
) (aref x i
)))))
155 (t (error "argument must be a matrix or a sequence"))))
160 Returns a list of the rows of M as vectors"
162 (let ((m (num-rows x
))
165 (declare (fixnum m n
))
168 (let ((row (make-array n
)))
171 (setf (aref row i
) (aref x k i
))))))
172 (dotimes (i m result
)
174 (setf result
(cons (get-row (- m i
1)) result
))))))
176 (defun column-list (x)
179 Returns a list of the columns of M as vectors"
181 (let ((m (num-rows x
))
184 (declare (fixnum m n
))
187 (let ((col (make-array m
)))
190 (setf (aref col i
) (aref x i k
))))))
191 (dotimes (i n result
)
193 (setf result
(cons (get-col (- n i
1)) result
))))))
195 (defun inner-product (x y
)
198 Returns inner product of sequences X and Y."
202 (cx (make-next-element x
))
203 (cy (make-next-element y
))
206 (if (/= n
(length y
)) (error "sequence lengths do not match"))
207 (dotimes (i n result
)
210 (+ result
(* (get-next-element cx i
)
211 (get-next-element cy i
)))))))
213 (defun outer-product (x y
&optional
(f #'*))
214 "Args: (x y &optional (fcn #'*))
216 Returns the generalized outer product of x and y, using fcn. Tat is, the result
217 is a matrix of dimension ((length x) (length y)) and the (i j) element of the
218 result is computed as (apply fcn (aref x i) (aref y j))."
220 (let* ((x (coerce x
'vector
))
221 (y (coerce y
'vector
))
224 (a (make-array (list m n
))))
225 (declare (fixnum m n
))
230 (setf (aref a i j
) (funcall f
(aref x i
) (aref y j
)))))))
232 (defun cross-product (x)
235 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
236 (inner-product X X)."
239 (let* ((n (num-rows x
))
241 (c (make-array (list p p
))))
242 (declare (fixnum n p
))
250 (incf val
(* (aref x k i
) (aref x k j
))))
251 (setf (aref c i j
) val
)
252 (setf (aref c j i
) val
))))))
254 (defun transpose-list (x)
255 (let ((m (length (first x
))))
257 (if (not (consp next
)) (error "not a list - ~a" x
))
258 (if (/= m
(length next
)) (error "sublists not the same length")))
259 (do* ((cx (copy-list x
))
260 (result (make-list m
))
261 (next result
(cdr next
)))
263 (setf (first next
) (mapcar #'first cx
))
264 (do ((next cx
(cdr next
)))
266 (setf (first next
) (rest (first next
)))))))
270 Returns the transpose of the matrix M."
272 ((consp x
) (transpose-list x
))
275 (let* ((m (num-rows x
))
277 (tx (make-array (list n m
))))
278 (declare (fixnum m n
))
283 (setf (aref tx j i
) (aref x i j
))))))))
285 (defun bind-columns (&rest args
)
288 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
289 along their columns. Example:
290 (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
))
307 (x (first args
) (first args
)))
311 (let ((cx (make-next-element x
)))
313 (setf (aref result i firstcol
) (get-next-element cx i
)))))
315 (let ((k (arg-cols x
)))
318 (setf (aref result i
(+ firstcol j
)) (aref x i j
)))))))
319 (incf firstcol
(arg-cols x
))))))
321 (defun bind-rows (&rest args
)
324 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
325 along their rows. Example:
326 (bind-rows #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2)(3 4)(5 6))"
328 (flet ((check-arg (x)
329 (if (not (or (typep x
'sequence
)
331 (error "bad argument type")))
332 (arg-cols (x) (if (typep x
'sequence
) (length x
) (num-cols x
)))
333 (arg-rows (x) (if (typep x
'sequence
) 1 (num-rows x
))))
334 (dolist (x args
) (check-arg x
))
335 (let ((m (arg-rows (first args
)))
336 (n (arg-cols (first args
))))
337 (declare (fixnum m n
))
338 (dolist (x (rest args
))
339 (if (/= n
(arg-cols x
)) (error "row lengths do not match"))
340 (incf m
(arg-rows x
)))
341 (do* ((result (make-array (list m n
)))
342 (args args
(rest args
))
344 (x (first args
) (first args
)))
348 (let ((cx (make-next-element x
)))
350 (setf (aref result firstrow i
) (get-next-element cx i
)))))
352 (let ((k (arg-rows x
)))
355 (setf (aref result
(+ firstrow j
) i
) (aref x j i
)))))))
356 (incf firstrow
(arg-rows x
))))))
360 ;;; Copying Functions
363 (defun copy-vector (x)
366 Returns a copy of the vector X"
369 (defun copy-array (a)
372 Returns a copy of the array A"
373 (vector-to-array (copy-seq (array-data-vector a
))
374 (array-dimensions a
)))
379 (:documentation
"methods for copying linaar algebra forms."))
381 (defmethod copy ((x vector
)))
383 (defmethod copy ((x matrix
)))