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