src/numerics/linalg.lisp

   1 ;;; -*- mode: lisp -*-
   2 ;;; Copyright (c) 2005--2008, 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
   7 (in-package #:lisp-matrix-stat-linalg)
   8
   9 ;;;;
  10 ;;;; Spline Interpolation
  11 ;;;;
  12
  13 (cffi:defcfun ("ccl_range_to_rseq" ccl-range-to-rseq)
  14     :int (x size-t) (y :pointer) (z size-t) (u :pointer))
  15 (defun la-range-to-rseq (x y z u)
  16   (ccl-range-to-rseq x y z u))
  17
  18 (cffi:defcfun ("ccl_spline_front" ccl-spline-front)
  19     :int (x size-t) (y :pointer) (z :pointer) (u size-t) (v :pointer) (w :pointer) (a :pointer))
  20 (defun spline-front (x y z u v w a)
  21   (ccl-spline-front x y z u v w a))
  22
  23 ;;;;
  24 ;;;; Kernel Density Estimators and Smoothers
  25 ;;;;
  26
  27 (cffi:defcfun ("ccl_kernel_dens_front" ccl-kernel-dens-front)
  28     :int (x :pointer) (y size-t) (z :double) (u :pointer) (v :pointer) (w size-t) (a :int))
  29 (defun kernel-dens-front (x y z u v w a)
  30   (ccl-kernel-dens-front x y (float z 1d0) u v w a))
  31
  32 (cffi:defcfun ("ccl_kernel_smooth_front" ccl-kernel-smooth-front)
  33     :int (x :pointer) (y :pointer) (z size-t) (u :double) (v :pointer) (w :pointer) (a size-t) (b :int))
  34 (defun kernel-smooth-front (x y z u v w a b)
  35   (ccl-kernel-smooth-front x y z (float u 1d0) v w a b))
  36
  37 ;;;;
  38 ;;;; Lowess Smoother Interface
  39 ;;;;
  40
  41 (cffi:defcfun ("ccl_base_lowess_front" ccl-base-lowess-front)
  42   :int (x :pointer) (y :pointer) (z size-t) (u :double) (v size-t) (w :double) (a :pointer) (b :pointer) (c :pointer))
  43 (defun base-lowess-front (x y z u v w a b c)
  44   (ccl-base-lowess-front x y z (float u 1d0) v (float w 1d0) a b c))
  45
  46 ;;;;
  47 ;;;; FFT
  48 ;;;;
  49
  50 (cffi:defcfun ("ccl_fft_front" ccl-fft-front)
  51     :int (x size-t) (y :pointer) (z :pointer) (u :int))
  52 (defun fft-front (x y z u)
  53   (ccl-fft-front x y z u))
  54
  55
  56
  57 ;;;;
  58 ;;;; Spline Interpolation
  59 ;;;;
  60
  61 (defun make-smoother-args (x y xvals)
  62   (check-sequence x)
  63   (check-real x)
  64   (when y
  65         (check-sequence y)
  66         (check-real y))
  67   (unless (integerp xvals)
  68           (check-sequence xvals)
  69           (check-real xvals))
  70   (let* ((n (length x))
  71          (ns (if (integerp xvals) xvals (length xvals)))
  72          (result (list (make-list ns) (make-list ns))))
  73     (if (and y (/= n (length y))) (error "sequences not the same length"))
  74     (list x y n (if (integerp xvals) 0 1) ns xvals result)))
  75
  76 (defun get-smoother-result (args) (seventh args))
  77
  78 (defmacro with-smoother-data ((x y xvals is-reg) &rest body)
  79   `(progn
  80      (check-sequence ,x)
  81      (check-real ,x)
  82      (when ,is-reg
  83            (check-sequence ,y)
  84            (check-real ,y))
  85      (unless (integerp ,xvals)
  86              (check-sequence ,xvals)
  87              (check-real ,xvals))
  88      (let* ((supplied (not (integerp ,xvals)))
  89             (n (length ,x))
  90             (ns (if supplied (length ,xvals) ,xvals))
  91             (result (list (make-list ns) (make-list ns))))
  92        (if (and ,is-reg (/= n (length ,y)))
  93            (error "sequences not the same length"))
  94        (if (and (not supplied) (< ns 2))
  95            (error "too few points for interpolation"))
  96        (let* ((px (la-data-to-vector ,x +mode-re+))
  97               (py (if ,is-reg (la-data-to-vector ,y +mode-re+)))
  98               (pxs (if supplied
  99                        (la-data-to-vector ,xvals +mode-re+)
 100                        (la-vector ns +mode-re+)))
 101               (pys (la-vector ns +mode-re+)))
 102          (unless supplied (la-range-to-rseq n px ns pxs))
 103          (unwind-protect
 104              (progn ,@body
 105                     (la-vector-to-data pxs ns +mode-re+ (first result))
 106                     (la-vector-to-data pys ns +mode-re+ (second result)))
 107            (la-free-vector px)
 108            (if ,is-reg (la-free-vector py))
 109            (la-free-vector pxs)
 110            (la-free-vector pys))
 111          result))))
 112
 113 (defun spline (x y &key (xvals 30))
 114 "Args: (x y &key xvals)
 115 Returns list of x and y values of natural cubic spline interpolation of (X,Y).
 116 X must be strictly increasing. XVALS can be an integer, the number of equally
 117 spaced points to use in the range of X, or it can be a sequence of points at
 118 which to interpolate."
 119   (with-smoother-data (x y xvals t)
 120     (let ((work (la-vector (* 2 n) +mode-re+))
 121           (error 0))
 122       (unwind-protect
 123           (setf error (spline-front n px py ns pxs pys work))
 124         (la-free-vector work))
 125       (if (/= error 0) (error "bad data for splines")))))
 126
 127 ;;;;
 128 ;;;; Kernel Density Estimators and Smoothers
 129 ;;;;
 130
 131 (defun kernel-type-code (type)
 132   (cond ((eq type 'u) 0)
 133         ((eq type 't) 1)
 134         ((eq type 'g) 2)
 135         (t 3)))
 136
 137 (defun kernel-dens (x &key (type 'b) (width -1.0) (xvals 30))
 138 "Args: (x &key xvals width type)
 139 Returns list of x and y values of kernel density estimate of X. XVALS can be an
 140 integer, the number of equally spaced points to use in the range of X, or it
 141 can be a sequence of points at which to interpolate. WIDTH specifies the
 142 window width. TYPE specifies the lernel and should be one of the symbols G, T,
 143 U or B for gaussian, triangular, uniform or bisquare. The default is B."
 144   (check-one-real width)
 145   (with-smoother-data (x nil xvals nil) ;; warning about deleting unreachable code is TRUE -- 2nd arg=nil!
 146     (let ((code (kernel-type-code type))
 147           (error 0))
 148       (setf error (kernel-dens-front px n width pxs pys ns code))
 149       (if (/= 0 error) (error "bad kernel density data")))))
 150
 151 (defun kernel-smooth (x y &key (type 'b) (width -1.0) (xvals 30))
 152 "Args: (x y &key xvals width type)
 153 Returns list of x and y values of kernel smooth of (X,Y). XVALS can be an
 154 integer, the number of equally spaced points to use in the range of X, or it
 155 can be a sequence of points at which to interpolate. WIDTH specifies the
 156 window width. TYPE specifies the lernel and should be one of the symbols G, T,
 157 U or B for Gaussian, triangular, uniform or bisquare. The default is B."
 158   (check-one-real width)
 159   (with-smoother-data (x y xvals t)
 160     (let ((code (kernel-type-code type))
 161           (error 0))
 162       (kernel-smooth-front px py n width pxs pys ns code)
 163       ;; if we get the Lisp version ported from C, uncomment below and
 164       ;; comment above.  (thanks to Carlos Ungil for the initial CFFI
 165       ;; work).
 166       ;;(kernel-smooth-Cport px py n width pxs pys ns code)
 167       (if (/= 0 error) (error "bad kernel density data")))))
 168
 169
 170 #|
 171 (defun kernel-smooth-Cport (px py n width ;;wts wds ;; see above for mismatch?
 172                             xs ys ns ktype)
 173   "Port of kernel_smooth (Lib/kernel.c) to Lisp.
 174 FIXME:kernel-smooth-Cport :  This is broken.
 175 Until this is fixed, we are using Luke's C code and CFFI as glue."
 176   (declare (ignore width xs))
 177   (cond ((< n 1) 1.0)
 178         ((and (< n 2) (<= width 0)) 1.0)
 179         (t (let* ((xmin (min px))
 180                   (xmax (max px))
 181                   (width (/ (- xmax xmin) (+ 1.0 (log n)))))
 182              (dotimes  (i (- ns 1))
 183                (setf (aref ys i)
 184                      (let ((wsum 0.0)
 185                            (ysum 0.0))
 186                        (dotimes (j (- n 1))   )
 187 ;;;possible nasty errors...
 188 ;;                     (let*
 189 ;;                           ((lwidth (if wds (* width (aref wds j)) width))
 190 ;;                            (lwt (* (kernel-Cport (aref xs i) (aref px j) lwidth ktype) ;; px?
 191 ;;                                    (if wts (aref wts j) 1.0))))
 192 ;;                         (setf wsum (+ wsum lwt))
 193 ;;                         (setf ysum (if py (+ ysum (* lwt (aref py j)))))) ;; py? y?
 194 ;;
 195 ;;; end of errors
 196                        (if py
 197                            (if (> wsum 0.0)
 198                                (/ ysum wsum)
 199                                0.0)
 200                            (/ wsum n)))))
 201              (values ys)))))
 202 |#
 203
 204
 205 (defun kernel-Cport (x y w ktype)
 206   "Port of kernel() (Lib/kernel.c) to Lisp.
 207 x,y,w are doubles, type is an integer"
 208   (if (<= w 0.0)
 209       0.0
 210       (let ((z (- x y)))
 211         (cond ((eq ktype "B")
 212                (let* ((w (* w 2.0))
 213                       (z (* z 0.5)))
 214                  (if (and (> z -0.5)
 215                           (< z 0.5))
 216                      (/ (/ (* 15.0 (* (- 1.0 (* 4 z z))  ;; k/w
 217                                       (- 1.0 (* 4 z z)))) ;; k/w
 218                            8.0)
 219                         w)
 220                      0)))
 221               ((eq ktype "G")
 222                (let* ((w (* w 0.25))
 223                       (z (* z 4.0))
 224                       (k (/ (exp (* -0.5 z z))
 225                             (sqrt (* 2 PI)))))
 226                  (/ k w)))
 227               ((eq ktype "U")
 228                (let* ((w (* 1.5 w))
 229                       (z (* z 0.75))
 230                       (k (if (< (abs z) 0.5)
 231                              1.0
 232                              0.0)))
 233                  (/ k w)))
 234               ((eq ktype "T")
 235                (cond ((and (> z -1.0)
 236                            (< z 0.0))
 237                       (+ 1.0 z))  ;; k
 238                      ((and (> z 0.0)
 239                            (< z 1.0))
 240                       (- 1.0 z))  ;; k
 241                      (t 0.0)))
 242               (t (values 0.0))))))
 243
 244
 245 ;;;;
 246 ;;;; Lowess Smoother Interface
 247 ;;;;
 248
 249 (defun |base-lowess| (s1 s2 f nsteps delta)
 250   (check-sequence s1)
 251   (check-sequence s2)
 252   (check-real s1)
 253   (check-real s2)
 254   (check-one-real f)
 255   (check-one-fixnum nsteps)
 256   (check-one-real delta)
 257   (let* ((n (length s1))
 258          (result (make-list n)))
 259     (if (/= n (length s2)) (error "sequences not the same length"))
 260     (let ((x (la-data-to-vector s1 +mode-re+))
 261           (y (la-data-to-vector s2 +mode-re+))
 262           (ys (la-vector n +mode-re+))
 263           (rw (la-vector n +mode-re+))
 264           (res (la-vector n +mode-re+))
 265           (error 0))
 266       (unwind-protect
 267           (progn
 268             (setf error (base-lowess-front x y n f nsteps delta ys rw res))
 269             (la-vector-to-data ys n +mode-re+ result))
 270         (la-free-vector x)
 271         (la-free-vector y)
 272         (la-free-vector ys)
 273         (la-free-vector rw)
 274         (la-free-vector res))
 275       (if (/= error 0) (error "bad data for lowess"))
 276       result)))
 277
 278 #|
 279 static LVAL add_contour_point(i, j, k, l, x, y, z, v, result)
 280         int i, j, k, l;
 281         RVector x, y;
 282         RMatrix z;
 283         double v;
 284         LVAL result;
 285 {
 286   LVAL pt;
 287   double p, q;
 288
 289   if ((z[i][j] <= v && v < z[k][l]) || (z[k][l] <= v && v < z[i][j])) {
 290     xlsave(pt);
 291     pt = mklist(2, NIL);
 292     p = (v - z[i][j]) / (z[k][l] - z[i][j]);
 293     q = 1.0 - p;
 294     rplaca(pt, cvflonum((FLOTYPE) (q * x[i] + p * x[k])));
 295     rplaca(cdr(pt), cvflonum((FLOTYPE) (q * y[j] + p * y[l])));
 296     result = cons(pt, result);
 297     xlpop();
 298   }
 299   return(result);
 300 }
 301
 302 LVAL xssurface_contour()
 303 {
 304   LVAL s1, s2, mat, result;
 305   RVector x, y;
 306   RMatrix z;
 307   double v;
 308   int i, j, n, m;
 309
 310   s1 = xsgetsequence();
 311   s2 = xsgetsequence();
 312   mat = xsgetmatrix();
 313   v = makedouble(xlgetarg());
 314   xllastarg();
 315
 316   n = seqlen(s1); m = seqlen(s2);
 317   if (n != numrows(mat) || m != numcols(mat)) xlfail("dimensions do not match");
 318   if (data_mode(s1) == CX || data_mode(s2) == CX || data_mode(mat) == CX)
 319     xlfail("data must be real");
 320
 321   x = (RVector) data_to_vector(s1, RE);
 322   y = (RVector) data_to_vector(s2, RE);
 323   z = (RMatrix) data_to_matrix(mat, RE);
 324
 325   xlsave1(result);
 326   result = NIL;
 327   for (i = 0; i < n - 1; i++) {
 328         for (j = 0; j < m - 1; j++) {
 329           result = add_contour_point(i, j, i, j+1, x, y, z, v, result);
 330           result = add_contour_point(i, j+1, i+1, j+1, x, y, z, v, result);
 331           result = add_contour_point(i+1, j+1, i+1, j, x, y, z, v, result);
 332           result = add_contour_point(i+1, j, i, j, x, y, z, v, result);
 333         }
 334   }
 335   xlpop();
 336
 337   free_vector(x);
 338   free_vector(y);
 339   free_matrix(z, n);
 340
 341   return(result);
 342 }
 343 |#
 344
 345 ;;;
 346 ;;; FFT
 347 ;;;
 348 ;;; FIXME:ajr
 349 ;;; ??replace with matlisp:fft  and matlisp:ifft (the latter for inverse mapping)
 350 ;;;
 351 (defun fft (x &optional inverse)
 352 "Args: (x &optional inverse)
 353 Returns unnormalized Fourier transform of X, or inverse transform if INVERSE
 354 is true."
 355   (check-sequence x)
 356   (let* ((n (length x))
 357          ;;(mode (la-data-mode x))
 358          (isign (if inverse -1 1))
 359          (result (if (consp x) (make-list n) (make-array n))))
 360     (let ((px (la-data-to-vector x +mode-cx+))
 361           (work (la-vector (+ (* 4 n) 15) +mode-re+)))
 362       (unwind-protect
 363           (progn
 364             (fft-front n px work isign)
 365             (la-vector-to-data px n +mode-cx+ result))
 366         (la-free-vector px)
 367         (la-free-vector work))
 368       result)))
 369
 370 ;;;
 371 ;;; SWEEP Operator: FIXME: use matlisp
 372 ;;;
 373
 374 (defun make-sweep-front (x y w n p mode has_w x_mean result)
 375   (declare (fixnum n p mode has_w))
 376   (let ((x_data nil)
 377         (result_data nil)
 378         (val 0.0)
 379         (dxi 0.0)
 380         (dyi 0.0)
 381         (dv 0.0)
 382         (dw 0.0)
 383         (sum_w 0.0)
 384         (dxik 0.0)
 385         (dxjk 0.0)
 386         (dyj 0.0)
 387         (dx_meani 0.0)
 388         (dx_meanj 0.0)
 389         (dy_mean 0.0)
 390         (has-w (if (/= 0 has_w) t nil))
 391         (RE 1))
 392     (declare (long-float val dxi dyi dv dw sum_w dxik dxjk dyj
 393                     dx_meani dx_meanj dy_mean)) ;; originally "declare-double" macro
 394
 395     (if (> mode RE) (error "not supported for complex data yet"))
 396
 397     (setf x_data (compound-data-seq x))
 398     (setf result_data (compound-data-seq result))
 399
 400     ;; find the mean of y
 401     (setf val 0.0)
 402     (setf sum_w 0.0)
 403     (dotimes (i n)
 404       (declare (fixnum i))
 405       (setf dyi (makedouble (aref y i)))
 406       (when has-w
 407         (setf dw (makedouble (aref w i)))
 408         (incf sum_w dw)
 409         (setf dyi (* dyi dw)))
 410       (incf val dyi))
 411     (if (not has-w) (setf sum_w (float n 0.0)))
 412     (if (<= sum_w 0.0) (error "non positive sum of weights"))
 413     (setf dy_mean (/ val sum_w))
 414
 415     ;; find the column means
 416     (dotimes (j p)
 417       (declare (fixnum j))
 418       (setf val 0.0)
 419       (dotimes (i n)
 420         (declare (fixnum i))
 421         (setf dxi (makedouble (aref x_data (+ (* p i) j))))
 422         (when has-w
 423           (setf dw (makedouble (aref w i)))
 424           (setf dxi (* dxi dw)))
 425         (incf val dxi))
 426       (setf (aref x_mean j) (/ val sum_w)))
 427
 428     ;; put 1/sum_w in topleft, means on left, minus means on top
 429     (setf (aref result_data 0) (/ 1.0 sum_w))
 430     (dotimes (i p)
 431       (declare (fixnum i))
 432       (setf dxi (makedouble (aref x_mean i)))
 433       (setf (aref result_data (+ i 1)) (- dxi))
 434       (setf (aref result_data (* (+ i 1) (+ p 2))) dxi))
 435     (setf (aref result_data (+ p 1)) (- dy_mean))
 436     (setf (aref result_data (* (+ p 1) (+ p 2))) dy_mean)
 437
 438     ;; put sums of adjusted cross products in body
 439     (dotimes (i p)
 440       (declare (fixnum i))
 441       (dotimes (j p)
 442         (declare (fixnum j))
 443         (setf val 0.0)
 444         (dotimes (k n)
 445           (declare (fixnum k))
 446           (setf dxik (makedouble (aref x_data (+ (* p k) i))))
 447           (setf dxjk (makedouble (aref x_data (+ (* p k) j))))
 448           (setf dx_meani (makedouble (aref x_mean i)))
 449           (setf dx_meanj (makedouble (aref x_mean j)))
 450           (setf dv (* (- dxik dx_meani) (- dxjk dx_meanj)))
 451           (when has-w
 452             (setf dw (makedouble (aref w k)))
 453             (setf dv (* dv dw)))
 454           (incf val dv))
 455         (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ j 1))) val)
 456         (setf (aref result_data (+ (* (+ j 1) (+ p 2)) (+ i 1))) val))
 457       (setf val 0.0)
 458       (dotimes (j n)
 459         (declare (fixnum j))
 460         (setf dxik (makedouble (aref x_data (+ (* p j) i))))
 461         (setf dyj (makedouble (aref y j)))
 462         (setf dx_meani (makedouble (aref x_mean i)))
 463         (setf dv (* (- dxik dx_meani) (- dyj dy_mean)))
 464         (when has-w
 465           (setf dw (makedouble (aref w j)))
 466           (setf dv (* dv dw)))
 467         (incf val dv))
 468       (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ p 1))) val)
 469       (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ i 1))) val))
 470     (setf val 0.0)
 471     (dotimes (j n)
 472       (declare (fixnum j))
 473       (setf dyj (makedouble (aref y j)))
 474       (setf dv (* (- dyj dy_mean) (- dyj dy_mean)))
 475       (when has-w
 476         (setf dw (makedouble (aref w j)))
 477         (setf dv (* dv dw)))
 478       (incf val dv))
 479     (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ p 1))) val)))
 480
 481 ;;; FIXME: use matlisp
 482 (defun sweep-in-place-front (a rows cols mode k tol)
 483   "Sweep algorithm for linear regression."
 484   (declare (long-float tol))
 485   (declare (fixnum rows cols mode k))
 486   (let ((data nil)
 487         (pivot 0.0)
 488         (aij 0.0)
 489         (aik 0.0)
 490         (akj 0.0)
 491         (akk 0.0)
 492         (RE 1))
 493     (declare (long-float pivot aij aik akj akk))
 494
 495     (if (> mode RE) (error "not supported for complex data yet"))
 496     (if (or (< k 0) (>= k rows) (>= k cols)) (error "index out of range"))
 497
 498     (setf tol (max tol machine-epsilon))
 499     (setf data (compound-data-seq a))
 500
 501     (setf pivot (makedouble (aref data (+ (* cols k) k))))
 502
 503     (cond
 504      ((or (> pivot tol) (< pivot (- tol)))
 505       (dotimes (i rows)
 506         (declare (fixnum i))
 507         (dotimes (j cols)
 508           (declare (fixnum j))
 509           (when (and (/= i k) (/= j k))
 510             (setf aij (makedouble (aref data (+ (* cols i) j))))
 511             (setf aik (makedouble (aref data (+ (* cols i) k))))
 512             (setf akj (makedouble (aref data (+ (* cols k) j))))
 513             (setf aij (- aij (/ (* aik akj) pivot)))
 514             (setf (aref data (+ (* cols i) j)) aij))))
 515
 516       (dotimes (i rows)
 517         (declare (fixnum i))
 518         (setf aik (makedouble (aref data (+ (* cols i) k))))
 519         (when (/= i k)
 520           (setf aik (/ aik pivot))
 521           (setf (aref data (+ (* cols i) k)) aik)))
 522
 523       (dotimes (j cols)
 524         (declare (fixnum j))
 525         (setf akj (makedouble (aref data (+ (* cols k) j))))
 526         (when (/= j k)
 527           (setf akj (- (/ akj pivot)))
 528           (setf (aref data (+ (* cols k) j)) akj)))
 529
 530       (setf akk (/ 1.0 pivot))
 531       (setf (aref data (+ (* cols k) k)) akk)
 532       1)
 533      (t 0))))
 534
 535 ;; FIXME: use matlisp
 536 (defun make-sweep-matrix (x y &optional w)
 537   "Args: (x y &optional weights)
 538  X is matrix-like, Y and WEIGHTS are vector-like. Returns the sweep matrix of the
 539  (weighted) regression of Y on X"
 540    (assert (typep x 'matrix-like))
 541    (assert (typep y 'vector-like))
 542    (if w (assert (typep w 'vector-like)))
 543    (let ((n (matrix-dimension x 0))
 544          (p (matrix-dimension x 1)))
 545      (if (/= n (length y)) (error "dimensions do not match"))
 546      (if (and w (/= n (length w))) (error "dimensions do not match"))
 547      (let ((mode (max (la-data-mode x)
 548                      (la-data-mode x)
 549                      (if w (la-data-mode w) 0)))
 550            (result (make-matrix (+ p 2) (+ p 2))))
 551        (x-mean (make-vector p))
 552        (has-w (if w 1 0))
 553        (make-sweep-front x y w n p mode has-w x-mean result)
 554        result)))
 555
 556  (defun sweep-in-place (a k tol)
 557   (assert (typep a 'matrix-like))
 558   (check-one-fixnum k)
 559   (check-one-real tol)
 560   (let ((rows (num-rows a))
 561         (cols (num-cols a))
 562         (mode (la-data-mode a)))
 563     (let ((swept (sweep-in-place-front
 564                   a
 565                   (matrix-dimensions a 0)
 566                   (matrix-dimensions a 1)
 567                   mode k tol)))
 568       (if (/= 0 swept) t nil))))
 569
 570 (defun sweep-operator (a columns &optional tolerances)
 571   "Args: (a indices &optional tolerances)
 572
 573 A is a matrix, INDICES a sequence of the column indices to be
 574 swept. Returns a list of the swept result and the list of the columns
 575 actually swept. (See MULTREG documentation.) If supplied, TOLERANCES
 576 should be a list of real numbers the same length as INDICES. An index
 577 will only be swept if its pivot element is larger than the
 578 corresponding element of TOLERANCES."
 579
 580   (check-matrix a)
 581   (if (not (typep columns 'sequence))
 582       (setf columns (list columns)))
 583   (check-sequence columns)
 584   (if tolerances
 585       (progn
 586         (if (not (typep tolerances 'sequence))
 587             (setf tolerances (list tolerances)))
 588         (check-sequence tolerances)))
 589
 590   (check-real a)
 591   (check-fixnum columns)
 592   (if tolerances (check-real tolerances))
 593   (do ((tol .0000001)
 594        (result (copy-array a))
 595        (swept-columns nil)
 596        (columns (coerce columns 'list) (cdr columns))
 597        (tolerances (if (consp tolerances) (coerce tolerances 'list))
 598                    (if (consp tolerances) (cdr tolerances))))
 599       ((null columns) (list result swept-columns))
 600       (let ((col (first columns))
 601             (tol (if (consp tolerances) (first tolerances) tol)))
 602         (if (sweep-in-place result col tol)
 603             (setf swept-columns (cons col swept-columns))))))
 604
 605
 606
 607
 608 (defun accumulate (f s)
 609 "Args: (f s)
 610 Accumulates elements of sequence S using binary function F.
 611 (accumulate #'+ x) returns the cumulative sum of x."
 612   (let* ((result (list (elt s 0)))
 613          (tail result))
 614     (flet ((acc (dummy x)
 615                 (rplacd tail (list (funcall f (first tail) x)))
 616                 (setf tail (cdr tail))))
 617       (reduce #'acc s))
 618     (if (vectorp s) (coerce result 'vector) result)))
 619
 620 (defun cumsum (x)
 621   "Args: (x)
 622 Returns the cumulative sum of X."
 623   (accumulate #'+ x))
 624
 625 (defun combine (&rest args)
 626   "Args (&rest args)
 627 Returns sequence of elements of all arguments."
 628   (copy-seq (element-seq args)))
 629
 630 (defun lowess (x y &key (f .25) (steps 2) (delta -1) sorted)
 631 "Args: (x y &key (f .25) (steps 2) delta sorted)
 632 Returns (list X YS) with YS the LOWESS fit. F is the fraction of data used for
 633 each point, STEPS is the number of robust iterations. Fits for points within
 634 DELTA of each other are interpolated linearly. If the X values setting SORTED
 635 to T speeds up the computation."
 636   (let ((x (if sorted x (sort-data x)))
 637         (y (if sorted y (select y (order x))))
 638         (delta (if (> delta 0.0) delta (/ (- (max x) (min x)) 50))))
 639     (list x y delta f steps)));; (|base-lowess| x y f steps delta))))