linalg.lsp

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