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    permute-array sum prod count-elements mean if-else sample
  31
  32    ;; matrices.lisp
  33    ;;   matrixp num-rows num-cols matmult identity-matrix diagonal row-list
  34    ;;   column-list inner-product outer-product cross-product transpose
  35    ;;   bind-columns bind-rows
  36
  37    ;;; linalg.lisp
  38    ;; chol-decomp lu-decomp lu-solve determinant inverse sv-decomp
  39    ;; qr-decomp rcondest make-rotation
  40    ;; fft make-sweep-matrix sweep-operator ax+y numgrad numhess
  41    ;; split-list eigenp
  42
  43    ;; in linalg.lisp, possibly not supported by matlisp
  44    spline kernel-dens kernel-smooth
  45
  46    ;; lispstat-macros
  47    make-rv-function make-rv-function-1
  48    ;; dists
  49    log-gamma uniform-rand normal-cdf normal-quant normal-dens
  50    normal-rand bivnorm-cdf cauchy-cdf cauchy-quant cauchy-dens
  51    cauchy-rand gamma-cdf gamma-quant gamma-dens gamma-rand
  52    chisq-cdf chisq-quant chisq-dens chisq-rand beta-cdf beta-quant
  53    beta-dens beta-rand t-cdf t-quant t-dens t-rand f-cdf f-quant
  54    f-dens f-rand poisson-cdf poisson-quant poisson-pmf poisson-rand
  55    binomial-cdf binomial-quant binomial-pmf binomial-rand
  56    ;;
  57    ))
  58
  59
  60 (in-package #:lisp-stat-basics)
  61
  62 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  63 ;;;;
  64 ;;;;               Subset Selection and Mutation Functions
  65 ;;;;
  66 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  67
  68
  69 (defun old-rowmajor-index (index indices dim olddim)
  70   "translate row major index in resulting subarray to row major index
  71 in the original array."
  72   (declare (fixnum index))
  73   (let ((rank (length dim))
  74         (face 1)
  75         (oldface 1)
  76         (oldindex 0))
  77     (declare (fixnum rank face oldface))
  78
  79     (dotimes (i rank)
  80       (declare (fixnum i))
  81       (setf face (* face (aref dim i)))
  82       (setf oldface (* oldface (aref olddim i))))
  83
  84     (dotimes (i rank)
  85       (declare (fixnum i))
  86       (setf face (/ face (aref dim i)))
  87       (setf oldface (/ oldface (aref olddim i)))
  88       (incf oldindex
  89             (* oldface (aref (aref indices i) (floor (/ index face))))) ;;*** is this floor really needed???
  90       (setf index (rem index face)))
  91     oldindex))
  92
  93 (defun subarray-select (a indexlist &optional (values nil set_values))
  94   "extract or set subarray for the indices from a displaced array."
  95   (let ((indices nil)
  96         (index)
  97         (dim)
  98         (vdim)
  99         (data)
 100         (result_data)
 101         (olddim)
 102         (result)
 103         (rank 0)
 104         (n 0)
 105         (k 0))
 106     (declare (fixnum rank n))
 107
 108     (if (or (sequencep a) (not (arrayp a))) (error "not an array - ~a" a))
 109     (if (not (listp indexlist))  (error "bad index list - ~a" indices))
 110     (if (/= (length indexlist)  (array-rank a))
 111         (error "wrong number of indices"))
 112
 113     (setf indices (coerce indexlist 'vector))
 114
 115     (setf olddim (coerce (array-dimensions a) 'vector))
 116
 117     ;; compute the result dimension vector and fix up the indices
 118     (setf rank (array-rank a))
 119     (setf dim (make-array rank))
 120     (dotimes (i rank)
 121       (declare (fixnum i))
 122       (setf index (aref indices i))
 123       (setf n (aref olddim i))
 124       (setf index (if (fixnump index) (vector index) (coerce index 'vector)))
 125       (setf k (length index))
 126       (dotimes (j k)
 127         (declare (fixnum j))
 128         (if (<= n (check-nonneg-fixnum (aref index j)))
 129           (error "index out of bounds - ~a" (aref index j)))
 130         (setf (aref indices i) index))
 131       (setf (aref dim i) (length index)))
 132
 133     ;; set up the result or check the values
 134     (let ((dim-list (coerce dim 'list)))
 135       (cond
 136        (set_values
 137         (cond
 138          ((compound-data-p values)
 139           (if (or (not (arrayp values)) (/= rank (array-rank values)))
 140             (error "bad values array - ~a" values))
 141           (setf vdim (coerce (array-dimensions values) 'vector))
 142           (dotimes (i rank)
 143             (declare (fixnum i))
 144             (if (/= (aref vdim i) (aref dim i))
 145               (error "bad value array dimensions - ~a" values)))
 146           (setf result values))
 147          (t (setf result (make-array dim-list :initial-element values)))))
 148        (t (setf result (make-array dim-list)))))
 149
 150     ;; compute the result or set the values
 151     (setf data (compound-data-seq a))
 152     (setf result_data (compound-data-seq result))
 153     (setf n (length result_data))
 154     (dotimes (i n)
 155       (declare (fixnum i))
 156       (setf k (old-rowmajor-index i indices dim olddim))
 157       (if (or (> 0 k) (>= k (length data))) (error "index out of range"))
 158       (if set_values
 159         (setf (aref data k) (aref result_data i))
 160         (setf (aref result_data i) (aref data k))))
 161
 162     result))
 163
 164
 165 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 166 ;;;;
 167 ;;;;                    Array Permutation Functions
 168 ;;;;
 169 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 170
 171 (defun permute-indices (x y perm check)
 172   "Args: (x y perm check).
 173 permute x into y using perm; all should be vectors; If check is TRUE
 174 the routine will check to make sure no indices are reused, but x
 175 will be destroyed."
 176   (let ((rank (length x)))
 177     (declare (fixnum rank))
 178     (dotimes (i rank)
 179       (declare (fixnum i))
 180       (let ((k (aref perm i)))
 181         (if (not (fixnump k)) (error "bad permutation sequence - ~a" perm))
 182         (if (or (< k 0) (>= k rank))
 183             (error "bad permutation sequence - ~a" perm))
 184         (setf (aref y i) (aref x k))
 185         ;; to insure dimensions are not re-used
 186         (if check (setf (aref x k) NIL))))))
 187
 188 (defun indices-from-rowmajor (a k result)
 189   "Args: (a k result).
 190 Compute indices in a from rowmajor index k, put in vector result."
 191   (declare (fixnum k))
 192
 193   (if (not (arrayp a)) (error "not an array - ~a" a))
 194   (if (or (> 0 k) (>= k (array-total-size a))) (error "index out of range"))
 195
 196   (let ((face 1)
 197         (rank (array-rank a))
 198         (dim (array-dimensions a)))
 199     (declare (fixnum face rank))
 200     (let ((cdim (make-next-element dim)))
 201       (dotimes (i rank)
 202         (declare (fixnum i))
 203         (setf face (* face (get-next-element cdim i)))))
 204     (let ((cdim (make-next-element dim)))
 205       (dotimes (i rank)
 206         (setf face (/ face (get-next-element cdim i)))
 207         (setf (aref result i) (floor (/ k face)))
 208         (setf k (rem k face))))))
 209
 210 (defun translate-index (i result x perm indices oldindices ilist)
 211   "Args: (i result x perm indices oldindices ilist).
 212 Translate row major index in original array to row major index in new
 213 array. Use indices vectors and ilist for temporary storage."
 214   (declare (fixnum i))
 215   (let ((rank (array-rank x)))
 216     (declare (fixnum rank))
 217     (indices-from-rowmajor x i oldindices)
 218     (permute-indices oldindices indices perm nil)
 219     (do ((next ilist (rest next))
 220          (k 0 (+ k 1)))
 221         ((not (and (< k rank) (consp next))))
 222       (setf (first next) (aref indices k)))
 223     (apply #'array-row-major-index result ilist)))
 224
 225 (defun permute-array (x perm)
 226   "Args: (a p)
 227 Returns a copy of the array A permuted according to the permutation P."
 228   (if (not (arrayp x)) (error "not an array - ~a" x))
 229   (check-sequence perm)
 230   (if (/= (length perm) (array-rank x))
 231     (error "bad permutation sequence - ~a" perm))
 232   (let* ((perm (coerce perm 'vector))
 233          (rank (array-rank x))
 234          (dim (make-array rank))
 235          (olddim (coerce (array-dimensions x) 'vector)))
 236     (declare (fixnum rank))
 237     ;; construct new dimension vector
 238     (permute-indices olddim dim perm t)
 239     ;; make result array and the index vectors and lists */
 240     (let* ((result (make-array (coerce dim 'list)))
 241           (indices (make-array rank))
 242           (oldindices (make-array rank))
 243           (ilist (make-list rank))
 244           (data (compound-data-seq x))
 245           (result_data (compound-data-seq result))
 246           (n (length data)))
 247       (declare (fixnum n))
 248       (dotimes (i rank)
 249         (declare (fixnum i))
 250         (setf (aref oldindices i) (list nil)))
 251       ;; fill in the result
 252       (if (/= n (length result_data)) (error "bad data"))
 253       (dotimes (i n result)
 254         (declare (fixnum i))
 255         (let ((k (translate-index i result x perm indices oldindices ilist)))
 256           (declare (fixnum k))
 257           (setf (aref result_data k) (aref data i)))))))
 258
 259 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 260 ;;;;
 261 ;;;;              SUM, PROD, COUNT-ELEMENTS, and MEAN Functions
 262 ;;;;
 263 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 264
 265 (defun sum-1 (x)
 266   (if (numberp x)
 267     x
 268     (let ((seq (compound-data-seq x))
 269           (sum 0))
 270       (if (consp seq)
 271         (dolist (x seq sum)
 272           (setf sum (+ sum (if (numberp x) x (sum-1 x)))))
 273         (let ((n (length seq)))
 274           (declare (fixnum n))
 275           (dotimes (i n sum)
 276             (declare (fixnum i))
 277             (let ((x (aref seq i)))
 278               (setf sum (+ sum (if (numberp x) x (sum-1 x)))))))))))
 279
 280 (defun sum (&rest args)
 281 "Args: (&rest number-data)
 282 Returns the sum of all the elements of its arguments. Returns 0 if there
 283 are no arguments. Vector reducing."
 284   (if args
 285     (sum-1 (if (rest args) args (first args)))
 286     0))
 287
 288 (defun prod-1 (x)
 289   (if (numberp x)
 290     x
 291     (let ((seq (compound-data-seq x))
 292           (prod 1))
 293       (if (consp seq)
 294         (dolist (x seq prod)
 295           (setf prod (* prod (if (numberp x) x (prod-1 x)))))
 296         (let ((n (length seq)))
 297           (declare (fixnum n))
 298           (dotimes (i n prod)
 299             (declare (fixnum i))
 300             (let ((x (aref seq i)))
 301               (setf prod (* prod (if (numberp x) x (prod-1 x)))))))))))
 302
 303 (defun prod (&rest args)
 304 "Args: (&rest number-data)
 305 Returns the product of all the elements of its arguments. Returns 1 if there
 306 are no arguments. Vector reducing."
 307   (if args
 308     (prod-1 (if (rest args) args (first args)))
 309     1))
 310
 311 (defun count-elements (x)
 312 "Args: (number &rest more-numbers)
 313 Returns the number of its arguments. Vector reducing"
 314   (if (compound-data-p x)
 315     (let ((seq (compound-data-seq x))
 316           (count 0))
 317       (if (consp seq)
 318         (dolist (x seq count)
 319           (incf count (if (compound-data-p x) (count-elements x) 1)))
 320         (let ((n (length seq)))
 321           (declare (fixnum n))
 322           (dotimes (i n count)
 323             (declare (fixnum i))
 324             (let ((x (aref seq i)))
 325               (incf count (if (compound-data-p x) (count-elements x) 1)))))))
 326     1))
 327
 328
 329 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 330 ;;;;
 331 ;;;;                    IF-ELSE Functions
 332 ;;;;
 333 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 334
 335 (defun if-else (a x y)
 336 "Args: (first x y)
 337 Takes simple or compound data items FIRST, X and Y and returns result of
 338 elementswise selecting from X if FIRST is not NIL and from Y otherwise."
 339   (flet ((base-if-else (a x y) (if a x y)))
 340     (recursive-map-elements #'base-if-else #'if-else a x y)))