2 ;;; Copyright (c) 2005--2007, by A.J. Rossini <blindglobe@gmail.com>
3 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
4 ;;; Since 1991, ANSI was finally finished. Edited for ANSI Common Lisp.
6 ;;; what this should do:
7 ;;; #1 - use CFFI (and possibly Verazanno) to import C/C++.
8 ;;; #2 - what to do for Fortran? Possibly: C <-> bridge, or CLapack?
9 ;;; problem: would be better to have access to Fortran. For
10 ;;; example, think of Doug Bates comment on reverse-calls (as
11 ;;; distinct from callbacks). It would be difficult if we don't
12 ;;; -- however, has anyone run Lapack or similar through F2CL?
13 ;;; Answer: yes, Matlisp does this.
15 ;;; #3 - Use a lisp-based matrix system drop-in? (matlisp, femlisp, clem, ...?)
19 ;;;; linalg -- Lisp-Stat interface to basic linear algebra routines.
21 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
22 ;;;; unrestricted use.
28 (defpackage :lisp-stat-linalg
36 :lisp-stat-compound-data
39 :lisp-stat-linalg-data
41 (:shadowing-import-from
:lisp-stat-math
42 expt
+ -
* / ** mod rem abs
1+ 1- log exp sqrt sin cos tan
43 asin acos atan sinh cosh tanh asinh acosh atanh float random
44 truncate floor ceiling round minusp zerop plusp evenp oddp
45 < <= = /= >= > complex conjugate realpart imagpart phase
46 min max logand logior logxor lognot ffloor fceiling
47 ftruncate fround signum cis
)
48 (:export chol-decomp lu-decomp lu-solve determinant inverse sv-decomp
49 qr-decomp rcondest make-rotation spline kernel-dens kernel-smooth
50 fft make-sweep-matrix sweep-operator ax
+y eigen
52 check-real
;; for optimize
55 (in-package #:lisp-stat-linalg
)
59 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
61 ;;; Lisp Interfaces to Linear Algebra Routines
63 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
66 ;;; Cholesky Decomposition
69 (cffi:defcfun
("ccl_chol_decomp_front" ccl-chol-decomp-front
)
70 :int
(x :pointer
) (y :int
) (z :pointer
))
71 (defun chol-decomp-front (x y z
)
72 (ccl-chol-decomp-front x y z
))
78 (cffi:defcfun
("ccl_lu_decomp_front" ccl-lu-decomp-front
)
79 :int
(x :pointer
) (y :int
) (z :pointer
) (u :int
) (v :pointer
))
80 (defun lu-decomp-front (x y z u v
)
81 (ccl-lu-decomp-front x y z u v
))
83 (cffi:defcfun
("ccl_lu_solve_front" ccl-lu-solve-front
)
84 :int
(x :pointer
) (y :int
) (z :pointer
) (u :pointer
) (v :int
))
85 (defun lu-solve-front (x y z u v
)
86 (ccl-lu-solve-front x y z u v
))
88 (cffi:defcfun
("ccl_lu_inverse_front" ccl-lu-inverse-front
)
89 :int
(x :pointer
) (y :int
) (z :pointer
) (u :pointer
) (v :int
) (w :pointer
))
90 (defun lu-inverse-front (x y z u v w
)
91 (ccl-lu-inverse-front x y z u v w
))
97 (cffi:defcfun
("ccl_sv_decomp_front" ccl-sv-decomp-front
)
98 :int
(x :pointer
) (y :int
) (z :int
) (u :pointer
) (v :pointer
))
99 (defun sv-decomp-front (x y z u v
)
100 (ccl-sv-decomp-front x y z u v
))
103 ;;;; QR Decomposition
106 (cffi:defcfun
("ccl_qr_decomp_front" ccl-qr-decomp-front
)
107 :int
(x :pointer
) (y :int
) (z :int
) (u :pointer
) (v :pointer
) (w :int
))
108 (defun qr-decomp-front (x y z u v w
)
109 (ccl-qr-decomp-front x y z u v w
))
112 ;;;; Estimate of Condition Number for Lower Triangular Matrix
115 (cffi:defcfun
("ccl_rcondest_front" ccl-rcondest-front
)
116 :double
(x :pointer
) (y :int
))
117 (defun rcondest-front (x y
)
118 (ccl-rcondest-front x y
))
121 ;;;; Make Rotation Matrix
124 (cffi:defcfun
("ccl_make_rotation_front" ccl-make-rotation-front
)
125 :int
(x :int
) (y :pointer
) (z :pointer
) (u :pointer
) (v :int
) (w :double
))
126 (defun make-rotation-front (x y z u v w
)
127 (ccl-make-rotation-front x y z u v
(float w
1d0
)))
130 ;;;; Eigenvalues and Eigenvectors
133 (cffi:defcfun
("ccl_eigen_front" ccl-eigen-front
)
134 :int
(x :pointer
) (y :int
) (z :pointer
) (u :pointer
) (v :pointer
))
135 (defun eigen-front (x y z u v
)
136 (ccl-eigen-front x y z u v
))
139 ;;;; Spline Interpolation
142 (cffi:defcfun
("ccl_range_to_rseq" ccl-range-to-rseq
)
143 :int
(x :int
) (y :pointer
) (z :int
) (u :pointer
))
144 (defun la-range-to-rseq (x y z u
)
145 (ccl-range-to-rseq x y z u
))
147 (cffi:defcfun
("ccl_spline_front" ccl-spline-front
)
148 :int
(x :int
) (y :pointer
) (z :pointer
) (u :int
) (v :pointer
) (w :pointer
) (a :pointer
))
149 (defun spline-front (x y z u v w a
)
150 (ccl-spline-front x y z u v w a
))
153 ;;;; Kernel Density Estimators and Smoothers
156 (cffi:defcfun
("ccl_kernel_dens_front" ccl-kernel-dens-front
)
157 :int
(x :pointer
) (y :int
) (z :double
) (u :pointer
) (v :pointer
) (w :int
) (a :int
))
158 (defun kernel-dens-front (x y z u v w a
)
159 (ccl-kernel-dens-front x y
(float z
1d0
) u v w a
))
161 (cffi:defcfun
("ccl_kernel_smooth_front" ccl-kernel-smooth-front
)
162 :int
(x :pointer
) (y :pointer
) (z :int
) (u :double
) (v :pointer
) (w :pointer
) (a :int
) (b :int
))
163 (defun kernel-smooth-front (x y z u v w a b
)
164 (ccl-kernel-smooth-front x y z
(float u
1d0
) v w a b
))
167 ;;;; Lowess Smoother Interface
170 (cffi:defcfun
("ccl_base_lowess_front" ccl-base-lowess-front
)
171 :int
(x :pointer
) (y :pointer
) (z :int
) (u :double
) (v :int
) (w :double
) (a :pointer
) (b :pointer
) (c :pointer
))
172 (defun base-lowess-front (x y z u v w a b c
)
173 (ccl-base-lowess-front x y z
(float u
1d0
) v
(float w
1d0
) a b c
))
179 (cffi:defcfun
("ccl_fft_front" ccl-fft-front
)
180 :int
(x :int
) (y :pointer
) (z :pointer
) (u :int
))
181 (defun fft-front (x y z u
)
182 (ccl-fft-front x y z u
))
187 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
189 ;;;; Lisp to C number conversion and checking
191 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
194 ;;;; Lisp to/from C sequence and matrix conversion and checking
198 "FIXME:AJR this is not used anywhere?"
201 (defun check-fixnum (a)
202 (if (/= 0 (la-data-mode a
)) (error "not an integer sequence - ~s" a
)))
204 (defun check-real (data)
205 (let ((data (compound-data-seq data
)))
208 (let ((n (length data
)))
212 (check-one-real (aref data i
)))))
213 ((consp data
) (dolist (x data
) (check-one-real x
)))
214 (t (error "bad sequence - ~s" data
)))))
216 (defun vec-assign (a i x
) (setf (aref a i
) x
))
218 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
219 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
221 ;;;; Lisp Interfaces to Linear Algebra Routines
223 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
224 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
226 ;;; FIXME: use dpbt[f2|rf], dpbstf, dpot[f2|rf]; dpptrf, zpbstf, zpbt[f2|rf]
227 ;;; remember: factorization = decomposition, depending on training.
229 (defun chol-decomp (a &optional
(maxoffl 0.0))
231 Modified Cholesky decomposition. A should be a square, symmetric matrix.
232 Computes lower triangular matrix L such that L L^T = A + D where D is a
233 diagonal matrix. If A is strictly positive definite D will be zero.
234 Otherwise, D is as small as possible to make A + D numerically strictly
235 positive definite. Returns a list (L (max D))."
236 (check-square-matrix a
)
238 (let* ((n (array-dimension a
0))
239 (result (make-array (list n n
)))
240 (dpars (list maxoffl
0.0)))
242 (let ((mat (la-data-to-matrix a
+mode-re
+))
243 (dp (la-data-to-vector dpars
+mode-re
+)))
246 (chol-decomp-front mat n dp
)
247 (la-matrix-to-data mat n n
+mode-re
+ result
)
248 (la-vector-to-data dp
2 +mode-re
+ dpars
))
249 (la-free-matrix mat n
)
250 (la-free-vector dp
)))
251 (list result
(second dpars
))))
256 ;;; i.e. result use by:
257 ;;; (setf (values (lu-out1 lu-out2 lu-out3)) (matlisp:lu my-matrix))
258 ;;; for solution, ...
260 ;;; (matlisp:gesv a b &opt ipivot)
264 A is a square matrix of numbers (real or complex). Computes the LU
265 decomposition of A and returns a list of the form (LU IV D FLAG), where
266 LU is a matrix with the L part in the lower triangle, the U part in the
267 upper triangle (the diagonal entries of L are taken to be 1), IV is a vector
268 describing the row permutation used, D is 1 if the number of permutations
269 is odd, -1 if even, and FLAG is T if A is numerically singular, NIL otherwise.
271 (check-square-matrix a
)
272 (let* ((n (array-dimension a
0))
273 (mode (max +mode-re
+ (la-data-mode a
)))
274 (result (list (make-array (list n n
)) (make-array n
) nil nil
)))
275 (let ((mat (la-data-to-matrix a mode
))
276 (iv (la-vector n
+mode-in
+))
277 (d (la-vector 1 +mode-re
+))
281 (setf singular
(lu-decomp-front mat n iv mode d
))
282 (la-matrix-to-data mat n n mode
(first result
))
283 (la-vector-to-data iv n
+mode-in
+ (second result
))
284 (setf (third result
) (la-get-double d
0))
285 (setf (fourth result
) (if (= singular
0.0) nil t
)))
286 (la-free-matrix mat n
)
291 (defun lu-solve (lu lb
)
293 LU is the result of (LU-DECOMP A) for a square matrix A, B is a sequence.
294 Returns the solution to the equation Ax = B. Signals an error if A is
296 (let ((la (first lu
))
298 (check-square-matrix la
)
299 (check-sequence lidx
)
302 (let* ((n (num-rows la
))
303 (result (make-sequence (if (consp lb
) 'list
'vector
) n
))
304 (a-mode (la-data-mode la
))
305 (b-mode (la-data-mode lb
)))
306 (if (/= n
(length lidx
)) (error "index sequence is wrong length"))
307 (if (/= n
(length lb
)) (error "right hand side is wrong length"))
308 (let* ((mode (max +mode-re
+ a-mode b-mode
))
309 (a (la-data-to-matrix la mode
))
310 (indx (la-data-to-vector lidx
+mode-in
+))
311 (b (la-data-to-vector lb mode
))
315 (setf singular
(lu-solve-front a n indx b mode
))
316 (la-vector-to-data b n mode result
))
318 (la-free-vector indx
)
320 (if (/= 0.0 singular
) (error "matrix is (numerically) singular"))
323 (defun determinant (a)
325 Returns the determinant of the square matrix M."
326 (let* ((lu (lu-decomp a
))
332 (flet ((fabs (x) (float (abs x
) 0.d0
)))
333 (dotimes (i n
(* d1
(exp d2
)))
335 (let* ((x (aref la i i
))
337 (if (= 0.0 magn
) (return 0.d0
))
338 (setf d1
(* d1
(/ x magn
)))
339 (setf d2
(+ d2
(log magn
))))))))
343 Returns the inverse of the the square matrix M; signals an error if M is ill
344 conditioned or singular"
345 (check-square-matrix a
)
346 (let ((n (num-rows a
))
347 (mode (max +mode-re
+ (la-data-mode a
))))
349 (let ((result (make-array (list n n
) :initial-element
0)))
352 (setf (aref result i i
) 1))
353 (let ((mat (la-data-to-matrix a mode
))
354 (inv (la-data-to-matrix result mode
))
355 (iv (la-vector n
+mode-in
+))
356 (v (la-vector n mode
))
360 (setf singular
(lu-inverse-front mat n iv v mode inv
))
361 (la-matrix-to-data inv n n mode result
))
362 (la-free-matrix mat n
)
363 (la-free-matrix inv n
)
366 (if (/= singular
0) (error "matrix is (numerically) singular"))
370 ;;;; SV Decomposition
375 A is a matrix of real numbers with at least as many rows as columns.
376 Computes the singular value decomposition of A and returns a list of the form
377 (U W V FLAG) where U and V are matrices whose columns are the left and right
378 singular vectors of A and W is the sequence of singular values of A. FLAG is T
379 if the algorithm converged, NIL otherwise."
381 (let* ((m (num-rows a
))
383 (mode (max +mode-re
+ (la-data-mode a
)))
384 (result (list (make-array (list m n
))
386 (make-array (list n n
))
388 (if (< m n
) (error "number of rows less than number of columns"))
389 (if (= mode
+mode-cx
+) (error "complex SVD not available yet"))
390 (let ((mat (la-data-to-matrix a mode
))
391 (w (la-vector n
+mode-re
+))
392 (v (la-matrix n n
+mode-re
+))
396 (setf converged
(sv-decomp-front mat m n w v
))
397 (la-matrix-to-data mat m n mode
(first result
))
398 (la-vector-to-data w n mode
(second result
))
399 (la-matrix-to-data v n n mode
(third result
))
400 (setf (fourth result
) (if (/= 0.0 converged
) t nil
)))
401 (la-free-matrix mat m
)
403 (la-free-matrix v n
))
408 ;;;; QR Decomposition
411 (defun qr-decomp (a &optional pivot
)
412 "Args: (a &optional pivot)
413 A is a matrix of real numbers with at least as many rows as columns. Computes
414 the QR factorization of A and returns the result in a list of the form (Q R).
415 If PIVOT is true the columns of X are first permuted to place insure the
416 absolute values of the diagonal elements of R are nonincreasing. In this case
417 the result includes a third element, a list of the indices of the columns in
418 the order in which they were used."
420 (let* ((m (num-rows a
))
422 (mode (max +mode-re
+ (la-data-mode a
)))
425 (list (make-array (list m n
))
426 (make-array (list n n
))
428 (list (make-array (list m n
)) (make-array (list n n
))))))
429 (if (< m n
) (error "number of rows less than number of columns"))
430 (if (= mode
+mode-cx
+) (error "complex QR decomposition not available yet"))
431 (let ((mat (la-data-to-matrix a mode
))
432 (v (la-matrix n n
+mode-re
+))
433 (jpvt (la-vector n
+mode-in
+)))
436 (qr-decomp-front mat m n v jpvt p
)
437 (la-matrix-to-data mat m n mode
(first result
))
438 (la-matrix-to-data v n n mode
(second result
))
439 (if pivot
(la-vector-to-data jpvt n
+mode-in
+ (third result
))))
440 (la-free-matrix mat m
)
442 (la-free-vector jpvt
))
446 ;;;; Estimate of Condition Number for Lower Triangular Matrix
451 Returns an estimate of the reciprocal of the L1 condition number of an upper
452 triangular matrix a."
453 (check-square-matrix a
)
454 (let ((mode (max +mode-re
+ (la-data-mode a
)))
456 (if (= mode
+mode-cx
+)
457 (error "complex condition estimate not available yet"))
458 (let ((mat (la-data-to-matrix a mode
))
461 (setf est
(rcondest-front mat n
))
462 (la-free-matrix mat n
))
466 ;;;; Make Rotation Matrix
469 (defun make-rotation (x y
&optional alpha
)
470 "Args: (x y &optional alpha)
471 Returns a rotation matrix for rotating from X to Y, or from X toward Y
472 by angle ALPHA, in radians. X and Y are sequences of the same length."
475 (if alpha
(check-one-real alpha
))
476 (let* ((n (length x
))
477 (mode (max +mode-re
+ (la-data-mode x
) (la-data-mode y
)))
478 (use-angle (if alpha
1 0))
479 (angle (if alpha
(float alpha
0.0) 0.0))
480 (result (make-array (list n n
))))
481 (if (/= n
(length y
)) (error "sequences not the same length"))
482 (if (= mode
+mode-cx
+) (error "complex data not supported yet"))
483 (let ((px (la-data-to-vector x
+mode-re
+))
484 (py (la-data-to-vector y
+mode-re
+))
485 (rot (la-matrix n n
+mode-re
+)))
488 (make-rotation-front n rot px py use-angle angle
)
489 (la-matrix-to-data rot n n
+mode-re
+ result
))
492 (la-free-matrix rot n
))
496 ;;;; Eigenvalues and Vectors
501 Returns list of list of eigenvalues and list of eigenvectors of square,
502 symmetric matrix A. Third element of result is NIL if algorithm converges.
503 If the algorithm does not converge, the third element is an integer I.
504 In this case the eigenvalues 0, ..., I are not reliable."
505 (check-square-matrix a
)
506 (let ((mode (max +mode-re
+ (la-data-mode a
)))
508 (if (= mode
+mode-cx
+) (error "matrix must be real and symmetric"))
509 (let ((evals (make-array n
))
510 (evecs (make-list (* n n
)))
511 (pa (la-data-to-vector (compound-data-seq a
) +mode-re
+))
512 (w (la-vector n
+mode-re
+))
513 (z (la-vector (* n n
) +mode-re
+))
514 (fv1 (la-vector n
+mode-re
+))
518 (setf ierr
(eigen-front pa n w z fv1
))
519 (la-vector-to-data w n
+mode-re
+ evals
)
520 (la-vector-to-data z
(* n n
) +mode-re
+ evecs
))
524 (la-free-vector fv1
))
525 (list (nreverse evals
)
526 (nreverse (mapcar #'(lambda (x) (coerce x
'vector
))
527 (split-list evecs n
)))
528 (if (/= 0 ierr
) (- n ierr
))))))
531 ;;;; Spline Interpolation
534 (defun make-smoother-args (x y xvals
)
540 (unless (integerp xvals
)
541 (check-sequence xvals
)
543 (let* ((n (length x
))
544 (ns (if (integerp xvals
) xvals
(length xvals
)))
545 (result (list (make-list ns
) (make-list ns
))))
546 (if (and y
(/= n
(length y
))) (error "sequences not the same length"))
547 (list x y n
(if (integerp xvals
) 0 1) ns xvals result
)))
549 (defun get-smoother-result (args) (seventh args
))
551 (defmacro with-smoother-data
((x y xvals is-reg
) &rest body
)
558 (unless (integerp ,xvals
)
559 (check-sequence ,xvals
)
561 (let* ((supplied (not (integerp ,xvals
)))
563 (ns (if supplied
(length ,xvals
) ,xvals
))
564 (result (list (make-list ns
) (make-list ns
))))
565 (if (and ,is-reg
(/= n
(length ,y
)))
566 (error "sequences not the same length"))
567 (if (and (not supplied
) (< ns
2))
568 (error "too few points for interpolation"))
569 (let* ((px (la-data-to-vector ,x
+mode-re
+))
570 (py (if ,is-reg
(la-data-to-vector ,y
+mode-re
+)))
572 (la-data-to-vector ,xvals
+mode-re
+)
573 (la-vector ns
+mode-re
+)))
574 (pys (la-vector ns
+mode-re
+)))
575 (unless supplied
(la-range-to-rseq n px ns pxs
))
578 (la-vector-to-data pxs ns
+mode-re
+ (first result
))
579 (la-vector-to-data pys ns
+mode-re
+ (second result
)))
581 (if ,is-reg
(la-free-vector py
))
583 (la-free-vector pys
))
586 (defun spline (x y
&key
(xvals 30))
587 "Args: (x y &key xvals)
588 Returns list of x and y values of natural cubic spline interpolation of (X,Y).
589 X must be strictly increasing. XVALS can be an integer, the number of equally
590 spaced points to use in the range of X, or it can be a sequence of points at
591 which to interpolate."
592 (with-smoother-data (x y xvals t
)
593 (let ((work (la-vector (* 2 n
) +mode-re
+))
596 (setf error
(spline-front n px py ns pxs pys work
))
597 (la-free-vector work
))
598 (if (/= error
0) (error "bad data for splines")))))
601 ;;;; Kernel Density Estimators and Smoothers
604 (defun kernel-type-code (type)
605 (cond ((eq type
'u
) 0)
610 (defun kernel-dens (x &key
(type 'b
) (width -
1.0) (xvals 30))
611 "Args: (x &key xvals width type)
612 Returns list of x and y values of kernel density estimate of X. XVALS can be an
613 integer, the number of equally spaced points to use in the range of X, or it
614 can be a sequence of points at which to interpolate. WIDTH specifies the
615 window width. TYPE specifies the lernel and should be one of the symbols G, T,
616 U or B for gaussian, triangular, uniform or bisquare. The default is B."
617 (check-one-real width
)
618 (with-smoother-data (x nil xvals nil
) ;; warning about deleting unreachable code is TRUE -- 2nd arg=nil!
619 (let ((code (kernel-type-code type
))
621 (setf error
(kernel-dens-front px n width pxs pys ns code
))
622 (if (/= 0 error
) (error "bad kernel density data")))))
624 (defun kernel-smooth (x y
&key
(type 'b
) (width -
1.0) (xvals 30))
625 "Args: (x y &key xvals width type)
626 Returns list of x and y values of kernel smooth of (X,Y). XVALS can be an
627 integer, the number of equally spaced points to use in the range of X, or it
628 can be a sequence of points at which to interpolate. WIDTH specifies the
629 window width. TYPE specifies the lernel and should be one of the symbols G, T,
630 U or B for Gaussian, triangular, uniform or bisquare. The default is B."
631 (check-one-real width
)
632 (with-smoother-data (x y xvals t
)
633 (let ((code (kernel-type-code type
))
635 (kernel-smooth-front px py n width pxs pys ns code
)
636 ;; if we get the Lisp version ported from C, uncomment below and
637 ;; comment above. (thanks to Carlos Ungil for the initial CFFI
639 ;;(kernel-smooth-Cport px py n width pxs pys ns code)
640 (if (/= 0 error
) (error "bad kernel density data")))))
644 (defun kernel-smooth-Cport (px py n width
;;wts wds ;; see above for mismatch?
646 "Port of kernel_smooth (Lib/kernel.c) to Lisp.
647 FIXME:kernel-smooth-Cport : This is broken.
648 Until this is fixed, we are using Luke's C code and CFFI as glue."
650 ((and (< n
2) (<= width
0)) 1.0)
651 (t (let* ((xmin (min px
))
653 (width (/ (- xmax xmin
) (+ 1.0 (log n
)))))
654 (dotimes (i (- ns
1))
658 (dotimes (j (- n
1)) )
659 ;;;possible nasty errors...
662 ((lwidth (if wds
(* width
(aref wds j
)) width
))
663 (lwt (* (kernel-Cport (aref xs i
) (aref px j
) lwidth ktype
) ;; px?
664 (if wts
(aref wts j
) 1.0))))
665 (setf wsum
(+ wsum lwt
))
666 (setf ysum
(if py
(+ ysum
(* lwt
(aref py j
)))))) ;; py? y?
676 (defun kernel-Cport (x y w ktype
)
677 "Port of kernel() (Lib/kernel.c) to Lisp.
678 x,y,w are doubles, type is an integer"
682 (cond ((eq ktype
"B")
687 (/ (/ (* 15.0 (* (- 1.0 (* 4 z z
)) ;; k/w
688 (- 1.0 (* 4 z z
)))) ;; k/w
693 (let* ((w (* w
0.25))
695 (k (/ (exp (* -
0.5 z z
))
701 (k (if (< (abs z
) 0.5)
706 (cond ((and (> z -
1.0)
717 ;;;; Lowess Smoother Interface
720 (defun |base-lowess|
(s1 s2 f nsteps delta
)
726 (check-one-fixnum nsteps
)
727 (check-one-real delta
)
728 (let* ((n (length s1
))
729 (result (make-list n
)))
730 (if (/= n
(length s2
)) (error "sequences not the same length"))
731 (let ((x (la-data-to-vector s1
+mode-re
+))
732 (y (la-data-to-vector s2
+mode-re
+))
733 (ys (la-vector n
+mode-re
+))
734 (rw (la-vector n
+mode-re
+))
735 (res (la-vector n
+mode-re
+))
739 (setf error
(base-lowess-front x y n f nsteps delta ys rw res
))
740 (la-vector-to-data ys n
+mode-re
+ result
))
745 (la-free-vector res
))
746 (if (/= error
0) (error "bad data for lowess"))
750 static LVAL add_contour_point
(i, j
, k
, l
, x
, y
, z
, v
, result
)
760 if
((z[i][j] <= v && v < z[k][l]) || (z[k][l] <= v && v < z[i][j])) {
763 p = (v - z[i][j]) / (z[k][l] - z[i][j]);
765 rplaca(pt, cvflonum((FLOTYPE) (q * x[i] + p * x[k])));
766 rplaca
(cdr(pt), cvflonum
((FLOTYPE) (q * y
[j] + p * y[l])));
767 result = cons(pt, result);
773 LVAL xssurface_contour()
775 LVAL s1, s2, mat, result;
781 s1 = xsgetsequence();
782 s2 = xsgetsequence();
784 v = makedouble(xlgetarg());
787 n = seqlen(s1); m = seqlen(s2);
788 if (n != numrows(mat) || m != numcols(mat)) xlfail("dimensions do not match");
789 if (data_mode(s1) == CX || data_mode(s2) == CX || data_mode(mat) == CX)
790 xlfail("data must be real");
792 x = (RVector) data_to_vector(s1, RE);
793 y = (RVector) data_to_vector(s2, RE);
794 z = (RMatrix) data_to_matrix(mat, RE);
798 for (i = 0; i < n - 1; i++) {
799 for (j = 0; j < m - 1; j++) {
800 result = add_contour_point(i, j, i, j+1, x, y, z, v, result);
801 result = add_contour_point(i, j+1, i+1, j+1, x, y, z, v, result);
802 result = add_contour_point(i+1, j+1, i+1, j, x, y, z, v, result);
803 result = add_contour_point(i+1, j, i, j, x, y, z, v, result);
820 ;;; ??replace with matlisp:fft and matlisp:ifft (the latter for inverse mapping)
822 (defun fft (x &optional inverse)
823 "Args: (x &optional inverse)
824 Returns unnormalized Fourier transform of X, or inverse transform if INVERSE
827 (let* ((n (length x))
828 (mode (la-data-mode x))
829 (isign (if inverse -1 1))
830 (result (if (consp x) (make-list n) (make-array n))))
831 (let ((px (la-data-to-vector x +mode-cx+))
832 (work (la-vector (+ (* 4 n) 15) +mode-re+)))
835 (fft-front n px work isign)
836 (la-vector-to-data px n +mode-cx+ result))
838 (la-free-vector work))
842 ;;; SWEEP Operator: FIXME: use matlisp
845 (defun make-sweep-front (x y w n p mode has_w x_mean result)
846 (declare (fixnum n p mode has_w))
861 (has-w (if (/= 0 has_w) t nil))
863 (declare (long-float val dxi dyi dv dw sum_w dxik dxjk dyj
864 dx_meani dx_meanj dy_mean))
865 ;; (declare-double val dxi dyi dv dw sum_w dxik dxjk dyj
866 ;; dx_meani dx_meanj dy_mean)
868 (if (> mode RE) (error "not supported for complex data yet"))
870 (setf x_data (compound-data-seq x))
871 (setf result_data (compound-data-seq result))
873 ;; find the mean of y
878 (setf dyi (makedouble (aref y i)))
880 (setf dw (makedouble (aref w i)))
882 (setf dyi (* dyi dw)))
884 (if (not has-w) (setf sum_w (float n 0.0)))
885 (if (<= sum_w 0.0) (error "non positive sum of weights"))
886 (setf dy_mean (/ val sum_w))
888 ;; find the column means
894 (setf dxi (makedouble (aref x_data (+ (* p i) j))))
896 (setf dw (makedouble (aref w i)))
897 (setf dxi (* dxi dw)))
899 (setf (aref x_mean j) (/ val sum_w)))
901 ;; put 1/sum_w in topleft, means on left, minus means on top
902 (setf (aref result_data 0) (/ 1.0 sum_w))
905 (setf dxi (makedouble (aref x_mean i)))
906 (setf (aref result_data (+ i 1)) (- dxi))
907 (setf (aref result_data (* (+ i 1) (+ p 2))) dxi))
908 (setf (aref result_data (+ p 1)) (- dy_mean))
909 (setf (aref result_data (* (+ p 1) (+ p 2))) dy_mean)
911 ;; put sums of adjusted cross products in body
919 (setf dxik (makedouble (aref x_data (+ (* p k) i))))
920 (setf dxjk (makedouble (aref x_data (+ (* p k) j))))
921 (setf dx_meani (makedouble (aref x_mean i)))
922 (setf dx_meanj (makedouble (aref x_mean j)))
923 (setf dv (* (- dxik dx_meani) (- dxjk dx_meanj)))
925 (setf dw (makedouble (aref w k)))
928 (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ j 1))) val)
929 (setf (aref result_data (+ (* (+ j 1) (+ p 2)) (+ i 1))) val))
933 (setf dxik (makedouble (aref x_data (+ (* p j) i))))
934 (setf dyj (makedouble (aref y j)))
935 (setf dx_meani (makedouble (aref x_mean i)))
936 (setf dv (* (- dxik dx_meani) (- dyj dy_mean)))
938 (setf dw (makedouble (aref w j)))
941 (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ p 1))) val)
942 (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ i 1))) val))
946 (setf dyj (makedouble (aref y j)))
947 (setf dv (* (- dyj dy_mean) (- dyj dy_mean)))
949 (setf dw (makedouble (aref w j)))
952 (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ p 1))) val)))
954 ;;; FIXME: use matlisp
955 (defun sweep-in-place-front (a rows cols mode k tol)
956 "Sweep algorithm for linear regression."
957 (declare (long-float tol))
958 (declare (fixnum rows cols mode k))
966 (declare (long-float pivot aij aik akj akk))
968 (if (> mode RE) (error "not supported for complex data yet"))
969 (if (or (< k 0) (>= k rows) (>= k cols)) (error "index out of range"))
971 (setf tol (max tol machine-epsilon))
972 (setf data (compound-data-seq a))
974 (setf pivot (makedouble (aref data (+ (* cols k) k))))
977 ((or (> pivot tol) (< pivot (- tol)))
982 (when (and (/= i k) (/= j k))
983 (setf aij (makedouble (aref data (+ (* cols i) j))))
984 (setf aik (makedouble (aref data (+ (* cols i) k))))
985 (setf akj (makedouble (aref data (+ (* cols k) j))))
986 (setf aij (- aij (/ (* aik akj) pivot)))
987 (setf (aref data (+ (* cols i) j)) aij))))
991 (setf aik (makedouble (aref data (+ (* cols i) k))))
993 (setf aik (/ aik pivot))
994 (setf (aref data (+ (* cols i) k)) aik)))
998 (setf akj (makedouble (aref data (+ (* cols k) j))))
1000 (setf akj (- (/ akj pivot)))
1001 (setf (aref data (+ (* cols k) j)) akj)))
1003 (setf akk (/ 1.0 pivot))
1004 (setf (aref data (+ (* cols k) k)) akk)
1008 ;; FIXME: use matlisp
1009 (defun make-sweep-matrix (x y &optional w)
1010 "Args: (x y &optional weights)
1011 X is matrix, Y and WEIGHTS are sequences. Returns the sweep matrix of the
1012 (weighted) regression of Y on X"
1015 (if w (check-sequence w))
1016 (let ((n (num-rows x))
1018 (if (/= n (length y)) (error "dimensions do not match"))
1019 (if (and w (/= n (length w))) (error "dimensions do not match"))
1020 (let ((mode (max (la-data-mode x)
1022 (if w (la-data-mode w) 0)))
1023 (result (make-array (list (+ p 2) (+ p 2))))
1024 (x-mean (make-array p))
1025 (y (coerce y 'vector))
1026 (w (if w (coerce w 'vector)))
1028 (make-sweep-front x y w n p mode has-w x-mean result)
1031 (defun sweep-in-place (a k tol)
1033 (check-one-fixnum k)
1034 (check-one-real tol)
1035 (let ((rows (num-rows a))
1037 (mode (la-data-mode a)))
1038 (let ((swept (sweep-in-place-front a rows cols mode k tol)))
1039 (if (/= 0 swept) t nil))))
1041 (defun sweep-operator (a columns &optional tolerances)
1042 "Args: (a indices &optional tolerances)
1043 A is a matrix, INDICES a sequence of the column indices to be swept. Returns
1044 a list of the swept result and the list of the columns actually swept. (See
1045 MULTREG documentation.) If supplied, TOLERANCES should be a list of real
1046 numbers the same length as INDICES. An index will only be swept if its pivot
1047 element is larger than the corresponding element of TOLERANCES."
1049 (check-sequence columns)
1050 (if tolerances (check-sequence tolerances))
1052 (check-fixnum columns)
1053 (if tolerances (check-real tolerances))
1055 (result (copy-array a))
1057 (columns (coerce columns 'list) (cdr columns))
1058 (tolerances (if (consp tolerances) (coerce tolerances 'list))
1059 (if (consp tolerances) (cdr tolerances))))
1060 ((null columns) (list result swept-columns))
1061 (let ((col (first columns))
1062 (tol (if (consp tolerances) (first tolerances) tol)))
1063 (if (sweep-in-place result col tol)
1064 (setf swept-columns (cons col swept-columns))))))
1073 (defun ax+y (a x y &optional lower)
1074 "Args (a x y &optional lower)
1075 Returns (+ (matmult A X) Y). If LOWER is not nil, A is taken to be lower
1077 This could probably be made more efficient."
1078 (check-square-matrix a)
1084 (let* ((n (num-rows a))
1085 (result (make-list n))
1086 (a (compound-data-seq a)))
1087 (declare (fixnum n))
1088 (if (or (/= n (length x)) (/= n (length y)))
1089 (error "dimensions do not match"))
1090 (do* ((tx (make-next-element x) (make-next-element x))
1091 (ty (make-next-element y))
1092 (tr (make-next-element result))
1094 (start 0 (+ start n))
1095 (end (if lower (+ i 1) n) (if lower (+ i 1) n)))
1097 (declare (fixnum i start end))
1098 (let ((val (get-next-element ty i)))
1100 (declare (fixnum j))
1101 (setf val (+ val (* (get-next-element tx j)
1102 (aref a (+ start j))))))
1103 (set-next-element tr i val)))))