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