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
)
20 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
22 ;;;; Array to Row-Major Data Vector Conversion Functions
24 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
26 (defun array-data-vector (a)
28 Displaces array A to a vector"
29 (make-array (array-total-size a
)
31 :element-type
(array-element-type a
)))
33 (defun vector-to-array (v dims
)
35 Displaces vector V to array with dimensions DIMS"
38 :element-type
(array-element-type v
)))
41 (defun check-square-matrix (a)
42 (if (and (check-matrix a
)
43 (/= (array-dimension a
0) (array-dimension a
1))
44 (error "matrix not square - ~s" a
))
47 (defun inner-product (x y
)
50 Returns inner product of sequences X and Y."
54 (cx (make-next-element x
))
55 (cy (make-next-element y
))
58 (if (/= n
(length y
)) (error "sequence lengths do not match"))
62 (+ result
(* (get-next-element cx i
)
63 (get-next-element cy i
)))))))
65 (defun outer-product (x y
&optional
(f #'*))
66 "Args: (x y &optional (fcn #'*))
68 Returns the generalized outer product of x and y, using fcn. Tat is, the result
69 is a matrix of dimension ((length x) (length y)) and the (i j) element of the
70 result is computed as (apply fcn (aref x i) (aref y j))."
72 (let* ((x (coerce x
'vector
))
73 (y (coerce y
'vector
))
76 (a (make-array (list m n
))))
77 (declare (fixnum m n
))
82 (setf (aref a i j
) (funcall f
(aref x i
) (aref y j
)))))))
84 (defun cross-product (x)
87 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
91 (let* ((n (num-rows x
))
93 (c (make-array (list p p
))))
94 (declare (fixnum n p
))
102 (incf val
(* (aref x k i
) (aref x k j
))))
103 (setf (aref c i j
) val
)
104 (setf (aref c j i
) val
))))))
106 (defun transpose-list (x)
107 (let ((m (length (first x
))))
109 (if (not (consp next
)) (error "not a list - ~a" x
))
110 (if (/= m
(length next
)) (error "sublists not the same length")))
111 (do* ((cx (copy-list x
))
112 (result (make-list m
))
113 (next result
(cdr next
)))
115 (setf (first next
) (mapcar #'first cx
))
116 (do ((next cx
(cdr next
)))
118 (setf (first next
) (rest (first next
)))))))
122 Returns the transpose of the matrix M."
124 ((consp x
) (transpose-list x
))
127 (let* ((m (num-rows x
))
129 (tx (make-array (list n m
))))
130 (declare (fixnum m n
))
135 (setf (aref tx j i
) (aref x i j
))))))))