start putting together unittesting framework for linear regression models.
[CommonLispStat.git] / external / cl-grnm / la.lisp
blob63e91a66153cacd953377ad6fc21ac60d88a7b63
1 (in-package :neldermead)
3 ;; Simple matrix-vector facility. Optimized for staying out of the way.
4 (defmacro with-matrix-dimensions (varforms &body body)
5 (let ((form `(progn ,@body)))
6 (dolist (vform varforms)
7 (setf form
8 `(destructuring-bind ,(car vform)
9 (array-dimensions ,(cadr vform))
10 (declare (fixnum ,@(car vform))
11 (ignorable ,@(car vform)))
12 ,form)))
13 form))
15 (defun ip (v w &optional (start 0))
16 (let ((l (length v))
17 (ax 0.0d0))
18 (dotimes (i (- l start))
19 (incf ax (* (aref v (+ i start))
20 (aref w (+ i start)))))
21 ax))
23 (defun norm (v &optional (start 0))
24 (sqrt (ip v v start)))
26 (defun make-matrix (m n)
27 (make-array (list m n)
28 :element-type 'double-float
29 :initial-element 0.0d0))
31 (defun make-vector (n &key (initial-element 0.0d0))
32 (make-array n
33 :element-type 'double-float
34 :initial-element initial-element))
36 (defun v*c (vector constant)
37 (let* ((n (length vector))
38 (res (make-vector n)))
39 (dotimes (i n)
40 (setf (aref res i)
41 (* constant (aref vector i))))
42 res))
44 (defun v+w*c (v w c)
45 (assert (= (length v) (length w)))
46 (let* ((n (length v))
47 (res (make-vector n)))
48 (dotimes (i n)
49 (setf (aref res i)
50 (+ (aref v i)
51 (* c (aref w i)))))
52 res))
54 (defun qr-factorization (mat &key (with-q t))
55 (declare (type (simple-array double-float (* *)) mat)
56 (optimize (speed 3) (safety 0)))
57 (with-matrix-dimensions (((m n) mat))
58 (let ((q (when with-q
59 (let ((pq (make-matrix m n)))
60 (declare (type (simple-array double-float (* *)) pq))
61 (dotimes (i (min m n))
62 (setf (aref pq i i) 1.0d0))
63 pq))))
65 (loop for i from 0 below (- m 1) do
66 (loop for j from (+ i 1) below m do
67 (let ((a (aref mat i i))
68 (b (aref mat j i)))
69 (when (/= b 0.0d0)
70 (let* ((r (sqrt (+ (* a a) (* b b))))
71 (c (/ a r))
72 (s (/ b r)))
74 (loop for k from i below m do
75 (let ((olda (aref mat i k))
76 (oldb (aref mat j k)))
78 (setf (aref mat i k)
79 (+ (* c olda) (* s oldb))
80 (aref mat j k)
81 (+ (* (- s) olda) (* c oldb)))))
83 (when with-q
84 (let ((q q))
85 (declare (type (simple-array double-float (* *)) q))
86 (loop for k from 0 below m do
87 (let ((olda (aref q k i))
88 (oldb (aref q k j)))
90 (setf (aref q k i)
91 (+ (* c olda) (* s oldb))
92 (aref q k j)
93 (+ (* (- s) olda) (* c oldb))))))))))))
95 (values mat q))))