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