crude hack to allow Carlos and my home directories. Feel free to add more.
[CommonLispStat.git] / linalg.lsp
blob13cbeadf356545b40f868b0dc0405f7aadf5cc76
1 ;;; -*- mode: lisp -*-
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.
5 ;;;
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.
14 ;;;
15 ;;; #3 - Use a lisp-based matrix system drop-in? (matlisp, femlisp, clem, ...?)
16 ;;;
19 ;;;; linalg -- Lisp-Stat interface to basic linear algebra routines.
20 ;;;;
21 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
22 ;;;; unrestricted use.
24 ;;;
25 ;;; Package Setup
26 ;;;
28 (in-package :cl-user)
30 (defpackage :lisp-stat-linalg
31 (:use :common-lisp
32 :cffi
33 :lisp-stat-ffi-int
34 :lisp-stat-math
35 :lisp-stat-types
36 :lisp-stat-float
37 :lisp-stat-compound-data
38 :lisp-stat-linalg-data
39 :lisp-stat-matrix)
40 (:shadowing-import-from :lisp-stat-math
41 expt + - * / ** mod rem abs 1+ 1- log exp sqrt sin cos tan
42 asin acos atan sinh cosh tanh asinh acosh atanh float random
43 truncate floor ceiling round minusp zerop plusp evenp oddp
44 < <= = /= >= > complex conjugate realpart imagpart phase
45 min max logand logior logxor lognot ffloor fceiling
46 ftruncate fround signum cis)
47 (:export chol-decomp lu-decomp lu-solve determinant inverse
48 sv-decomp qr-decomp rcondest make-rotation spline
49 kernel-dens kernel-smooth
50 fft make-sweep-matrix sweep-operator ax+y eigen
52 check-real ;; for optimize
54 covariance-matrix matrix print-matrix solve
55 backsolve eigenvalues eigenvectors accumulate cumsum combine
56 lowess))
58 (in-package #:lisp-stat-linalg)
60 ;;; CFFI Support
62 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
63 ;;;
64 ;;; Lisp Interfaces to Linear Algebra Routines
65 ;;;
66 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
68 ;;;
69 ;;; Cholesky Decomposition
70 ;;;
72 (cffi:defcfun ("ccl_chol_decomp_front" ccl-chol-decomp-front)
73 :int (x :pointer) (y :int) (z :pointer))
74 (defun chol-decomp-front (x y z)
75 (ccl-chol-decomp-front x y z))
77 ;;;;
78 ;;;; LU Decomposition
79 ;;;;
81 (cffi:defcfun ("ccl_lu_decomp_front" ccl-lu-decomp-front)
82 :int (x :pointer) (y :int) (z :pointer) (u :int) (v :pointer))
83 (defun lu-decomp-front (x y z u v)
84 (ccl-lu-decomp-front x y z u v))
86 (cffi:defcfun ("ccl_lu_solve_front" ccl-lu-solve-front)
87 :int (x :pointer) (y :int) (z :pointer) (u :pointer) (v :int))
88 (defun lu-solve-front (x y z u v)
89 (ccl-lu-solve-front x y z u v))
91 (cffi:defcfun ("ccl_lu_inverse_front" ccl-lu-inverse-front)
92 :int (x :pointer) (y :int) (z :pointer) (u :pointer) (v :int) (w :pointer))
93 (defun lu-inverse-front (x y z u v w)
94 (ccl-lu-inverse-front x y z u v w))
96 ;;;;
97 ;;;; SV Decomposition
98 ;;;;
100 (cffi:defcfun ("ccl_sv_decomp_front" ccl-sv-decomp-front)
101 :int (x :pointer) (y :int) (z :int) (u :pointer) (v :pointer))
102 (defun sv-decomp-front (x y z u v)
103 (ccl-sv-decomp-front x y z u v))
105 ;;;;
106 ;;;; QR Decomposition
107 ;;;;
109 (cffi:defcfun ("ccl_qr_decomp_front" ccl-qr-decomp-front)
110 :int (x :pointer) (y :int) (z :int) (u :pointer) (v :pointer) (w :int))
111 (defun qr-decomp-front (x y z u v w)
112 (ccl-qr-decomp-front x y z u v w))
115 ;;; Estimate of Condition Number for Lower Triangular Matrix
118 (cffi:defcfun ("ccl_rcondest_front" ccl-rcondest-front)
119 :double (x :pointer) (y :int))
120 (defun rcondest-front (x y)
121 (ccl-rcondest-front x y))
125 (defun rcondest-lisp-front (mat)
126 "Lisp-only version of rcondest."
127 (if (and (check-square-matrix mat)
128 (not (is-complex mat))) ;; Complex condition estimate not available
129 (if (= 0.0 (diag mat))
131 (let ((est (aref mat 0 0)))
132 (dotimes (j 0 n)
134 /* Set est to reciprocal of L1 matrix norm of A */
135 est = fabs(a[0][0]);
136 for (j = 1; j < n; j++) {
137 for (i = 0, temp = fabs(a[j][j]); i < j; i++)
138 temp += fabs(a[i][j]);
139 est = Max(est, temp);
141 est = 1.0 / est;
143 /* Solve A^Tx = e, selecting e as you proceed */
144 x[0] = 1.0 / a[0][0];
145 for (i = 1; i < n; i++) p[i] = a[0][i] * x[0];
146 for (j = 1; j < n; j++) {
147 /* select ej and calculate x[j] */
148 xp = ( 1.0 - p[j]) / a[j][j];
149 xm = (-1.0 - p[j]) / a[j][j];
150 temp = fabs(xp);
151 tempm = fabs(xm);
152 for (i = j + 1; i < n; i++) {
153 pm[i] = p[i] + a[j][i] * xm;
154 tempm += fabs(pm[i] / a[i][i]);
155 p[i] += a[j][i] * xp;
156 temp += fabs(p[i] / a[i][i]);
158 if (temp >= tempm) x[j] = xp;
159 else {
160 x[j] = xm;
161 for (i = j + 1; i < n; i++) p[i] = pm[i];
165 for (j = 0, xnorm = 0.0; j < n; j++) xnorm += fabs(x[j]);
166 est = est * xnorm;
167 backsolve(a, x, n, RE);
168 for (j = 0, xnorm = 0.0; j < n; j++) xnorm += fabs(x[j]);
169 if (xnorm > 0) est = est / xnorm;
174 ;;;;
175 ;;;; Make Rotation Matrix
176 ;;;;
178 (cffi:defcfun ("ccl_make_rotation_front" ccl-make-rotation-front)
179 :int (x :int) (y :pointer) (z :pointer) (u :pointer) (v :int) (w :double))
180 (defun make-rotation-front (x y z u v w)
181 (ccl-make-rotation-front x y z u v (float w 1d0)))
183 ;;;;
184 ;;;; Eigenvalues and Eigenvectors
185 ;;;;
187 (cffi:defcfun ("ccl_eigen_front" ccl-eigen-front)
188 :int (x :pointer) (y :int) (z :pointer) (u :pointer) (v :pointer))
189 (defun eigen-front (x y z u v)
190 (ccl-eigen-front x y z u v))
192 ;;;;
193 ;;;; Spline Interpolation
194 ;;;;
196 (cffi:defcfun ("ccl_range_to_rseq" ccl-range-to-rseq)
197 :int (x :int) (y :pointer) (z :int) (u :pointer))
198 (defun la-range-to-rseq (x y z u)
199 (ccl-range-to-rseq x y z u))
201 (cffi:defcfun ("ccl_spline_front" ccl-spline-front)
202 :int (x :int) (y :pointer) (z :pointer) (u :int) (v :pointer) (w :pointer) (a :pointer))
203 (defun spline-front (x y z u v w a)
204 (ccl-spline-front x y z u v w a))
206 ;;;;
207 ;;;; Kernel Density Estimators and Smoothers
208 ;;;;
210 (cffi:defcfun ("ccl_kernel_dens_front" ccl-kernel-dens-front)
211 :int (x :pointer) (y :int) (z :double) (u :pointer) (v :pointer) (w :int) (a :int))
212 (defun kernel-dens-front (x y z u v w a)
213 (ccl-kernel-dens-front x y (float z 1d0) u v w a))
215 (cffi:defcfun ("ccl_kernel_smooth_front" ccl-kernel-smooth-front)
216 :int (x :pointer) (y :pointer) (z :int) (u :double) (v :pointer) (w :pointer) (a :int) (b :int))
217 (defun kernel-smooth-front (x y z u v w a b)
218 (ccl-kernel-smooth-front x y z (float u 1d0) v w a b))
220 ;;;;
221 ;;;; Lowess Smoother Interface
222 ;;;;
224 (cffi:defcfun ("ccl_base_lowess_front" ccl-base-lowess-front)
225 :int (x :pointer) (y :pointer) (z :int) (u :double) (v :int) (w :double) (a :pointer) (b :pointer) (c :pointer))
226 (defun base-lowess-front (x y z u v w a b c)
227 (ccl-base-lowess-front x y z (float u 1d0) v (float w 1d0) a b c))
229 ;;;;
230 ;;;; FFT
231 ;;;;
233 (cffi:defcfun ("ccl_fft_front" ccl-fft-front)
234 :int (x :int) (y :pointer) (z :pointer) (u :int))
235 (defun fft-front (x y z u)
236 (ccl-fft-front x y z u))
241 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
242 ;;;;
243 ;;;; Lisp to C number conversion and checking
244 ;;;;
245 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
247 ;;;;
248 ;;;; Lisp to/from C sequence and matrix conversion and checking
249 ;;;;
251 (defun is-cons (a)
252 "FIXME:AJR this is not used anywhere?"
253 (if (consp a) 1 0))
255 (defun check-fixnum (a)
256 (if (/= 0 (la-data-mode a)) (error "not an integer sequence - ~s" a)))
258 (defun check-real (data)
259 (let ((data (compound-data-seq data)))
260 (cond
261 ((vectorp data)
262 (let ((n (length data)))
263 (declare (fixnum n))
264 (dotimes (i n)
265 (declare (fixnum i))
266 (check-one-real (aref data i)))))
267 ((consp data) (dolist (x data) (check-one-real x)))
268 (t (error "bad sequence - ~s" data)))))
270 (defun vec-assign (a i x) (setf (aref a i) x))
272 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
273 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
274 ;;;;
275 ;;;; Lisp Interfaces to Linear Algebra Routines
276 ;;;;
277 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
278 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
280 ;;; FIXME: use dpbt[f2|rf], dpbstf, dpot[f2|rf]; dpptrf, zpbstf, zpbt[f2|rf]
281 ;;; remember: factorization = decomposition, depending on training.
283 (defun chol-decomp (a &optional (maxoffl 0.0))
284 "Args: (a)
285 Modified Cholesky decomposition. A should be a square, symmetric matrix.
286 Computes lower triangular matrix L such that L L^T = A + D where D is a
287 diagonal matrix. If A is strictly positive definite D will be zero.
288 Otherwise, D is as small as possible to make A + D numerically strictly
289 positive definite. Returns a list (L (max D))."
290 (check-square-matrix a)
291 (check-real a)
292 (let* ((n (array-dimension a 0))
293 (result (make-array (list n n)))
294 (dpars (list maxoffl 0.0)))
295 (check-real dpars)
296 (let ((mat (la-data-to-matrix a +mode-re+))
297 (dp (la-data-to-vector dpars +mode-re+)))
298 (unwind-protect
299 (progn
300 (chol-decomp-front mat n dp)
301 (la-matrix-to-data mat n n +mode-re+ result)
302 (la-vector-to-data dp 2 +mode-re+ dpars))
303 (la-free-matrix mat n)
304 (la-free-vector dp)))
305 (list result (second dpars))))
308 ;;; REPLACE with
309 ;;; (matlisp:lu M)
310 ;;; i.e. result use by:
311 ;;; (setf (values (lu-out1 lu-out2 lu-out3)) (matlisp:lu my-matrix))
312 ;;; for solution, ...
313 ;;; for lu-solve:
314 ;;; (matlisp:gesv a b &opt ipivot)
316 (defun lu-decomp (a)
317 "Args: (a)
318 A is a square matrix of numbers (real or complex). Computes the LU
319 decomposition of A and returns a list of the form (LU IV D FLAG), where
320 LU is a matrix with the L part in the lower triangle, the U part in the
321 upper triangle (the diagonal entries of L are taken to be 1), IV is a vector
322 describing the row permutation used, D is 1 if the number of permutations
323 is odd, -1 if even, and FLAG is T if A is numerically singular, NIL otherwise.
324 Used bu LU-SOLVE."
325 (check-square-matrix a)
326 (let* ((n (array-dimension a 0))
327 (mode (max +mode-re+ (la-data-mode a)))
328 (result (list (make-array (list n n)) (make-array n) nil nil)))
329 (let ((mat (la-data-to-matrix a mode))
330 (iv (la-vector n +mode-in+))
331 (d (la-vector 1 +mode-re+))
332 (singular 0))
333 (unwind-protect
334 (progn
335 (setf singular (lu-decomp-front mat n iv mode d))
336 (la-matrix-to-data mat n n mode (first result))
337 (la-vector-to-data iv n +mode-in+ (second result))
338 (setf (third result) (la-get-double d 0))
339 (setf (fourth result) (if (= singular 0.0) nil t)))
340 (la-free-matrix mat n)
341 (la-free-vector iv)
342 (la-free-vector d)))
343 result))
345 (defun lu-solve (lu lb)
346 "Args: (lu b)
347 LU is the result of (LU-DECOMP A) for a square matrix A, B is a sequence.
348 Returns the solution to the equation Ax = B. Signals an error if A is
349 singular."
350 (let ((la (first lu))
351 (lidx (second lu)))
352 (check-square-matrix la)
353 (check-sequence lidx)
354 (check-sequence lb)
355 (check-fixnum lidx)
356 (let* ((n (num-rows la))
357 (result (make-sequence (if (consp lb) 'list 'vector) n))
358 (a-mode (la-data-mode la))
359 (b-mode (la-data-mode lb)))
360 (if (/= n (length lidx)) (error "index sequence is wrong length"))
361 (if (/= n (length lb)) (error "right hand side is wrong length"))
362 (let* ((mode (max +mode-re+ a-mode b-mode))
363 (a (la-data-to-matrix la mode))
364 (indx (la-data-to-vector lidx +mode-in+))
365 (b (la-data-to-vector lb mode))
366 (singular 0))
367 (unwind-protect
368 (progn
369 (setf singular (lu-solve-front a n indx b mode))
370 (la-vector-to-data b n mode result))
371 (la-free-matrix a n)
372 (la-free-vector indx)
373 (la-free-vector b))
374 (if (/= 0.0 singular) (error "matrix is (numerically) singular"))
375 result))))
377 (defun determinant (a)
378 "Args: (m)
379 Returns the determinant of the square matrix M."
380 (let* ((lu (lu-decomp a))
381 (la (first lu))
382 (n (num-rows a))
383 (d1 (third lu))
384 (d2 0.d0))
385 (declare (fixnum n))
386 (flet ((fabs (x) (float (abs x) 0.d0)))
387 (dotimes (i n (* d1 (exp d2)))
388 (declare (fixnum i))
389 (let* ((x (aref la i i))
390 (magn (fabs x)))
391 (if (= 0.0 magn) (return 0.d0))
392 (setf d1 (* d1 (/ x magn)))
393 (setf d2 (+ d2 (log magn))))))))
395 (defun inverse (a)
396 "Args: (m)
397 Returns the inverse of the the square matrix M; signals an error if M is ill
398 conditioned or singular"
399 (check-square-matrix a)
400 (let ((n (num-rows a))
401 (mode (max +mode-re+ (la-data-mode a))))
402 (declare (fixnum n))
403 (let ((result (make-array (list n n) :initial-element 0)))
404 (dotimes (i n)
405 (declare (fixnum i))
406 (setf (aref result i i) 1))
407 (let ((mat (la-data-to-matrix a mode))
408 (inv (la-data-to-matrix result mode))
409 (iv (la-vector n +mode-in+))
410 (v (la-vector n mode))
411 (singular 0))
412 (unwind-protect
413 (progn
414 (setf singular (lu-inverse-front mat n iv v mode inv))
415 (la-matrix-to-data inv n n mode result))
416 (la-free-matrix mat n)
417 (la-free-matrix inv n)
418 (la-free-vector iv)
419 (la-free-vector v))
420 (if (/= singular 0) (error "matrix is (numerically) singular"))
421 result))))
423 ;;;;
424 ;;;; SV Decomposition
425 ;;;;
427 (defun sv-decomp (a)
428 "Args: (a)
429 A is a matrix of real numbers with at least as many rows as columns.
430 Computes the singular value decomposition of A and returns a list of the form
431 (U W V FLAG) where U and V are matrices whose columns are the left and right
432 singular vectors of A and W is the sequence of singular values of A. FLAG is T
433 if the algorithm converged, NIL otherwise."
434 (check-matrix a)
435 (let* ((m (num-rows a))
436 (n (num-cols a))
437 (mode (max +mode-re+ (la-data-mode a)))
438 (result (list (make-array (list m n))
439 (make-array n)
440 (make-array (list n n))
441 nil)))
442 (if (< m n) (error "number of rows less than number of columns"))
443 (if (= mode +mode-cx+) (error "complex SVD not available yet"))
444 (let ((mat (la-data-to-matrix a mode))
445 (w (la-vector n +mode-re+))
446 (v (la-matrix n n +mode-re+))
447 (converged 0))
448 (unwind-protect
449 (progn
450 (setf converged (sv-decomp-front mat m n w v))
451 (la-matrix-to-data mat m n mode (first result))
452 (la-vector-to-data w n mode (second result))
453 (la-matrix-to-data v n n mode (third result))
454 (setf (fourth result) (if (/= 0.0 converged) t nil)))
455 (la-free-matrix mat m)
456 (la-free-vector w)
457 (la-free-matrix v n))
458 result)))
461 ;;;;
462 ;;;; QR Decomposition
463 ;;;;
465 (defun qr-decomp (a &optional pivot)
466 "Args: (a &optional pivot)
467 A is a matrix of real numbers with at least as many rows as columns. Computes
468 the QR factorization of A and returns the result in a list of the form (Q R).
469 If PIVOT is true the columns of X are first permuted to place insure the
470 absolute values of the diagonal elements of R are nonincreasing. In this case
471 the result includes a third element, a list of the indices of the columns in
472 the order in which they were used."
473 (check-matrix a)
474 (let* ((m (num-rows a))
475 (n (num-cols a))
476 (mode (max +mode-re+ (la-data-mode a)))
477 (p (if pivot 1 0))
478 (result (if pivot
479 (list (make-array (list m n))
480 (make-array (list n n))
481 (make-array n))
482 (list (make-array (list m n)) (make-array (list n n))))))
483 (if (< m n) (error "number of rows less than number of columns"))
484 (if (= mode +mode-cx+) (error "complex QR decomposition not available yet"))
485 (let ((mat (la-data-to-matrix a mode))
486 (v (la-matrix n n +mode-re+))
487 (jpvt (la-vector n +mode-in+)))
488 (unwind-protect
489 (progn
490 (qr-decomp-front mat m n v jpvt p)
491 (la-matrix-to-data mat m n mode (first result))
492 (la-matrix-to-data v n n mode (second result))
493 (if pivot (la-vector-to-data jpvt n +mode-in+ (third result))))
494 (la-free-matrix mat m)
495 (la-free-matrix v n)
496 (la-free-vector jpvt))
497 result)))
499 ;;;;
500 ;;;; Estimate of Condition Number for Lower Triangular Matrix
501 ;;;;
503 (defun rcondest (a)
504 "Args: (a)
505 Returns an estimate of the reciprocal of the L1 condition number of an upper
506 triangular matrix a."
507 (check-square-matrix a)
508 (let ((mode (max +mode-re+ (la-data-mode a)))
509 (n (num-rows a)))
510 (if (= mode +mode-cx+)
511 (error "complex condition estimate not available yet"))
512 (let ((mat (la-data-to-matrix a mode))
513 (est 0.0))
514 (unwind-protect
515 (setf est (rcondest-front mat n))
516 (la-free-matrix mat n))
517 est)))
519 ;;;;
520 ;;;; Make Rotation Matrix
521 ;;;;
523 (defun make-rotation (x y &optional alpha)
524 "Args: (x y &optional alpha)
525 Returns a rotation matrix for rotating from X to Y, or from X toward Y
526 by angle ALPHA, in radians. X and Y are sequences of the same length."
527 (check-sequence x)
528 (check-sequence y)
529 (if alpha (check-one-real alpha))
530 (let* ((n (length x))
531 (mode (max +mode-re+ (la-data-mode x) (la-data-mode y)))
532 (use-angle (if alpha 1 0))
533 (angle (if alpha (float alpha 0.0) 0.0))
534 (result (make-array (list n n))))
535 (if (/= n (length y)) (error "sequences not the same length"))
536 (if (= mode +mode-cx+) (error "complex data not supported yet"))
537 (let ((px (la-data-to-vector x +mode-re+))
538 (py (la-data-to-vector y +mode-re+))
539 (rot (la-matrix n n +mode-re+)))
540 (unwind-protect
541 (progn
542 (make-rotation-front n rot px py use-angle angle)
543 (la-matrix-to-data rot n n +mode-re+ result))
544 (la-free-vector px)
545 (la-free-vector py)
546 (la-free-matrix rot n))
547 result)))
549 ;;;;
550 ;;;; Eigenvalues and Vectors
551 ;;;;
553 (defun eigen (a)
554 "Args: (a)
555 Returns list of list of eigenvalues and list of eigenvectors of square,
556 symmetric matrix A. Third element of result is NIL if algorithm converges.
557 If the algorithm does not converge, the third element is an integer I.
558 In this case the eigenvalues 0, ..., I are not reliable."
559 (check-square-matrix a)
560 (let ((mode (max +mode-re+ (la-data-mode a)))
561 (n (num-rows a)))
562 (if (= mode +mode-cx+) (error "matrix must be real and symmetric"))
563 (let ((evals (make-array n))
564 (evecs (make-list (* n n)))
565 (pa (la-data-to-vector (compound-data-seq a) +mode-re+))
566 (w (la-vector n +mode-re+))
567 (z (la-vector (* n n) +mode-re+))
568 (fv1 (la-vector n +mode-re+))
569 (ierr 0))
570 (unwind-protect
571 (progn
572 (setf ierr (eigen-front pa n w z fv1))
573 (la-vector-to-data w n +mode-re+ evals)
574 (la-vector-to-data z (* n n) +mode-re+ evecs))
575 (la-free-vector pa)
576 (la-free-vector z)
577 (la-free-vector w)
578 (la-free-vector fv1))
579 (list (nreverse evals)
580 (nreverse (mapcar #'(lambda (x) (coerce x 'vector))
581 (split-list evecs n)))
582 (if (/= 0 ierr) (- n ierr))))))
584 ;;;;
585 ;;;; Spline Interpolation
586 ;;;;
588 (defun make-smoother-args (x y xvals)
589 (check-sequence x)
590 (check-real x)
591 (when y
592 (check-sequence y)
593 (check-real y))
594 (unless (integerp xvals)
595 (check-sequence xvals)
596 (check-real xvals))
597 (let* ((n (length x))
598 (ns (if (integerp xvals) xvals (length xvals)))
599 (result (list (make-list ns) (make-list ns))))
600 (if (and y (/= n (length y))) (error "sequences not the same length"))
601 (list x y n (if (integerp xvals) 0 1) ns xvals result)))
603 (defun get-smoother-result (args) (seventh args))
605 (defmacro with-smoother-data ((x y xvals is-reg) &rest body)
606 `(progn
607 (check-sequence ,x)
608 (check-real ,x)
609 (when ,is-reg
610 (check-sequence ,y)
611 (check-real ,y))
612 (unless (integerp ,xvals)
613 (check-sequence ,xvals)
614 (check-real ,xvals))
615 (let* ((supplied (not (integerp ,xvals)))
616 (n (length ,x))
617 (ns (if supplied (length ,xvals) ,xvals))
618 (result (list (make-list ns) (make-list ns))))
619 (if (and ,is-reg (/= n (length ,y)))
620 (error "sequences not the same length"))
621 (if (and (not supplied) (< ns 2))
622 (error "too few points for interpolation"))
623 (let* ((px (la-data-to-vector ,x +mode-re+))
624 (py (if ,is-reg (la-data-to-vector ,y +mode-re+)))
625 (pxs (if supplied
626 (la-data-to-vector ,xvals +mode-re+)
627 (la-vector ns +mode-re+)))
628 (pys (la-vector ns +mode-re+)))
629 (unless supplied (la-range-to-rseq n px ns pxs))
630 (unwind-protect
631 (progn ,@body
632 (la-vector-to-data pxs ns +mode-re+ (first result))
633 (la-vector-to-data pys ns +mode-re+ (second result)))
634 (la-free-vector px)
635 (if ,is-reg (la-free-vector py))
636 (la-free-vector pxs)
637 (la-free-vector pys))
638 result))))
640 (defun spline (x y &key (xvals 30))
641 "Args: (x y &key xvals)
642 Returns list of x and y values of natural cubic spline interpolation of (X,Y).
643 X must be strictly increasing. XVALS can be an integer, the number of equally
644 spaced points to use in the range of X, or it can be a sequence of points at
645 which to interpolate."
646 (with-smoother-data (x y xvals t)
647 (let ((work (la-vector (* 2 n) +mode-re+))
648 (error 0))
649 (unwind-protect
650 (setf error (spline-front n px py ns pxs pys work))
651 (la-free-vector work))
652 (if (/= error 0) (error "bad data for splines")))))
654 ;;;;
655 ;;;; Kernel Density Estimators and Smoothers
656 ;;;;
658 (defun kernel-type-code (type)
659 (cond ((eq type 'u) 0)
660 ((eq type 't) 1)
661 ((eq type 'g) 2)
662 (t 3)))
664 (defun kernel-dens (x &key (type 'b) (width -1.0) (xvals 30))
665 "Args: (x &key xvals width type)
666 Returns list of x and y values of kernel density estimate of X. XVALS can be an
667 integer, the number of equally spaced points to use in the range of X, or it
668 can be a sequence of points at which to interpolate. WIDTH specifies the
669 window width. TYPE specifies the lernel and should be one of the symbols G, T,
670 U or B for gaussian, triangular, uniform or bisquare. The default is B."
671 (check-one-real width)
672 (with-smoother-data (x nil xvals nil) ;; warning about deleting unreachable code is TRUE -- 2nd arg=nil!
673 (let ((code (kernel-type-code type))
674 (error 0))
675 (setf error (kernel-dens-front px n width pxs pys ns code))
676 (if (/= 0 error) (error "bad kernel density data")))))
678 (defun kernel-smooth (x y &key (type 'b) (width -1.0) (xvals 30))
679 "Args: (x y &key xvals width type)
680 Returns list of x and y values of kernel smooth of (X,Y). XVALS can be an
681 integer, the number of equally spaced points to use in the range of X, or it
682 can be a sequence of points at which to interpolate. WIDTH specifies the
683 window width. TYPE specifies the lernel and should be one of the symbols G, T,
684 U or B for Gaussian, triangular, uniform or bisquare. The default is B."
685 (check-one-real width)
686 (with-smoother-data (x y xvals t)
687 (let ((code (kernel-type-code type))
688 (error 0))
689 (kernel-smooth-front px py n width pxs pys ns code)
690 ;; if we get the Lisp version ported from C, uncomment below and
691 ;; comment above. (thanks to Carlos Ungil for the initial CFFI
692 ;; work).
693 ;;(kernel-smooth-Cport px py n width pxs pys ns code)
694 (if (/= 0 error) (error "bad kernel density data")))))
698 (defun kernel-smooth-Cport (px py n width ;;wts wds ;; see above for mismatch?
699 xs ys ns ktype)
700 "Port of kernel_smooth (Lib/kernel.c) to Lisp.
701 FIXME:kernel-smooth-Cport : This is broken.
702 Until this is fixed, we are using Luke's C code and CFFI as glue."
703 (declare (ignore width xs))
704 (cond ((< n 1) 1.0)
705 ((and (< n 2) (<= width 0)) 1.0)
706 (t (let* ((xmin (min px))
707 (xmax (max px))
708 (width (/ (- xmax xmin) (+ 1.0 (log n)))))
709 (dotimes (i (- ns 1))
710 (setf (aref ys i)
711 (let ((wsum 0.0)
712 (ysum 0.0))
713 (dotimes (j (- n 1)) )
714 ;;;possible nasty errors...
715 ;; (let*
716 ;; ((lwidth (if wds (* width (aref wds j)) width))
717 ;; (lwt (* (kernel-Cport (aref xs i) (aref px j) lwidth ktype) ;; px?
718 ;; (if wts (aref wts j) 1.0))))
719 ;; (setf wsum (+ wsum lwt))
720 ;; (setf ysum (if py (+ ysum (* lwt (aref py j)))))) ;; py? y?
722 ;;; end of errors
723 (if py
724 (if (> wsum 0.0)
725 (/ ysum wsum)
726 0.0)
727 (/ wsum n)))))
728 (values ys)))))
732 (defun kernel-Cport (x y w ktype)
733 "Port of kernel() (Lib/kernel.c) to Lisp.
734 x,y,w are doubles, type is an integer"
735 (if (<= w 0.0)
737 (let ((z (- x y)))
738 (cond ((eq ktype "B")
739 (let* ((w (* w 2.0))
740 (z (* z 0.5)))
741 (if (and (> z -0.5)
742 (< z 0.5))
743 (/ (/ (* 15.0 (* (- 1.0 (* 4 z z)) ;; k/w
744 (- 1.0 (* 4 z z)))) ;; k/w
745 8.0)
747 0)))
748 ((eq ktype "G")
749 (let* ((w (* w 0.25))
750 (z (* z 4.0))
751 (k (/ (exp (* -0.5 z z))
752 (sqrt (* 2 PI)))))
753 (/ k w)))
754 ((eq ktype "U")
755 (let* ((w (* 1.5 w))
756 (z (* z 0.75))
757 (k (if (< (abs z) 0.5)
759 0.0)))
760 (/ k w)))
761 ((eq ktype "T")
762 (cond ((and (> z -1.0)
763 (< z 0.0))
764 (+ 1.0 z)) ;; k
765 ((and (> z 0.0)
766 (< z 1.0))
767 (- 1.0 z)) ;; k
768 (t 0.0)))
769 (t (values 0.0))))))
772 ;;;;
773 ;;;; Lowess Smoother Interface
774 ;;;;
776 (defun |base-lowess| (s1 s2 f nsteps delta)
777 (check-sequence s1)
778 (check-sequence s2)
779 (check-real s1)
780 (check-real s2)
781 (check-one-real f)
782 (check-one-fixnum nsteps)
783 (check-one-real delta)
784 (let* ((n (length s1))
785 (result (make-list n)))
786 (if (/= n (length s2)) (error "sequences not the same length"))
787 (let ((x (la-data-to-vector s1 +mode-re+))
788 (y (la-data-to-vector s2 +mode-re+))
789 (ys (la-vector n +mode-re+))
790 (rw (la-vector n +mode-re+))
791 (res (la-vector n +mode-re+))
792 (error 0))
793 (unwind-protect
794 (progn
795 (setf error (base-lowess-front x y n f nsteps delta ys rw res))
796 (la-vector-to-data ys n +mode-re+ result))
797 (la-free-vector x)
798 (la-free-vector y)
799 (la-free-vector ys)
800 (la-free-vector rw)
801 (la-free-vector res))
802 (if (/= error 0) (error "bad data for lowess"))
803 result)))
806 static LVAL add_contour_point(i, j, k, l, x, y, z, v, result)
807 int i, j, k, l;
808 RVector x, y;
809 RMatrix z;
810 double v;
811 LVAL result;
813 LVAL pt;
814 double p, q;
816 if ((z[i][j] <= v && v < z[k][l]) || (z[k][l] <= v && v < z[i][j])) {
817 xlsave(pt);
818 pt = mklist(2, NIL);
819 p = (v - z[i][j]) / (z[k][l] - z[i][j]);
820 q = 1.0 - p;
821 rplaca(pt, cvflonum((FLOTYPE) (q * x[i] + p * x[k])));
822 rplaca(cdr(pt), cvflonum((FLOTYPE) (q * y[j] + p * y[l])));
823 result = cons(pt, result);
824 xlpop();
826 return(result);
829 LVAL xssurface_contour()
831 LVAL s1, s2, mat, result;
832 RVector x, y;
833 RMatrix z;
834 double v;
835 int i, j, n, m;
837 s1 = xsgetsequence();
838 s2 = xsgetsequence();
839 mat = xsgetmatrix();
840 v = makedouble(xlgetarg());
841 xllastarg();
843 n = seqlen(s1); m = seqlen(s2);
844 if (n != numrows(mat) || m != numcols(mat)) xlfail("dimensions do not match");
845 if (data_mode(s1) == CX || data_mode(s2) == CX || data_mode(mat) == CX)
846 xlfail("data must be real");
848 x = (RVector) data_to_vector(s1, RE);
849 y = (RVector) data_to_vector(s2, RE);
850 z = (RMatrix) data_to_matrix(mat, RE);
852 xlsave1(result);
853 result = NIL;
854 for (i = 0; i < n - 1; i++) {
855 for (j = 0; j < m - 1; j++) {
856 result = add_contour_point(i, j, i, j+1, x, y, z, v, result);
857 result = add_contour_point(i, j+1, i+1, j+1, x, y, z, v, result);
858 result = add_contour_point(i+1, j+1, i+1, j, x, y, z, v, result);
859 result = add_contour_point(i+1, j, i, j, x, y, z, v, result);
862 xlpop();
864 free_vector(x);
865 free_vector(y);
866 free_matrix(z, n);
868 return(result);
873 ;;; FFT
875 ;;; FIXME:ajr
876 ;;; ??replace with matlisp:fft and matlisp:ifft (the latter for inverse mapping)
878 (defun fft (x &optional inverse)
879 "Args: (x &optional inverse)
880 Returns unnormalized Fourier transform of X, or inverse transform if INVERSE
881 is true."
882 (check-sequence x)
883 (let* ((n (length x))
884 (mode (la-data-mode x))
885 (isign (if inverse -1 1))
886 (result (if (consp x) (make-list n) (make-array n))))
887 (let ((px (la-data-to-vector x +mode-cx+))
888 (work (la-vector (+ (* 4 n) 15) +mode-re+)))
889 (unwind-protect
890 (progn
891 (fft-front n px work isign)
892 (la-vector-to-data px n +mode-cx+ result))
893 (la-free-vector px)
894 (la-free-vector work))
895 result)))
898 ;;; SWEEP Operator: FIXME: use matlisp
901 (defun make-sweep-front (x y w n p mode has_w x_mean result)
902 (declare (fixnum n p mode has_w))
903 (let ((x_data nil)
904 (result_data nil)
905 (val 0.0)
906 (dxi 0.0)
907 (dyi 0.0)
908 (dv 0.0)
909 (dw 0.0)
910 (sum_w 0.0)
911 (dxik 0.0)
912 (dxjk 0.0)
913 (dyj 0.0)
914 (dx_meani 0.0)
915 (dx_meanj 0.0)
916 (dy_mean 0.0)
917 (has-w (if (/= 0 has_w) t nil))
918 (RE 1))
919 (declare (long-float val dxi dyi dv dw sum_w dxik dxjk dyj
920 dx_meani dx_meanj dy_mean))
921 ;; (declare-double val dxi dyi dv dw sum_w dxik dxjk dyj
922 ;; dx_meani dx_meanj dy_mean)
924 (if (> mode RE) (error "not supported for complex data yet"))
926 (setf x_data (compound-data-seq x))
927 (setf result_data (compound-data-seq result))
929 ;; find the mean of y
930 (setf val 0.0)
931 (setf sum_w 0.0)
932 (dotimes (i n)
933 (declare (fixnum i))
934 (setf dyi (makedouble (aref y i)))
935 (when has-w
936 (setf dw (makedouble (aref w i)))
937 (incf sum_w dw)
938 (setf dyi (* dyi dw)))
939 (incf val dyi))
940 (if (not has-w) (setf sum_w (float n 0.0)))
941 (if (<= sum_w 0.0) (error "non positive sum of weights"))
942 (setf dy_mean (/ val sum_w))
944 ;; find the column means
945 (dotimes (j p)
946 (declare (fixnum j))
947 (setf val 0.0)
948 (dotimes (i n)
949 (declare (fixnum i))
950 (setf dxi (makedouble (aref x_data (+ (* p i) j))))
951 (when has-w
952 (setf dw (makedouble (aref w i)))
953 (setf dxi (* dxi dw)))
954 (incf val dxi))
955 (setf (aref x_mean j) (/ val sum_w)))
957 ;; put 1/sum_w in topleft, means on left, minus means on top
958 (setf (aref result_data 0) (/ 1.0 sum_w))
959 (dotimes (i p)
960 (declare (fixnum i))
961 (setf dxi (makedouble (aref x_mean i)))
962 (setf (aref result_data (+ i 1)) (- dxi))
963 (setf (aref result_data (* (+ i 1) (+ p 2))) dxi))
964 (setf (aref result_data (+ p 1)) (- dy_mean))
965 (setf (aref result_data (* (+ p 1) (+ p 2))) dy_mean)
967 ;; put sums of adjusted cross products in body
968 (dotimes (i p)
969 (declare (fixnum i))
970 (dotimes (j p)
971 (declare (fixnum j))
972 (setf val 0.0)
973 (dotimes (k n)
974 (declare (fixnum k))
975 (setf dxik (makedouble (aref x_data (+ (* p k) i))))
976 (setf dxjk (makedouble (aref x_data (+ (* p k) j))))
977 (setf dx_meani (makedouble (aref x_mean i)))
978 (setf dx_meanj (makedouble (aref x_mean j)))
979 (setf dv (* (- dxik dx_meani) (- dxjk dx_meanj)))
980 (when has-w
981 (setf dw (makedouble (aref w k)))
982 (setf dv (* dv dw)))
983 (incf val dv))
984 (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ j 1))) val)
985 (setf (aref result_data (+ (* (+ j 1) (+ p 2)) (+ i 1))) val))
986 (setf val 0.0)
987 (dotimes (j n)
988 (declare (fixnum j))
989 (setf dxik (makedouble (aref x_data (+ (* p j) i))))
990 (setf dyj (makedouble (aref y j)))
991 (setf dx_meani (makedouble (aref x_mean i)))
992 (setf dv (* (- dxik dx_meani) (- dyj dy_mean)))
993 (when has-w
994 (setf dw (makedouble (aref w j)))
995 (setf dv (* dv dw)))
996 (incf val dv))
997 (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ p 1))) val)
998 (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ i 1))) val))
999 (setf val 0.0)
1000 (dotimes (j n)
1001 (declare (fixnum j))
1002 (setf dyj (makedouble (aref y j)))
1003 (setf dv (* (- dyj dy_mean) (- dyj dy_mean)))
1004 (when has-w
1005 (setf dw (makedouble (aref w j)))
1006 (setf dv (* dv dw)))
1007 (incf val dv))
1008 (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ p 1))) val)))
1010 ;;; FIXME: use matlisp
1011 (defun sweep-in-place-front (a rows cols mode k tol)
1012 "Sweep algorithm for linear regression."
1013 (declare (long-float tol))
1014 (declare (fixnum rows cols mode k))
1015 (let ((data nil)
1016 (pivot 0.0)
1017 (aij 0.0)
1018 (aik 0.0)
1019 (akj 0.0)
1020 (akk 0.0)
1021 (RE 1))
1022 (declare (long-float pivot aij aik akj akk))
1024 (if (> mode RE) (error "not supported for complex data yet"))
1025 (if (or (< k 0) (>= k rows) (>= k cols)) (error "index out of range"))
1027 (setf tol (max tol machine-epsilon))
1028 (setf data (compound-data-seq a))
1030 (setf pivot (makedouble (aref data (+ (* cols k) k))))
1032 (cond
1033 ((or (> pivot tol) (< pivot (- tol)))
1034 (dotimes (i rows)
1035 (declare (fixnum i))
1036 (dotimes (j cols)
1037 (declare (fixnum j))
1038 (when (and (/= i k) (/= j k))
1039 (setf aij (makedouble (aref data (+ (* cols i) j))))
1040 (setf aik (makedouble (aref data (+ (* cols i) k))))
1041 (setf akj (makedouble (aref data (+ (* cols k) j))))
1042 (setf aij (- aij (/ (* aik akj) pivot)))
1043 (setf (aref data (+ (* cols i) j)) aij))))
1045 (dotimes (i rows)
1046 (declare (fixnum i))
1047 (setf aik (makedouble (aref data (+ (* cols i) k))))
1048 (when (/= i k)
1049 (setf aik (/ aik pivot))
1050 (setf (aref data (+ (* cols i) k)) aik)))
1052 (dotimes (j cols)
1053 (declare (fixnum j))
1054 (setf akj (makedouble (aref data (+ (* cols k) j))))
1055 (when (/= j k)
1056 (setf akj (- (/ akj pivot)))
1057 (setf (aref data (+ (* cols k) j)) akj)))
1059 (setf akk (/ 1.0 pivot))
1060 (setf (aref data (+ (* cols k) k)) akk)
1062 (t 0))))
1064 ;; FIXME: use matlisp
1065 (defun make-sweep-matrix (x y &optional w)
1066 "Args: (x y &optional weights)
1067 X is matrix, Y and WEIGHTS are sequences. Returns the sweep matrix of the
1068 (weighted) regression of Y on X"
1069 (check-matrix x)
1070 (check-sequence y)
1071 (if w (check-sequence w))
1072 (let ((n (num-rows x))
1073 (p (num-cols x)))
1074 (if (/= n (length y)) (error "dimensions do not match"))
1075 (if (and w (/= n (length w))) (error "dimensions do not match"))
1076 (let ((mode (max (la-data-mode x)
1077 (la-data-mode x)
1078 (if w (la-data-mode w) 0)))
1079 (result (make-array (list (+ p 2) (+ p 2))))
1080 (x-mean (make-array p))
1081 (y (coerce y 'vector))
1082 (w (if w (coerce w 'vector)))
1083 (has-w (if w 1 0)))
1084 (make-sweep-front x y w n p mode has-w x-mean result)
1085 result)))
1087 (defun sweep-in-place (a k tol)
1088 (check-matrix a)
1089 (check-one-fixnum k)
1090 (check-one-real tol)
1091 (let ((rows (num-rows a))
1092 (cols (num-cols a))
1093 (mode (la-data-mode a)))
1094 (let ((swept (sweep-in-place-front a rows cols mode k tol)))
1095 (if (/= 0 swept) t nil))))
1097 (defun sweep-operator (a columns &optional tolerances)
1098 "Args: (a indices &optional tolerances)
1099 A is a matrix, INDICES a sequence of the column indices to be swept. Returns
1100 a list of the swept result and the list of the columns actually swept. (See
1101 MULTREG documentation.) If supplied, TOLERANCES should be a list of real
1102 numbers the same length as INDICES. An index will only be swept if its pivot
1103 element is larger than the corresponding element of TOLERANCES."
1104 (check-matrix a)
1105 (check-sequence columns)
1106 (if tolerances (check-sequence tolerances))
1107 (check-real a)
1108 (check-fixnum columns)
1109 (if tolerances (check-real tolerances))
1110 (do ((tol .0000001)
1111 (result (copy-array a))
1112 (swept-columns nil)
1113 (columns (coerce columns 'list) (cdr columns))
1114 (tolerances (if (consp tolerances) (coerce tolerances 'list))
1115 (if (consp tolerances) (cdr tolerances))))
1116 ((null columns) (list result swept-columns))
1117 (let ((col (first columns))
1118 (tol (if (consp tolerances) (first tolerances) tol)))
1119 (if (sweep-in-place result col tol)
1120 (setf swept-columns (cons col swept-columns))))))
1124 ;;; AX+Y
1127 ;;; matlisp:axpy
1129 (defun ax+y (a x y &optional lower)
1130 "Args (a x y &optional lower)
1131 Returns (+ (matmult A X) Y). If LOWER is not nil, A is taken to be lower
1132 triangular.
1133 This could probably be made more efficient."
1134 (check-square-matrix a)
1135 (check-sequence x)
1136 (check-sequence y)
1137 (check-real a)
1138 (check-real x)
1139 (check-real y)
1140 (let* ((n (num-rows a))
1141 (result (make-list n))
1142 (a (compound-data-seq a)))
1143 (declare (fixnum n))
1144 (if (or (/= n (length x)) (/= n (length y)))
1145 (error "dimensions do not match"))
1146 (do* ((tx (make-next-element x) (make-next-element x))
1147 (ty (make-next-element y))
1148 (tr (make-next-element result))
1149 (i 0 (+ i 1))
1150 (start 0 (+ start n))
1151 (end (if lower (+ i 1) n) (if lower (+ i 1) n)))
1152 ((<= n i) result)
1153 (declare (fixnum i start end))
1154 (let ((val (get-next-element ty i)))
1155 (dotimes (j end)
1156 (declare (fixnum j))
1157 (setf val (+ val (* (get-next-element tx j)
1158 (aref a (+ start j))))))
1159 (set-next-element tr i val)))))
1165 ;;;;
1166 ;;;; Linear Algebra Functions
1167 ;;;;
1169 (defun matrix (dim data)
1170 "Args: (dim data)
1171 returns a matrix of dimensions DIM initialized using sequence DATA
1172 in row major order."
1173 (let ((dim (coerce dim 'list))
1174 (data (coerce data 'list)))
1175 (make-array dim :initial-contents (split-list data (nth 1 dim)))))
1177 (defun flatsize (x)
1178 (length x)) ;; FIXME: defined badly!!
1180 (defun print-matrix (a &optional (stream *standard-output*))
1181 "Args: (matrix &optional stream)
1182 Prints MATRIX to STREAM in a nice form that is still machine readable"
1183 (unless (matrixp a) (error "not a matrix - ~a" a))
1184 (let ((size (min 15 (max (map-elements #'flatsize a)))))
1185 (format stream "#2a(~%")
1186 (dolist (x (row-list a))
1187 (format stream " (")
1188 (let ((n (length x)))
1189 (dotimes (i n)
1190 (let ((y (aref x i)))
1191 (cond
1192 ((integerp y) (format stream "~vd" size y))
1193 ((floatp y) (format stream "~vg" size y))
1194 (t (format stream "~va" size y))))
1195 (if (< i (- n 1)) (format stream " "))))
1196 (format stream ")~%"))
1197 (format stream " )~%")
1198 nil))
1200 (defun solve (a b)
1201 "Args: (a b)
1202 Solves A x = B using LU decomposition and backsolving. B can be a sequence
1203 or a matrix."
1204 (let ((lu (lu-decomp a)))
1205 (if (matrixp b)
1206 (apply #'bind-columns
1207 (mapcar #'(lambda (x) (lu-solve lu x)) (column-list b)))
1208 (lu-solve lu b))))
1210 (defun backsolve (a b)
1211 "Args: (a b)
1212 Solves A x = B by backsolving, assuming A is upper triangular. B must be a
1213 sequence. For use with qr-decomp."
1214 (let* ((n (length b))
1215 (sol (make-array n)))
1216 (dotimes (i n)
1217 (let* ((k (- n i 1))
1218 (val (elt b k)))
1219 (dotimes (j i)
1220 (let ((l (- n j 1)))
1221 (setq val (- val (* (aref sol l) (aref a k l))))))
1222 (setf (aref sol k) (/ val (aref a k k)))))
1223 (if (listp b) (coerce sol 'list) sol)))
1225 (defun eigenvalues (a)
1226 "Args: (a)
1227 Returns list of eigenvalues of square, symmetric matrix A"
1228 (first (eigen a)))
1230 (defun eigenvectors (a)
1231 "Args: (a)
1232 Returns list of eigenvectors of square, symmetric matrix A"
1233 (second (eigen a)))
1235 (defun accumulate (f s)
1236 "Args: (f s)
1237 Accumulates elements of sequence S using binary function F.
1238 (accumulate #'+ x) returns the cumulative sum of x."
1239 (let* ((result (list (elt s 0)))
1240 (tail result))
1241 (flet ((acc (dummy x)
1242 (rplacd tail (list (funcall f (first tail) x)))
1243 (setf tail (cdr tail))))
1244 (reduce #'acc s))
1245 (if (vectorp s) (coerce result 'vector) result)))
1247 (defun cumsum (x)
1248 "Args: (x)
1249 Returns the cumulative sum of X."
1250 (accumulate #'+ x))
1252 (defun combine (&rest args)
1253 "Args (&rest args)
1254 Returns sequence of elements of all arguments."
1255 (copy-seq (element-seq args)))
1257 (defun lowess (x y &key (f .25) (steps 2) (delta -1) sorted)
1258 "Args: (x y &key (f .25) (steps 2) delta sorted)
1259 Returns (list X YS) with YS the LOWESS fit. F is the fraction of data used for
1260 each point, STEPS is the number of robust iterations. Fits for points within
1261 DELTA of each other are interpolated linearly. If the X values setting SORTED
1262 to T speeds up the computation."
1263 (let ((x (if sorted x (sort-data x)))
1264 (y (if sorted y (select y (order x))))
1265 (delta (if (> delta 0.0) delta (/ (- (max x) (min x)) 50))))
1266 (list x)));; (|base-lowess| x y f steps delta))))