lsbasics.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 ;;;; lsbasics -- Low level Lisp-Stat functions
   7 ;;;;
   8 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
   9 ;;;; unrestricted use.
  10
  11 ;;; Package Setup
  12
  13 (defpackage :lisp-stat-basics
  14     (:use :common-lisp
  15           :lisp-stat-object-system
  16           :lisp-stat-types
  17           :lisp-stat-fastmap
  18           :lisp-stat-float
  19           :lisp-stat-macros
  20           :lisp-stat-compound-data
  21           :lisp-stat-sequence
  22           :lisp-stat-matrix
  23           :lisp-stat-linalg
  24           :lisp-stat-probability)
  25   (:shadowing-import-from :lisp-stat-object-system
  26                           slot-value call-method call-next-method)
  27   (:shadowing-import-from :lisp-stat-sequence
  28                           check-sequence)
  29   (:export
  30    ;; lsbasics.lisp
  31    copy-vector copy-array which repeat
  32    permute-array sum prod count-elements mean if-else sample
  33    select
  34
  35    ;; matrices.lisp
  36    ;;   matrixp num-rows num-cols matmult identity-matrix diagonal row-list
  37    ;;   column-list inner-product outer-product cross-product transpose
  38    ;;   bind-columns bind-rows
  39
  40    ;;; linalg.lisp
  41    ;; chol-decomp lu-decomp lu-solve determinant inverse sv-decomp
  42    ;; qr-decomp rcondest make-rotation
  43    ;; fft make-sweep-matrix sweep-operator ax+y numgrad numhess
  44    ;; split-list eigenp
  45
  46    ;; in linalg.lisp, possibly not supported by matlisp
  47    spline kernel-dens kernel-smooth
  48    ;; lispstat-macros
  49    make-rv-function make-rv-function-1
  50    ;; dists
  51    log-gamma uniform-rand normal-cdf normal-quant normal-dens
  52    normal-rand bivnorm-cdf cauchy-cdf cauchy-quant cauchy-dens
  53    cauchy-rand gamma-cdf gamma-quant gamma-dens gamma-rand
  54    chisq-cdf chisq-quant chisq-dens chisq-rand beta-cdf beta-quant
  55    beta-dens beta-rand t-cdf t-quant t-dens t-rand f-cdf f-quant
  56    f-dens f-rand poisson-cdf poisson-quant poisson-pmf poisson-rand
  57    binomial-cdf binomial-quant binomial-pmf binomial-rand
  58    ;;
  59    ))
  60
  61
  62 (in-package #:lisp-stat-basics)
  63
  64
  65 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  66 ;;;;
  67 ;;;;                           Copying Functions
  68 ;;;;
  69 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  70
  71 ;;;
  72 ;;; COPY-VECTOR function
  73 ;;;
  74
  75 (defun copy-vector (x)
  76 "Args: (x)
  77 Returns a copy of the vector X"
  78   (copy-seq x))
  79
  80 ;;;
  81 ;;; COPY-ARRAY function
  82 ;;;
  83
  84 (defun copy-array (a)
  85 "Args: (a)
  86 Returns a copy of the array A"
  87   (vector-to-array (copy-seq (array-data-vector a))
  88                    (array-dimensions a)))
  89
  90 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  91 ;;;;
  92 ;;;;                         Sequence Functions
  93 ;;;;
  94 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  95
  96
  97 ;;;;
  98 ;;;; WHICH function
  99 ;;;;
 100
 101 (defun which (x)
 102 "Args: (x)
 103 Returns a list of the indices where elements of sequence X are not NIL."
 104   (let ((x (list (compound-data-seq x)))
 105         (result nil)
 106         (tail nil))
 107     (flet ((add-result (x)
 108              (if result (setf (rest tail) (list x)) (setf result (list x)))
 109              (setf tail (if tail (rest tail) result)))
 110            (get-next-element (seq-list i)
 111              (cond ((consp (first seq-list))
 112                     (let ((elem (first (first seq-list))))
 113                       (setf (first seq-list) (rest (first seq-list)))
 114                       elem))
 115                    (t (aref (first seq-list) i)))))
 116           (let ((n (length (first x))))
 117             (dotimes (i n result)
 118                      (if (get-next-element x i) (add-result i)))))))
 119
 120 ;;;;
 121 ;;;; REPEAT function
 122 ;;;;
 123
 124 (defun repeat (a b)
 125 "Args: (vals times)
 126 Repeats VALS. If TIMES is a number and VALS is a non-null, non-array atom,
 127 a list of length TIMES with all elements eq to VALS is returned. If VALS
 128 is a list and TIMES is a number then VALS is appended TIMES times. If
 129 TIMES is a list of numbers then VALS must be a list of equal length and
 130 the simpler version of repeat is mapped down the two lists.
 131 Examples: (repeat 2 5)                 returns (2 2 2 2 2)
 132           (repeat '(1 2) 3)            returns (1 2 1 2 1 2)
 133           (repeat '(4 5 6) '(1 2 3))   returns (4 5 5 6 6 6)
 134           (repeat '((4) (5 6)) '(2 3)) returns (4 4 5 6 5 6 5 6)"
 135   (cond ((compound-data-p b)
 136          (let* ((reps (coerce (compound-data-seq (map-elements #'repeat a b))
 137                               'list))
 138                 (result (first reps))
 139                 (tail (last (first reps))))
 140            (dolist (next (rest reps) result)
 141                    (when next
 142                          (setf (rest tail) next)
 143                          (setf tail (last next))))))
 144         (t (let* ((a (if (compound-data-p a)
 145                          (coerce (compound-data-seq a) 'list)
 146                          (list a)))
 147                   (result nil))
 148              (dotimes (i b result)
 149                       (let ((next (copy-list a)))
 150                         (if result (setf (rest (last next)) result))
 151                         (setf result next)))))))
 152
 153 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 154 ;;;;
 155 ;;;;               Subset Selection and Mutation Functions
 156 ;;;;
 157 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 158
 159 ;;;; is x an ordered sequence of nonnegative positive integers?
 160 (defun ordered-nneg-seq(x)
 161   (if (sequencep x)
 162       (let ((n (length x))
 163             (cx (make-next-element x))
 164             (m 0))
 165         (dotimes (i n t)
 166           (let ((elem (check-nonneg-fixnum (get-next-element cx i))))
 167             (if (> m elem) (return nil) (setf m elem)))))))
 168
 169 ;;;; select or set the subsequence corresponding to the specified indices
 170 (defun sequence-select(x indices &optional (values nil set-values))
 171   (let ((rlen 0)
 172         (dlen 0)
 173         (vlen 0)
 174         (data nil)
 175         (result nil))
 176     (declare (fixnum rlen dlen vlen))
 177
 178     ;; Check the input data
 179     (check-sequence x)
 180     (check-sequence indices)
 181     (if set-values (check-sequence values))
 182
 183     ;; Find the data sizes
 184     (setf data (if (ordered-nneg-seq indices) x (coerce x 'vector)))
 185     (setf dlen (length data))
 186     (setf rlen (length indices))
 187     (when set-values
 188       (setf vlen (length values))
 189       (if (/= vlen rlen) (error "value and index sequences do not match")))
 190
 191     ;; set up the result/value sequence
 192     (setf result
 193           (if set-values
 194               values
 195               (make-sequence (if (listp x) 'list 'vector) rlen)))
 196
 197     ;; get or set the sequence elements
 198     (if set-values
 199       (do ((nextx x)
 200            (cr (make-next-element result))
 201            (ci (make-next-element indices))
 202            (i 0 (+ i 1))
 203            (j 0)
 204            (index 0))
 205           ((>= i rlen))
 206         (declare (fixnum i j index))
 207         (setf index (get-next-element ci i))
 208         (if (<= dlen index) (error "index out of range - ~a" index))
 209         (let ((elem (get-next-element cr i)))
 210           (cond
 211            ((listp x)
 212             (when (> j index)
 213               (setf j 0)
 214               (setf nextx x))
 215             (do ()
 216                 ((not (and (< j index) (consp nextx))))
 217               (incf j 1)
 218               (setf nextx (rest nextx)))
 219             (setf (first nextx) elem))
 220            (t (setf (aref x index) elem)))))
 221       (do ((nextx data)
 222            (cr (make-next-element result))
 223            (ci (make-next-element indices))
 224            (i 0 (+ i 1))
 225            (j 0)
 226            (index 0)
 227            (elem nil))
 228           ((>= i rlen))
 229         (declare (fixnum i j index))
 230         (setf index (get-next-element ci i))
 231         (if (<= dlen index) (error "index out of range - ~a" index))
 232         (cond
 233          ((listp data) ;; indices must be ordered
 234           (do ()
 235               ((not (and (< j index) (consp nextx))))
 236             (incf j 1)
 237             (setf nextx (rest nextx)))
 238           (setf elem (first nextx)))
 239          (t (setf elem (aref data index))))
 240         (set-next-element cr i elem)))
 241
 242     result))
 243
 244 (defun old-rowmajor-index (index indices dim olddim)
 245   "translate row major index in resulting subarray to row major index
 246 in the original array."
 247   (declare (fixnum index))
 248   (let ((rank (length dim))
 249         (face 1)
 250         (oldface 1)
 251         (oldindex 0))
 252     (declare (fixnum rank face oldface))
 253
 254     (dotimes (i rank)
 255       (declare (fixnum i))
 256       (setf face (* face (aref dim i)))
 257       (setf oldface (* oldface (aref olddim i))))
 258
 259     (dotimes (i rank)
 260       (declare (fixnum i))
 261       (setf face (/ face (aref dim i)))
 262       (setf oldface (/ oldface (aref olddim i)))
 263       (incf oldindex
 264             (* oldface (aref (aref indices i) (floor (/ index face))))) ;;*** is this floor really needed???
 265       (setf index (rem index face)))
 266     oldindex))
 267
 268 (defun subarray-select (a indexlist &optional (values nil set_values))
 269   "extract or set subarray for the indices from a displaced array."
 270   (let ((indices nil)
 271         (index)
 272         (dim)
 273         (vdim)
 274         (data)
 275         (result_data)
 276         (olddim)
 277         (result)
 278         (rank 0)
 279         (n 0)
 280         (k 0))
 281     (declare (fixnum rank n))
 282
 283     (if (or (sequencep a) (not (arrayp a))) (error "not an array - ~a" a))
 284     (if (not (listp indexlist))  (error "bad index list - ~a" indices))
 285     (if (/= (length indexlist)  (array-rank a))
 286         (error "wrong number of indices"))
 287
 288     (setf indices (coerce indexlist 'vector))
 289
 290     (setf olddim (coerce (array-dimensions a) 'vector))
 291
 292     ;; compute the result dimension vector and fix up the indices
 293     (setf rank (array-rank a))
 294     (setf dim (make-array rank))
 295     (dotimes (i rank)
 296       (declare (fixnum i))
 297       (setf index (aref indices i))
 298       (setf n (aref olddim i))
 299       (setf index (if (fixnump index) (vector index) (coerce index 'vector)))
 300       (setf k (length index))
 301       (dotimes (j k)
 302         (declare (fixnum j))
 303         (if (<= n (check-nonneg-fixnum (aref index j)))
 304           (error "index out of bounds - ~a" (aref index j)))
 305         (setf (aref indices i) index))
 306       (setf (aref dim i) (length index)))
 307
 308     ;; set up the result or check the values
 309     (let ((dim-list (coerce dim 'list)))
 310       (cond
 311        (set_values
 312         (cond
 313          ((compound-data-p values)
 314           (if (or (not (arrayp values)) (/= rank (array-rank values)))
 315             (error "bad values array - ~a" values))
 316           (setf vdim (coerce (array-dimensions values) 'vector))
 317           (dotimes (i rank)
 318             (declare (fixnum i))
 319             (if (/= (aref vdim i) (aref dim i))
 320               (error "bad value array dimensions - ~a" values)))
 321           (setf result values))
 322          (t (setf result (make-array dim-list :initial-element values)))))
 323        (t (setf result (make-array dim-list)))))
 324
 325     ;; compute the result or set the values
 326     (setf data (compound-data-seq a))
 327     (setf result_data (compound-data-seq result))
 328     (setf n (length result_data))
 329     (dotimes (i n)
 330       (declare (fixnum i))
 331       (setf k (old-rowmajor-index i indices dim olddim))
 332       (if (or (> 0 k) (>= k (length data))) (error "index out of range"))
 333       (if set_values
 334         (setf (aref data k) (aref result_data i))
 335         (setf (aref result_data i) (aref data k))))
 336
 337     result))
 338
 339 ;;;;
 340 ;;;; SELECT function
 341 ;;;;
 342
 343 (defun select (x &rest args)
 344 "Args: (a &rest indices)
 345 A can be a list or an array. If A is a list and INDICES is a single number
 346 then the appropriate element of A is returned. If  is a list and INDICES is
 347 a list of numbers then the sublist of the corresponding elements is returned.
 348 If A in an array then the number of INDICES must match the ARRAY-RANK of A.
 349 If each index is a number then the appropriate array element is returned.
 350 Otherwise the INDICES must all be lists of numbers and the corresponding
 351 submatrix of A is returned. SELECT can be used in setf."
 352   (cond
 353    ((every #'fixnump args)
 354     (if (listp x) (nth (first args) x) (apply #'aref x args)))
 355    ((sequencep x) (sequence-select x (first args)))
 356    (t (subarray-select x args))))
 357
 358
 359 ;; Built in SET-SELECT (SETF method for SELECT)
 360 (defun set-select (x &rest args)
 361   (let ((indices (butlast args))
 362         (values (first (last args))))
 363     (cond
 364      ((sequencep x)
 365       (if (not (consp indices)) (error "bad indices - ~a" indices))
 366       (let* ((indices (first indices))
 367              (i-list (if (fixnump indices) (list indices) indices))
 368              (v-list (if (fixnump indices) (list values) values)))
 369         (sequence-select x i-list v-list)))
 370      ((arrayp x)
 371       (subarray-select x indices values))
 372      (t (error "bad argument type - ~a" x)))
 373     values))
 374
 375 (defsetf select set-select)
 376
 377
 378 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 379 ;;;;
 380 ;;;;                    Array Permutation Functions
 381 ;;;;
 382 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 383
 384 (defun permute-indices (x y perm check)
 385   "Args: (x y perm check).
 386 permute x into y using perm; all should be vectors; If check is TRUE
 387 the routine will check to make sure no indices are reused, but x
 388 will be destroyed."
 389   (let ((rank (length x)))
 390     (declare (fixnum rank))
 391     (dotimes (i rank)
 392       (declare (fixnum i))
 393       (let ((k (aref perm i)))
 394         (if (not (fixnump k)) (error "bad permutation sequence - ~a" perm))
 395         (if (or (< k 0) (>= k rank))
 396             (error "bad permutation sequence - ~a" perm))
 397         (setf (aref y i) (aref x k))
 398         ;; to insure dimensions are not re-used
 399         (if check (setf (aref x k) NIL))))))
 400
 401 (defun indices-from-rowmajor (a k result)
 402   "Args: (a k result).
 403 Compute indices in a from rowmajor index k, put in vector result."
 404   (declare (fixnum k))
 405
 406   (if (not (arrayp a)) (error "not an array - ~a" a))
 407   (if (or (> 0 k) (>= k (array-total-size a))) (error "index out of range"))
 408
 409   (let ((face 1)
 410         (rank (array-rank a))
 411         (dim (array-dimensions a)))
 412     (declare (fixnum face rank))
 413     (let ((cdim (make-next-element dim)))
 414       (dotimes (i rank)
 415         (declare (fixnum i))
 416         (setf face (* face (get-next-element cdim i)))))
 417     (let ((cdim (make-next-element dim)))
 418       (dotimes (i rank)
 419         (setf face (/ face (get-next-element cdim i)))
 420         (setf (aref result i) (floor (/ k face)))
 421         (setf k (rem k face))))))
 422
 423 (defun translate-index (i result x perm indices oldindices ilist)
 424   "Args: (i result x perm indices oldindices ilist).
 425 Translate row major index in original array to row major index in new
 426 array. Use indices vectors and ilist for temporary storage."
 427   (declare (fixnum i))
 428   (let ((rank (array-rank x)))
 429     (declare (fixnum rank))
 430     (indices-from-rowmajor x i oldindices)
 431     (permute-indices oldindices indices perm nil)
 432     (do ((next ilist (rest next))
 433          (k 0 (+ k 1)))
 434         ((not (and (< k rank) (consp next))))
 435       (setf (first next) (aref indices k)))
 436     (apply #'array-row-major-index result ilist)))
 437
 438 (defun permute-array (x perm)
 439   "Args: (a p)
 440 Returns a copy of the array A permuted according to the permutation P."
 441   (if (not (arrayp x)) (error "not an array - ~a" x))
 442   (check-sequence perm)
 443   (if (/= (length perm) (array-rank x))
 444     (error "bad permutation sequence - ~a" perm))
 445   (let* ((perm (coerce perm 'vector))
 446          (rank (array-rank x))
 447          (dim (make-array rank))
 448          (olddim (coerce (array-dimensions x) 'vector)))
 449     (declare (fixnum rank))
 450     ;; construct new dimension vector
 451     (permute-indices olddim dim perm t)
 452     ;; make result array and the index vectors and lists */
 453     (let* ((result (make-array (coerce dim 'list)))
 454           (indices (make-array rank))
 455           (oldindices (make-array rank))
 456           (ilist (make-list rank))
 457           (data (compound-data-seq x))
 458           (result_data (compound-data-seq result))
 459           (n (length data)))
 460       (declare (fixnum n))
 461       (dotimes (i rank)
 462         (declare (fixnum i))
 463         (setf (aref oldindices i) (list nil)))
 464       ;; fill in the result
 465       (if (/= n (length result_data)) (error "bad data"))
 466       (dotimes (i n result)
 467         (declare (fixnum i))
 468         (let ((k (translate-index i result x perm indices oldindices ilist)))
 469           (declare (fixnum k))
 470           (setf (aref result_data k) (aref data i)))))))
 471
 472 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 473 ;;;;
 474 ;;;;              SUM, PROD, COUNT-ELEMENTS, and MEAN Functions
 475 ;;;;
 476 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 477
 478 (defun sum-1 (x)
 479   (if (numberp x)
 480     x
 481     (let ((seq (compound-data-seq x))
 482           (sum 0))
 483       (if (consp seq)
 484         (dolist (x seq sum)
 485           (setf sum (+ sum (if (numberp x) x (sum-1 x)))))
 486         (let ((n (length seq)))
 487           (declare (fixnum n))
 488           (dotimes (i n sum)
 489             (declare (fixnum i))
 490             (let ((x (aref seq i)))
 491               (setf sum (+ sum (if (numberp x) x (sum-1 x)))))))))))
 492
 493 (defun sum (&rest args)
 494 "Args: (&rest number-data)
 495 Returns the sum of all the elements of its arguments. Returns 0 if there
 496 are no arguments. Vector reducing."
 497   (if args
 498     (sum-1 (if (rest args) args (first args)))
 499     0))
 500
 501 (defun prod-1 (x)
 502   (if (numberp x)
 503     x
 504     (let ((seq (compound-data-seq x))
 505           (prod 1))
 506       (if (consp seq)
 507         (dolist (x seq prod)
 508           (setf prod (* prod (if (numberp x) x (prod-1 x)))))
 509         (let ((n (length seq)))
 510           (declare (fixnum n))
 511           (dotimes (i n prod)
 512             (declare (fixnum i))
 513             (let ((x (aref seq i)))
 514               (setf prod (* prod (if (numberp x) x (prod-1 x)))))))))))
 515
 516 (defun prod (&rest args)
 517 "Args: (&rest number-data)
 518 Returns the product of all the elements of its arguments. Returns 1 if there
 519 are no arguments. Vector reducing."
 520   (if args
 521     (prod-1 (if (rest args) args (first args)))
 522     1))
 523
 524 (defun count-elements (x)
 525 "Args: (number &rest more-numbers)
 526 Returns the number of its arguments. Vector reducing"
 527   (if (compound-data-p x)
 528     (let ((seq (compound-data-seq x))
 529           (count 0))
 530       (if (consp seq)
 531         (dolist (x seq count)
 532           (incf count (if (compound-data-p x) (count-elements x) 1)))
 533         (let ((n (length seq)))
 534           (declare (fixnum n))
 535           (dotimes (i n count)
 536             (declare (fixnum i))
 537             (let ((x (aref seq i)))
 538               (incf count (if (compound-data-p x) (count-elements x) 1)))))))
 539     1))
 540
 541
 542 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 543 ;;;;
 544 ;;;;                    IF-ELSE Functions
 545 ;;;;
 546 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 547
 548 (defun if-else (a x y)
 549 "Args: (first x y)
 550 Takes simple or compound data items FIRST, X and Y and returns result of
 551 elementswise selecting from X if FIRST is not NIL and from Y otherwise."
 552   (flet ((base-if-else (a x y) (if a x y)))
 553     (recursive-map-elements #'base-if-else #'if-else a x y)))