linalg.lsp

   1 ;;;; linalg -- Lisp-Stat interface to basic linear algebra routines.
   2 ;;;;
   3 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
   4 ;;;; unrestricted use.
   5
   6 (provide "linalg")
   7
   8 ;;;;
   9 ;;;; Package Setup
  10 ;;;;
  11
  12 #+:CLtL2
  13 (in-package lisp-stat-basics)
  14 #-:CLtL2
  15 (in-package 'lisp-stat-basics)
  16
  17 (export '(chol-decomp lu-decomp lu-solve determinant inverse sv-decomp
  18           qr-decomp rcondest make-rotation spline kernel-dens kernel-smooth
  19           fft make-sweep-matrix sweep-operator ax+y numgrad numhess
  20           split-list eigen))
  21
  22 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  23 ;;;;
  24 ;;;;                 Lisp to C number conversion and checking
  25 ;;;;
  26 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  27
  28 ;;;;
  29 ;;;; Lisp to/from C sequence and matrix conversion and checking
  30 ;;;;
  31
  32 (defun is-cons (a) (if (consp a) 1 0))
  33
  34 (defun check-fixnum (a)
  35   (if (/= 0 (la-data-mode a)) (error "not an integer sequence - ~s" a)))
  36
  37 (defun check-real (data)
  38   (let ((data (compound-data-seq data)))
  39     (cond
  40      ((vectorp data)
  41       (let ((n (length data)))
  42         (declare (fixnum n))
  43         (dotimes (i n)
  44                  (declare (fixnum i))
  45                  (check-one-real (aref data i)))))
  46      ((consp data) (dolist (x data) (check-one-real x)))
  47      (t (error "bad sequence - ~s" data)))))
  48
  49 (defun vec-assign (a i x) (setf (aref a i) x))
  50
  51 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  52 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  53 ;;;;
  54 ;;;;             Lisp Interfaces to Linear Algebra Routines
  55 ;;;;
  56 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  57 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  58
  59 ;;;;
  60 ;;;; Cholesky Decomposition
  61 ;;;;
  62
  63 (defun chol-decomp (a &optional (maxoffl 0.0))
  64 "Args: (a)
  65 Modified Cholesky decomposition. A should be a square, symmetric matrix.
  66 Computes lower triangular matrix L such that L L^T = A + D where D is a
  67 diagonal matrix. If A is strictly positive definite D will be zero.
  68 Otherwise, D is as small as possible to make A + D numerically strictly
  69 positive definite. Returns a list (L (max D))."
  70   (check-square-matrix a)
  71   (check-real a)
  72   (let* ((n (array-dimension a 0))
  73          (result (make-array (list n n)))
  74          (dpars (list maxoffl 0.0)))
  75     (check-real dpars)
  76     (let ((mat (la-data-to-matrix a mode-re))
  77           (dp (la-data-to-vector dpars mode-re)))
  78       (unwind-protect
  79           (progn
  80             (chol-decomp-front mat n dp)
  81             (la-matrix-to-data mat n n mode-re result)
  82             (la-vector-to-data dp 2 mode-re dpars))
  83         (la-free-matrix mat n)
  84         (la-free-vector dp)))
  85     (list result (second dpars))))
  86
  87 ;;;;
  88 ;;;; LU Decomposition
  89 ;;;;
  90
  91 (defun lu-decomp (a)
  92 "Args: (a)
  93 A is a square matrix of numbers (real or complex). Computes the LU
  94 decomposition of A and returns a list of the form (LU IV D FLAG), where
  95 LU is a matrix with the L part in the lower triangle, the U part in the
  96 upper triangle (the diagonal entries of L are taken to be 1), IV is a vector
  97 describing the row permutation used, D is 1 if the number of permutations
  98 is odd, -1 if even, and FLAG is T if A is numerically singular, NIL otherwise.
  99 Used bu LU-SOLVE."
 100   (check-square-matrix a)
 101   (let* ((n (array-dimension a 0))
 102          (mode (max mode-re (la-data-mode a)))
 103          (result (list (make-array (list n n)) (make-array n) nil nil)))
 104     (let ((mat (la-data-to-matrix a mode))
 105           (iv (la-vector n mode-in))
 106           (d (la-vector 1 mode-re))
 107           (singular 0))
 108       (unwind-protect
 109           (progn
 110             (setf singular (lu-decomp-front mat n iv mode d))
 111             (la-matrix-to-data mat n n mode (first result))
 112             (la-vector-to-data iv n mode-in (second result))
 113             (setf (third result) (la-get-double d 0))
 114             (setf (fourth result) (if (= singular 0) nil t)))
 115         (la-free-matrix mat n)
 116         (la-free-vector iv)
 117         (la-free-vector d)))
 118     result))
 119
 120 (defun lu-solve (lu lb)
 121 "Args: (lu b)
 122 LU is the result of (LU-DECOMP A) for a square matrix A, B is a sequence.
 123 Returns the solution to the equation Ax = B. Signals an error if A is
 124 singular."
 125   (let ((la (first lu))
 126         (lidx (second lu)))
 127     (check-square-matrix la)
 128     (check-sequence lidx)
 129     (check-sequence lb)
 130     (check-fixnum lidx)
 131     (let* ((n (num-rows la))
 132            (result (make-sequence (if (consp lb) 'list 'vector) n))
 133            (a-mode (la-data-mode la))
 134            (b-mode (la-data-mode lb)))
 135       (if (/= n (length lidx)) (error "index sequence is wrong length"))
 136       (if (/= n (length lb)) (error "right hand side is wrong length"))
 137       (let* ((mode (max mode-re a-mode b-mode))
 138              (a (la-data-to-matrix la mode))
 139              (indx (la-data-to-vector lidx mode-in))
 140              (b (la-data-to-vector lb mode))
 141              (singular 0))
 142         (unwind-protect
 143             (progn
 144               (setf singular (lu-solve-front a n indx b mode))
 145               (la-vector-to-data b n mode result))
 146           (la-free-matrix a n)
 147           (la-free-vector indx)
 148           (la-free-vector b))
 149         (if (/= 0 singular) (error "matrix is (numerically) singular"))
 150         result))))
 151
 152 (defun determinant (a)
 153 "Args: (m)
 154 Returns the determinant of the square matrix M."
 155   (let* ((lu (lu-decomp a))
 156          (la (first lu))
 157          (n (num-rows a))
 158          (d1 (third lu))
 159          (d2 0.d0))
 160     (declare (fixnum n))
 161     (flet ((fabs (x) (float (abs x) 0.d0)))
 162       (dotimes (i n (* d1 (exp d2)))
 163         (declare (fixnum i))
 164         (let* ((x (aref la i i))
 165                (magn (fabs x)))
 166           (if (= 0.0 magn) (return 0.d0))
 167           (setf d1 (* d1 (/ x magn)))
 168           (setf d2 (+ d2 (log magn))))))))
 169
 170 (defun inverse (a)
 171 "Args: (m)
 172 Returns the inverse of the the square matrix M; signals an error if M is ill
 173 conditioned or singular"
 174   (check-square-matrix a)
 175   (let ((n (num-rows a))
 176         (mode (max mode-re (la-data-mode a))))
 177     (declare (fixnum n))
 178     (let ((result (make-array (list n n) :initial-element 0)))
 179       (dotimes (i n)
 180                (declare (fixnum i))
 181                (setf (aref result i i) 1))
 182       (let ((mat (la-data-to-matrix a mode))
 183             (inv (la-data-to-matrix result mode))
 184             (iv (la-vector n mode-in))
 185             (v (la-vector n mode))
 186             (singular 0))
 187         (unwind-protect
 188             (progn
 189               (setf singular (lu-inverse-front mat n iv v mode inv))
 190               (la-matrix-to-data inv n n mode result))
 191           (la-free-matrix mat n)
 192           (la-free-matrix inv n)
 193           (la-free-vector iv)
 194           (la-free-vector v))
 195         (if (/= singular 0) (error "matrix is (numerically) singular"))
 196         result))))
 197
 198 ;;;;
 199 ;;;; SV Decomposition
 200 ;;;;
 201
 202 (defun sv-decomp (a)
 203 "Args: (a)
 204 A is a matrix of real numbers with at least as many rows as columns.
 205 Computes the singular value decomposition of A and returns a list of the form
 206 (U W V FLAG) where U and V are matrices whose columns are the left and right
 207 singular vectors of A and W is the sequence of singular values of A. FLAG is T
 208 if the algorithm converged, NIL otherwise."
 209   (check-matrix a)
 210   (let* ((m (num-rows a))
 211          (n (num-cols a))
 212          (mode (max mode-re (la-data-mode a)))
 213          (result (list (make-array (list m n))
 214                        (make-array n)
 215                        (make-array (list n n))
 216                        nil)))
 217     (if (< m n) (error "number of rows less than number of columns"))
 218     (if (= mode mode-cx) (error "complex SVD not available yet"))
 219     (let ((mat (la-data-to-matrix a mode))
 220           (w (la-vector n mode-re))
 221           (v (la-matrix n n mode-re))
 222           (converged 0))
 223       (unwind-protect
 224           (progn
 225             (setf converged (sv-decomp-front mat m n w v))
 226             (la-matrix-to-data mat m n mode (first result))
 227             (la-vector-to-data w n mode (second result))
 228             (la-matrix-to-data v n n mode (third result))
 229             (setf (fourth result) (if (/= 0 converged) t nil)))
 230         (la-free-matrix mat m)
 231         (la-free-vector w)
 232         (la-free-matrix v n))
 233       result)))
 234
 235
 236 ;;;;
 237 ;;;; QR Decomposition
 238 ;;;;
 239
 240 (defun qr-decomp (a &optional pivot)
 241 "Args: (a &optional pivot)
 242 A is a matrix of real numbers with at least as many rows as columns. Computes
 243 the QR factorization of A and returns the result in a list of the form (Q R).
 244 If PIVOT is true the columns of X are first permuted to place insure the
 245 absolute values of the diagonal elements of R are nonincreasing. In this case
 246 the result includes a third element, a list of the indices of the columns in
 247 the order in which they were used."
 248   (check-matrix a)
 249   (let* ((m (num-rows a))
 250          (n (num-cols a))
 251          (mode (max mode-re (la-data-mode a)))
 252          (p (if pivot 1 0))
 253          (result (if pivot
 254                      (list (make-array (list m n))
 255                            (make-array (list n n))
 256                            (make-array n))
 257                      (list (make-array (list m n)) (make-array (list n n))))))
 258     (if (< m n) (error "number of rows less than number of columns"))
 259     (if (= mode mode-cx) (error "complex QR decomposition not available yet"))
 260     (let ((mat (la-data-to-matrix a mode))
 261           (v (la-matrix n n mode-re))
 262           (jpvt (la-vector n mode-in)))
 263       (unwind-protect
 264           (progn
 265             (qr-decomp-front mat m n v jpvt p)
 266             (la-matrix-to-data mat m n mode (first result))
 267             (la-matrix-to-data v n n mode (second result))
 268             (if pivot (la-vector-to-data jpvt n mode-in (third result))))
 269         (la-free-matrix mat m)
 270         (la-free-matrix v n)
 271         (la-free-vector jpvt))
 272       result)))
 273
 274 ;;;;
 275 ;;;; Estimate of Condition Number for Lower Triangular Matrix
 276 ;;;;
 277
 278 (defun rcondest (a)
 279 "Args: (a)
 280 Returns an estimate of the reciprocal of the L1 condition number of an upper
 281 triangular matrix a."
 282   (check-square-matrix a)
 283   (let ((mode (max mode-re (la-data-mode a)))
 284         (n (num-rows a)))
 285     (if (= mode mode-cx)
 286         (error "complex condition estimate not available yet"))
 287     (let ((mat (la-data-to-matrix a mode))
 288           (est 0.0))
 289       (unwind-protect
 290           (setf est (rcondest-front mat n))
 291         (la-free-matrix mat n))
 292       est)))
 293
 294 ;;;;
 295 ;;;; Make Rotation Matrix
 296 ;;;;
 297
 298 (defun make-rotation (x y &optional alpha)
 299 "Args: (x y &optional alpha)
 300 Returns a rotation matrix for rotating from X to Y, or from X toward Y
 301 by angle ALPHA, in radians. X and Y are sequences of the same length."
 302   (check-sequence x)
 303   (check-sequence y)
 304   (if alpha (check-one-real alpha))
 305   (let* ((n (length x))
 306          (mode (max mode-re (la-data-mode x) (la-data-mode y)))
 307          (use-angle (if alpha 1 0))
 308          (angle (if alpha (float alpha 0.0) 0.0))
 309          (result (make-array (list n n))))
 310     (if (/= n (length y)) (error "sequences not the same length"))
 311     (if (= mode mode-cx) (error "complex data not supported yet"))
 312     (let ((px (la-data-to-vector x mode-re))
 313           (py (la-data-to-vector y mode-re))
 314           (rot (la-matrix n n mode-re)))
 315       (unwind-protect
 316           (progn
 317             (make-rotation-front n rot px py use-angle angle)
 318             (la-matrix-to-data rot n n mode-re result))
 319         (la-free-vector px)
 320         (la-free-vector py)
 321         (la-free-matrix rot n))
 322       result)))
 323
 324 ;;;;
 325 ;;;; Eigenvalues and Vectors
 326 ;;;;
 327
 328 (defun eigen (a)
 329 "Args: (a)
 330 Returns list of list of eigenvalues and list of eigenvectors of square,
 331 symmetric matrix A. Third element of result is NIL if algorithm converges.
 332 If the algorithm does not converge, the third element is an integer I.
 333 In this case the eigenvalues 0, ..., I are not reliable."
 334   (check-square-matrix a)
 335   (let ((mode (max mode-re (la-data-mode a)))
 336         (n (num-rows a)))
 337     (if (= mode mode-cx) (error "matrix must be real and symmetric"))
 338     (let ((evals (make-array n))
 339           (evecs (make-list (* n n)))
 340           (pa (la-data-to-vector (compound-data-seq a) mode-re))
 341           (w (la-vector n mode-re))
 342           (z (la-vector (* n n) mode-re))
 343           (fv1 (la-vector n mode-re))
 344           (ierr 0))
 345       (unwind-protect
 346           (progn
 347             (setf ierr (eigen-front pa n w z fv1))
 348             (la-vector-to-data w n mode-re evals)
 349             (la-vector-to-data z (* n n) mode-re evecs))
 350         (la-free-vector pa)
 351         (la-free-vector z)
 352         (la-free-vector w)
 353         (la-free-vector fv1))
 354       (list (nreverse evals)
 355             (nreverse (mapcar #'(lambda (x) (coerce x 'vector))
 356                               (split-list evecs n)))
 357             (if (/= 0 ierr) (- n ierr))))))
 358
 359 ;;;;
 360 ;;;; Spline Interpolation
 361 ;;;;
 362
 363 (defun make-smoother-args (x y xvals)
 364   (check-sequence x)
 365   (check-real x)
 366   (when y
 367         (check-sequence y)
 368         (check-real y))
 369   (unless (integerp xvals)
 370           (check-sequence xvals)
 371           (check-real xvals))
 372   (let* ((n (length x))
 373          (ns (if (integerp xvals) xvals (length xvals)))
 374          (result (list (make-list ns) (make-list ns))))
 375     (if (and y (/= n (length y))) (error "sequences not the same length"))
 376     (list x y n (if (integerp xvals) 0 1) ns xvals result)))
 377
 378 (defun get-smoother-result (args) (seventh args))
 379
 380 (defmacro with-smoother-data ((x y xvals is-reg) &rest body)
 381   `(progn
 382      (check-sequence ,x)
 383      (check-real ,x)
 384      (when ,is-reg
 385            (check-sequence ,y)
 386            (check-real ,y))
 387      (unless (integerp ,xvals)
 388              (check-sequence ,xvals)
 389              (check-real ,xvals))
 390      (let* ((supplied (not (integerp ,xvals)))
 391             (n (length ,x))
 392             (ns (if supplied (length ,xvals) ,xvals))
 393             (result (list (make-list ns) (make-list ns))))
 394        (if (and ,is-reg (/= n (length ,y)))
 395            (error "sequences not the same length"))
 396        (if (and (not supplied) (< ns 2))
 397            (error "too few points for interpolation"))
 398        (let* ((px (la-data-to-vector ,x mode-re))
 399               (py (if ,is-reg (la-data-to-vector ,y mode-re)))
 400               (pxs (if supplied
 401                        (la-data-to-vector ,xvals mode-re)
 402                        (la-vector ns mode-re)))
 403               (pys (la-vector ns mode-re)))
 404          (unless supplied (la-range-to-rseq n px ns pxs))
 405          (unwind-protect
 406              (progn ,@body
 407                     (la-vector-to-data pxs ns mode-re (first result))
 408                     (la-vector-to-data pys ns mode-re (second result)))
 409            (la-free-vector px)
 410            (if ,is-reg (la-free-vector py))
 411            (la-free-vector pxs)
 412            (la-free-vector pys))
 413          result))))
 414
 415 (defun spline (x y &key (xvals 30))
 416 "Args: (x y &key xvals)
 417 Returns list of x and y values of natural cubic spline interpolation of (X,Y).
 418 X must be strictly increasing. XVALS can be an integer, the number of equally
 419 spaced points to use in the range of X, or it can be a sequence of points at
 420 which to interpolate."
 421   (with-smoother-data (x y xvals t)
 422     (let ((work (la-vector (* 2 n) mode-re))
 423           (error 0))
 424       (unwind-protect
 425           (setf error (spline-front n px py ns pxs pys work))
 426         (la-free-vector work))
 427       (if (/= error 0) (error "bad data for splines")))))
 428
 429 ;;;;
 430 ;;;; Kernel Density Estimators and Smoothers
 431 ;;;;
 432
 433 (defun kernel-type-code (type)
 434   (cond ((eq type 'u) 0)
 435         ((eq type 't) 1)
 436         ((eq type 'g) 2)
 437         (t 3)))
 438
 439 (defun kernel-dens (x &key (type 'b) (width -1.0) (xvals 30))
 440 "Args: (x &key xvals width type)
 441 Returns list of x and y values of kernel density estimate of X. XVALS can be an
 442 integer, the number of equally spaced points to use in the range of X, or it
 443 can be a sequence of points at which to interpolate. WIDTH specifies the
 444 window width. TYPE specifies the lernel and should be one of the symbols G, T,
 445 U or B for gaussian, triangular, uniform or bisquare. The default is B."
 446   (check-one-real width)
 447   (with-smoother-data (x nil xvals nil)
 448     (let ((code (kernel-type-code type))
 449           (error 0))
 450       (setf error (kernel-dens-front px n width pxs pys ns code))
 451       (if (/= 0 error) (error "bad kernel density data")))))
 452
 453 (defun kernel-smooth (x y &key (type 'b) (width -1.0) (xvals 30))
 454 "Args: (x y &key xvals width type)
 455 Returns list of x and y values of kernel smooth of (X,Y). XVALS can be an
 456 integer, the number of equally spaced points to use in the range of X, or it
 457 can be a sequence of points at which to interpolate. WIDTH specifies the
 458 window width. TYPE specifies the lernel and should be one of the symbols G, T,
 459 U or B for Gaussian, triangular, uniform or bisquare. The default is B."
 460   (check-one-real width)
 461   (with-smoother-data (x y xvals t)
 462     (let ((code (kernel-type-code type))
 463           (error 0))
 464       (kernel-smooth-front px py n width pxs pys ns code)
 465       (if (/= 0 error) (error "bad kernel density data")))))
 466
 467 ;;;;
 468 ;;;; Lowess Smoother Interface
 469 ;;;;
 470
 471 (defun |base-lowess| (s1 s2 f nsteps delta)
 472   (check-sequence s1)
 473   (check-sequence s2)
 474   (check-real s1)
 475   (check-real s2)
 476   (check-one-real f)
 477   (check-one-fixnum nsteps)
 478   (check-one-real delta)
 479   (let* ((n (length s1))
 480          (result (make-list n)))
 481     (if (/= n (length s2)) (error "sequences not the same length"))
 482     (let ((x (la-data-to-vector s1 mode-re))
 483           (y (la-data-to-vector s2 mode-re))
 484           (ys (la-vector n mode-re))
 485           (rw (la-vector n mode-re))
 486           (res (la-vector n mode-re))
 487           (error 0))
 488       (unwind-protect
 489           (progn
 490             (setf error (base-lowess-front x y n f nsteps delta ys rw res))
 491             (la-vector-to-data ys n mode-re result))
 492         (la-free-vector x)
 493         (la-free-vector y)
 494         (la-free-vector ys)
 495         (la-free-vector rw)
 496         (la-free-vector res))
 497       (if (/= error 0) (error "bad data for lowess"))
 498       result)))
 499
 500 #|
 501 static LVAL add_contour_point(i, j, k, l, x, y, z, v, result)
 502         int i, j, k, l;
 503         RVector x, y;
 504         RMatrix z;
 505         double v;
 506         LVAL result;
 507 {
 508   LVAL pt;
 509   double p, q;
 510
 511   if ((z[i][j] <= v && v < z[k][l]) || (z[k][l] <= v && v < z[i][j])) {
 512     xlsave(pt);
 513     pt = mklist(2, NIL);
 514     p = (v - z[i][j]) / (z[k][l] - z[i][j]);
 515     q = 1.0 - p;
 516     rplaca(pt, cvflonum((FLOTYPE) (q * x[i] + p * x[k])));
 517     rplaca(cdr(pt), cvflonum((FLOTYPE) (q * y[j] + p * y[l])));
 518     result = cons(pt, result);
 519     xlpop();
 520   }
 521   return(result);
 522 }
 523
 524 LVAL xssurface_contour()
 525 {
 526   LVAL s1, s2, mat, result;
 527   RVector x, y;
 528   RMatrix z;
 529   double v;
 530   int i, j, n, m;
 531
 532   s1 = xsgetsequence();
 533   s2 = xsgetsequence();
 534   mat = xsgetmatrix();
 535   v = makedouble(xlgetarg());
 536   xllastarg();
 537
 538   n = seqlen(s1); m = seqlen(s2);
 539   if (n != numrows(mat) || m != numcols(mat)) xlfail("dimensions do not match");
 540   if (data_mode(s1) == CX || data_mode(s2) == CX || data_mode(mat) == CX)
 541     xlfail("data must be real");
 542
 543   x = (RVector) data_to_vector(s1, RE);
 544   y = (RVector) data_to_vector(s2, RE);
 545   z = (RMatrix) data_to_matrix(mat, RE);
 546
 547   xlsave1(result);
 548   result = NIL;
 549   for (i = 0; i < n - 1; i++) {
 550         for (j = 0; j < m - 1; j++) {
 551           result = add_contour_point(i, j, i, j+1, x, y, z, v, result);
 552           result = add_contour_point(i, j+1, i+1, j+1, x, y, z, v, result);
 553           result = add_contour_point(i+1, j+1, i+1, j, x, y, z, v, result);
 554           result = add_contour_point(i+1, j, i, j, x, y, z, v, result);
 555         }
 556   }
 557   xlpop();
 558
 559   free_vector(x);
 560   free_vector(y);
 561   free_matrix(z, n);
 562
 563   return(result);
 564 }
 565 |#
 566
 567 ;;;;
 568 ;;;; FFT
 569 ;;;;
 570
 571 (defun fft (x &optional inverse)
 572 "Args: (x &optional inverse)
 573 Returns unnormalized Fourier transform of X, or inverse transform if INVERSE
 574 is true."
 575   (check-sequence x)
 576   (let* ((n (length x))
 577          (mode (la-data-mode x))
 578          (isign (if inverse -1 1))
 579          (result (if (consp x) (make-list n) (make-array n))))
 580     (let ((px (la-data-to-vector x mode-cx))
 581           (work (la-vector (+ (* 4 n) 15) mode-re)))
 582       (unwind-protect
 583           (progn
 584             (fft-front n px work isign)
 585             (la-vector-to-data px n mode-cx result))
 586         (la-free-vector px)
 587         (la-free-vector work))
 588       result)))
 589
 590 ;;;;
 591 ;;;; SWEEP Operator
 592 ;;;;
 593
 594 (defun make-sweep-front (x y w n p mode has_w x_mean result)
 595   (declare (fixnum n p mode has_w))
 596   (let ((x_data nil)
 597         (result_data nil)
 598         (val 0.0)
 599         (dxi 0.0)
 600         (dyi 0.0)
 601         (dv 0.0)
 602         (dw 0.0)
 603         (sum_w 0.0)
 604         (dxik 0.0)
 605         (dxjk 0.0)
 606         (dyj 0.0)
 607         (dx_meani 0.0)
 608         (dx_meanj 0.0)
 609         (dy_mean 0.0)
 610         (has-w (if (/= 0 has_w) t nil))
 611         (RE 1))
 612     (declare-double val dxi dyi dv dw sum_w dxik dxjk dyj
 613                     dx_meani dx_meanj dy_mean)
 614
 615     (if (> mode RE) (error "not supported for complex data yet"))
 616
 617     (setf x_data (compound-data-seq x))
 618     (setf result_data (compound-data-seq result))
 619
 620     ;; find the mean of y
 621     (setf val 0.0)
 622     (setf sum_w 0.0)
 623     (dotimes (i n)
 624       (declare (fixnum i))
 625       (setf dyi (makedouble (aref y i)))
 626       (when has-w
 627         (setf dw (makedouble (aref w i)))
 628         (incf sum_w dw)
 629         (setf dyi (* dyi dw)))
 630       (incf val dyi))
 631     (if (not has-w) (setf sum_w (float n 0.0)))
 632     (if (<= sum_w 0.0) (error "non positive sum of weights"))
 633     (setf dy_mean (/ val sum_w))
 634
 635     ;; find the column means
 636     (dotimes (j p)
 637       (declare (fixnum j))
 638       (setf val 0.0)
 639       (dotimes (i n)
 640         (declare (fixnum i))
 641         (setf dxi (makedouble (aref x_data (+ (* p i) j))))
 642         (when has-w
 643           (setf dw (makedouble (aref w i)))
 644           (setf dxi (* dxi dw)))
 645         (incf val dxi))
 646       (setf (aref x_mean j) (/ val sum_w)))
 647
 648     ;; put 1/sum_w in topleft, means on left, minus means on top
 649     (setf (aref result_data 0) (/ 1.0 sum_w))
 650     (dotimes (i p)
 651       (declare (fixnum i))
 652       (setf dxi (makedouble (aref x_mean i)))
 653       (setf (aref result_data (+ i 1)) (- dxi))
 654       (setf (aref result_data (* (+ i 1) (+ p 2))) dxi))
 655     (setf (aref result_data (+ p 1)) (- dy_mean))
 656     (setf (aref result_data (* (+ p 1) (+ p 2))) dy_mean)
 657
 658     ;; put sums of adjusted cross products in body
 659     (dotimes (i p)
 660       (declare (fixnum i))
 661       (dotimes (j p)
 662         (declare (fixnum j))
 663         (setf val 0.0)
 664         (dotimes (k n)
 665           (declare (fixnum k))
 666           (setf dxik (makedouble (aref x_data (+ (* p k) i))))
 667           (setf dxjk (makedouble (aref x_data (+ (* p k) j))))
 668           (setf dx_meani (makedouble (aref x_mean i)))
 669           (setf dx_meanj (makedouble (aref x_mean j)))
 670           (setf dv (* (- dxik dx_meani) (- dxjk dx_meanj)))
 671           (when has-w
 672             (setf dw (makedouble (aref w k)))
 673             (setf dv (* dv dw)))
 674           (incf val dv))
 675         (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ j 1))) val)
 676         (setf (aref result_data (+ (* (+ j 1) (+ p 2)) (+ i 1))) val))
 677       (setf val 0.0)
 678       (dotimes (j n)
 679         (declare (fixnum j))
 680         (setf dxik (makedouble (aref x_data (+ (* p j) i))))
 681         (setf dyj (makedouble (aref y j)))
 682         (setf dx_meani (makedouble (aref x_mean i)))
 683         (setf dv (* (- dxik dx_meani) (- dyj dy_mean)))
 684         (when has-w
 685           (setf dw (makedouble (aref w j)))
 686           (setf dv (* dv dw)))
 687         (incf val dv))
 688       (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ p 1))) val)
 689       (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ i 1))) val))
 690     (setf val 0.0)
 691     (dotimes (j n)
 692       (declare (fixnum j))
 693       (setf dyj (makedouble (aref y j)))
 694       (setf dv (* (- dyj dy_mean) (- dyj dy_mean)))
 695       (when has-w
 696         (setf dw (makedouble (aref w j)))
 697         (setf dv (* dv dw)))
 698       (incf val dv))
 699     (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ p 1))) val)))
 700
 701 (defun sweep-in-place-front (a rows cols mode k tol)
 702   (declare-double tol)
 703   (declare (fixnum rows cols mode k))
 704   (let ((data nil)
 705         (pivot 0.0)
 706         (aij 0.0)
 707         (aik 0.0)
 708         (akj 0.0)
 709         (akk 0.0)
 710         (RE 1))
 711     (declare-double pivot aij aik akj akk)
 712
 713     (if (> mode RE) (error "not supported for complex data yet"))
 714     (if (or (< k 0) (>= k rows) (>= k cols)) (error "index out of range"))
 715
 716     (setf tol (max tol machine-epsilon))
 717     (setf data (compound-data-seq a))
 718
 719     (setf pivot (makedouble (aref data (+ (* cols k) k))))
 720
 721     (cond
 722      ((or (> pivot tol) (< pivot (- tol)))
 723       (dotimes (i rows)
 724         (declare (fixnum i))
 725         (dotimes (j cols)
 726           (declare (fixnum j))
 727           (when (and (/= i k) (/= j k))
 728             (setf aij (makedouble (aref data (+ (* cols i) j))))
 729             (setf aik (makedouble (aref data (+ (* cols i) k))))
 730             (setf akj (makedouble (aref data (+ (* cols k) j))))
 731             (setf aij (- aij (/ (* aik akj) pivot)))
 732             (setf (aref data (+ (* cols i) j)) aij))))
 733
 734       (dotimes (i rows)
 735         (declare (fixnum i))
 736         (setf aik (makedouble (aref data (+ (* cols i) k))))
 737         (when (/= i k)
 738           (setf aik (/ aik pivot))
 739           (setf (aref data (+ (* cols i) k)) aik)))
 740
 741       (dotimes (j cols)
 742         (declare (fixnum j))
 743         (setf akj (makedouble (aref data (+ (* cols k) j))))
 744         (when (/= j k)
 745           (setf akj (- (/ akj pivot)))
 746           (setf (aref data (+ (* cols k) j)) akj)))
 747
 748       (setf akk (/ 1.0 pivot))
 749       (setf (aref data (+ (* cols k) k)) akk)
 750       1)
 751      (t 0))))
 752
 753 (defun make-sweep-matrix (x y &optional w)
 754 "Args: (x y &optional weights)
 755 X is a matrix, Y and WEIGHTS are sequences. Returns the sweep matrix for the
 756 (possibly weighted) regression of Y on X."
 757   (check-matrix x)
 758   (check-sequence y)
 759   (if w (check-sequence w))
 760   (let ((n (num-rows x))
 761         (p (num-cols x)))
 762     (if (/= n (length y)) (error "dimensions do not match"))
 763     (if (and w (/= n (length w))) (error "dimensions do not match"))
 764     (let ((mode (max (la-data-mode x)
 765                      (la-data-mode x)
 766                      (if w (la-data-mode w) 0)))
 767           (result (make-array (list (+ p 2) (+ p 2))))
 768           (x-mean (make-array p))
 769           (y (coerce y 'vector))
 770           (w (if w (coerce w 'vector)))
 771           (has-w (if w 1 0)))
 772       (make-sweep-front x y w n p mode has-w x-mean result)
 773       result)))
 774
 775 (defun sweep-in-place (a k tol)
 776   (check-matrix a)
 777   (check-one-fixnum k)
 778   (check-one-real tol)
 779   (let ((rows (num-rows a))
 780         (cols (num-cols a))
 781         (mode (la-data-mode a)))
 782     (let ((swept (sweep-in-place-front a rows cols mode k tol)))
 783       (if (/= 0 swept) t nil))))
 784
 785 (defun sweep-operator (a columns &optional tolerances)
 786 "Args: (a indices &optional tolerances)
 787 A is a matrix, INDICES a sequence of the column indices to be swept. Returns
 788 a list of the swept result and the list of the columns actually swept. (See
 789 MULTREG documentation.) If supplied, TOLERANCES should be a list of real
 790 numbers the same length as INDICES. An index will only be swept if its pivot
 791 element is larger than the corresponding element of TOLERANCES."
 792   (check-matrix a)
 793   (check-sequence columns)
 794   (if tolerances (check-sequence tolerances))
 795   (check-real a)
 796   (check-fixnum columns)
 797   (if tolerances (check-real tolerances))
 798   (do ((tol .0000001)
 799        (result (copy-array a))
 800        (swept-columns nil)
 801        (columns (coerce columns 'list) (cdr columns))
 802        (tolerances (if (consp tolerances) (coerce tolerances 'list))
 803                    (if (consp tolerances) (cdr tolerances))))
 804       ((null columns) (list result swept-columns))
 805       (let ((col (first columns))
 806             (tol (if (consp tolerances) (first tolerances) tol)))
 807         (if (sweep-in-place result col tol)
 808             (setf swept-columns (cons col swept-columns))))))
 809
 810
 811 ;;;;
 812 ;;;; AX+Y
 813 ;;;;
 814
 815 ;;;*** this could probably be made more efficient ***
 816 (defun ax+y (a x y &optional lower)
 817 "Args (a x y &optional lower)
 818 Returns (+ (matmult A X) Y). If LOWER is not nil, A is taken to be lower
 819 triangular."
 820   (check-square-matrix a)
 821   (check-sequence x)
 822   (check-sequence y)
 823   (check-real a)
 824   (check-real x)
 825   (check-real y)
 826   (let* ((n (num-rows a))
 827          (result (make-list n))
 828          (a (compound-data-seq a)))
 829     (declare (fixnum n))
 830     (if (or (/= n (length x)) (/= n (length y)))
 831         (error "dimensions do not match"))
 832     (do* ((tx (make-next-element x) (make-next-element x))
 833           (ty (make-next-element y))
 834           (tr (make-next-element result))
 835           (i 0 (+ i 1))
 836           (start 0 (+ start n))
 837           (end (if lower (+ i 1) n) (if lower (+ i 1) n)))
 838          ((<= n i) result)
 839        (declare (fixnum i start end))
 840        (let ((val (get-next-element ty i)))
 841          (dotimes (j end)
 842            (declare (fixnum j))
 843            (setf val (+ val (* (get-next-element tx j)
 844                                (aref a (+ start j))))))
 845          (set-next-element tr i val)))))
 846
 847 ;;;;
 848 ;;;; Maximization and Numerical Derivatives
 849 ;;;;
 850
 851 (defvar *maximize-callback-function* nil)
 852 (defvar *maximize-callback-arg* nil)
 853
 854 (defun data2double (n data ptr)
 855   (declare (fixnum n))
 856   (let* ((seq (compound-data-seq data))
 857          (elem (make-next-element seq)))
 858     (if (/= (length seq) n) (error "bad data size"))
 859     (dotimes (i n)
 860       (declare (fixnum i))
 861       (la-put-double ptr i (get-next-element elem i)))))
 862
 863 (defun maximize-callback (n px pfval pgrad phess pderivs)
 864   (la-vector-to-data px n mode-re *maximize-callback-arg*)
 865   (let* ((val (funcall *maximize-callback-function* *maximize-callback-arg*))
 866          (derivs (if (consp val) (- (length val) 1) 0)))
 867     (la-put-integer pderivs 0 derivs)
 868     (la-put-double pfval 0 (if (consp val) (first val) val))
 869     (if (<= 1 derivs) (data2double n (second val) pgrad))
 870     (if (<= 2 derivs) (data2double (* n n) (third val) phess))))
 871
 872 (defun numgrad (f x &optional scale (h -1.0))
 873 "Args: (f x &optional scale derivstep)
 874 Computes the numerical gradient of F at X."
 875   (check-sequence x)
 876   (check-real x)
 877   (when scale
 878         (check-sequence scale)
 879         (check-real scale))
 880   (check-one-real h)
 881   (let* ((n (length x))
 882          (result (make-list n)))
 883     (if (and scale (/= n (length scale)))
 884         (error "scale not the same length as x"))
 885     (let ((*maximize-callback-function* f)
 886           (*maximize-callback-arg* (make-list n)))
 887       (let ((px (la-data-to-vector x mode-re))
 888             (pgrad (la-vector n mode-re))
 889             (pscale (la-data-to-vector
 890                      (if scale scale (make-list n :initial-element 1.0))
 891                      mode-re)))
 892         (unwind-protect
 893             (progn
 894               (numgrad-front n px pgrad h pscale)
 895               (la-vector-to-data pgrad n mode-re result))
 896           (la-free-vector px)
 897           (la-free-vector pgrad)
 898           (la-free-vector pscale))))
 899       result))
 900
 901 (defun numhess (f x &optional scale (h -1.0) all)
 902 "Args: (f x &optional scale derivstep)
 903 Computes the numerical Hessian matrix of F at X."
 904   (check-sequence x)
 905   (check-real x)
 906   (when scale
 907         (check-sequence scale)
 908         (check-real scale))
 909   (check-one-real h)
 910   (let* ((n (length x))
 911          (result (if all
 912                      (list nil (make-list n) (make-array (list n n)))
 913                      (make-array (list n n)))))
 914     (if (and scale (/= n (length scale)))
 915         (error "scale not the same length as x"))
 916     (let ((*maximize-callback-function* f)
 917           (*maximize-callback-arg* (make-list n)))
 918       (let ((hess-data (compound-data-seq (if all (third result) result)))
 919             (px (la-data-to-vector x mode-re))
 920             (pf (la-vector 1 mode-re))
 921             (pgrad (la-vector n mode-re))
 922             (phess (la-vector (* n n) mode-re))
 923             (pscale (la-data-to-vector
 924                      (if scale scale (make-list n :initial-element 1.0))
 925                      mode-re)))
 926         (unwind-protect
 927             (progn
 928               (numhess-front n px pf pgrad phess h pscale)
 929               (when all
 930                     (setf (first result) (la-get-double pf 0))
 931                     (la-vector-to-data pgrad n mode-re (second result)))
 932               (la-vector-to-data phess (* n n) mode-re hess-data))
 933           (la-free-vector pf)
 934           (la-free-vector px)
 935           (la-free-vector pgrad)
 936           (la-free-vector phess)
 937           (la-free-vector pscale))))
 938     result))
 939
 940 (defun init-minfo-ipar-values (n ipars)
 941   (let* ((TRUE 1)
 942          (FALSE 0)
 943          (k 0)
 944          (m 0)
 945          (itnlimit -1)
 946          (backtrack TRUE)
 947          (verbose 0)
 948          (vals_suppl FALSE)
 949          (exptilt TRUE)
 950          (count 0)
 951          (termcode 0))
 952     (setf (aref ipars 0) n)
 953     (setf (aref ipars 1) m)
 954     (setf (aref ipars 2) k)
 955     (setf (aref ipars 3) itnlimit)
 956     (setf (aref ipars 4) backtrack)
 957     (setf (aref ipars 5) verbose)
 958     (setf (aref ipars 6) vals_suppl)
 959     (setf (aref ipars 7) exptilt)
 960     (setf (aref ipars 8) count)
 961     (setf (aref ipars 9) termcode)))
 962
 963 (defun init-minfo-dpar-values (h dpars)
 964   (let ((typf 1.0)
 965         (gradtol -1.0)
 966         (steptol -1.0)
 967         (maxstep -1.0)
 968         (dflt 0.0)
 969         (tilt 0.0)
 970         (newtilt 0.0)
 971         (hessadd 0.0))
 972     (setf (aref dpars 0) typf)
 973     (setf (aref dpars 1) h)
 974     (setf (aref dpars 2) gradtol)
 975     (setf (aref dpars 3) steptol)
 976     (setf (aref dpars 4) maxstep)
 977     (setf (aref dpars 5) dflt)
 978     (setf (aref dpars 6) tilt)
 979     (setf (aref dpars 7) newtilt)
 980     (setf (aref dpars 8) hessadd)))
 981
 982 (defun init-minfo-internals (n h internals)
 983   (let ((ipars (aref internals 8))
 984         (dpars (aref internals 9)))
 985     (init-minfo-ipar-values n ipars)
 986     (init-minfo-dpar-values h dpars)))
 987
 988 (defun new-minfo-internals (f x &key scale ((:derivstep h) -1.0))
 989   (check-sequence x)
 990   (check-real x)
 991   (check-one-real h)
 992   (let ((n (length x)))
 993     (when scale
 994           (check-sequence scale)
 995           (check-real scale)
 996           (if (/= n (length scale)) (error "scale and x not the same length")))
 997     (let ((internals (make-array 12)))
 998       (setf (aref internals 0) f)
 999       (setf (aref internals 3) (if (consp x) (copy-list x) (coerce x 'list)))
1000       (setf (aref internals 4)
1001             (if scale (copy-seq scale) (make-array n :initial-element 1.0)))
1002       (setf (aref internals 5) (make-list (+ 1 n (* n n))))
1003       (setf (aref internals 8) (make-array 10))
1004       (setf (aref internals 9) (make-array 9))
1005       (init-minfo-internals n h internals)
1006       internals)))
1007
1008 (defun minfo-maximize (internals &optional verbose)
1009   (let* ((f (aref internals 0))
1010          (x (aref internals 3))
1011          (fvals (aref internals 5))
1012          (n (length x))
1013          (v (if verbose (if (integerp verbose) verbose 1) -1)))
1014     (setf (aref internals 3) (copy-list x))
1015     (setf (aref internals 5) (copy-list fvals))
1016     (let ((*maximize-callback-function* f)
1017           (*maximize-callback-arg* (make-list n)))
1018       (let* ((x (aref internals 3))
1019              (scale (aref internals 4))
1020              (fvals (aref internals 5))
1021              (ip (aref internals 8))
1022              (dp (aref internals 9))
1023              (px (la-data-to-vector x mode-re))
1024              (pscale (la-data-to-vector scale mode-re))
1025              (pfvals (la-vector (length fvals) mode-re))
1026              (pip (la-data-to-vector ip mode-in))
1027              (pdp (la-data-to-vector dp mode-re)))
1028         (unwind-protect
1029             (progn
1030               (base-minfo-maximize px pfvals pscale pip pdp v))
1031           (la-vector-to-data px n mode-re x)
1032           (la-vector-to-data pfvals (+ 1 n (* n n)) mode-re fvals)
1033           (la-vector-to-data pip (length ip) mode-in ip)
1034           (la-vector-to-data pdp (length dp) mode-re dp))
1035         (get-buf)))))
1036
1037 ;;;;
1038 ;;;; Miscellaneous Routines
1039 ;;;;
1040
1041
1042 (defun split-list (x n)
1043 "Args: (list cols)
1044 Returns a list of COLS lists of equal length of the elements of LIST.
1045 Example: (split-list '(1 2 3 4 5 6) 2) returns ((1 2 3) (4 5 6))"
1046   (check-one-fixnum n)
1047   (if (/= (rem (length x) n) 0) (error "length not divisible by ~a" n))
1048   (flet ((next-split ()
1049            (let ((result nil)
1050                  (end nil))
1051              (dotimes (i n result)
1052                (declare (fixnum i))
1053                (let ((c-elem (list (first x))))
1054                  (cond ((null result)
1055                         (setf result c-elem)
1056                         (setf end result))
1057                        (t
1058                         (setf (rest end) c-elem)
1059                         (setf end (rest end)))))
1060                (setf x (rest x))))))
1061     (let ((result nil)
1062           (end nil)
1063           (k (/ (length x) n)))
1064       (declare (fixnum k))
1065       (dotimes (i k result)
1066         (declare (fixnum i))
1067         (let ((c-sub (list (next-split))))
1068           (cond ((null result)
1069                  (setf result c-sub)
1070                  (setf end result))
1071                 (t
1072                  (setf (rest end) c-sub)
1073                  (setf end (rest end)))))))))
1074