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-matrix
  22           ;;:lisp-stat-linalg
  23           :lisp-stat-probability)
  24   (:shadowing-import-from :lisp-stat-object-system
  25                           slot-value call-method call-next-method)
  26   (:export
  27    permute-array sum prod count-elements mean if-else sample
  28
  29    ;; matrices.lisp
  30    ;;   matrixp num-rows num-cols matmult identity-matrix diagonal row-list
  31    ;;   column-list inner-product outer-product cross-product transpose
  32    ;;   bind-columns bind-rows
  33
  34    ;;; linalg.lisp
  35    ;; chol-decomp lu-decomp lu-solve determinant inverse sv-decomp
  36    ;; qr-decomp rcondest make-rotation
  37    ;; fft make-sweep-matrix sweep-operator ax+y numgrad numhess
  38    ;; split-list eigenp
  39
  40    ;; in linalg.lisp, possibly not supported by matlisp
  41    spline kernel-dens kernel-smooth
  42
  43    ;; lispstat-macros
  44    make-rv-function make-rv-function-1
  45    ;; dists
  46    log-gamma uniform-rand normal-cdf normal-quant normal-dens
  47    normal-rand bivnorm-cdf cauchy-cdf cauchy-quant cauchy-dens
  48    cauchy-rand gamma-cdf gamma-quant gamma-dens gamma-rand
  49    chisq-cdf chisq-quant chisq-dens chisq-rand beta-cdf beta-quant
  50    beta-dens beta-rand t-cdf t-quant t-dens t-rand f-cdf f-quant
  51    f-dens f-rand poisson-cdf poisson-quant poisson-pmf poisson-rand
  52    binomial-cdf binomial-quant binomial-pmf binomial-rand
  53    ;;
  54    ))
  55
  56
  57 (in-package #:lisp-stat-basics)
  58
  59
  60
  61
  62 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  63 ;;;;
  64 ;;;;                    Array Permutation Functions
  65 ;;;;
  66 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  67
  68 (defun permute-indices (x y perm check)
  69   "Args: (x y perm check).
  70 permute x into y using perm; all should be vectors; If check is TRUE
  71 the routine will check to make sure no indices are reused, but x
  72 will be destroyed."
  73   (let ((rank (length x)))
  74     (declare (fixnum rank))
  75     (dotimes (i rank)
  76       (declare (fixnum i))
  77       (let ((k (aref perm i)))
  78         (if (not (fixnump k)) (error "bad permutation sequence - ~a" perm))
  79         (if (or (< k 0) (>= k rank))
  80             (error "bad permutation sequence - ~a" perm))
  81         (setf (aref y i) (aref x k))
  82         ;; to insure dimensions are not re-used
  83         (if check (setf (aref x k) NIL))))))
  84
  85 (defun indices-from-rowmajor (a k result)
  86   "Args: (a k result).
  87 Compute indices in a from rowmajor index k, put in vector result."
  88   (declare (fixnum k))
  89
  90   (if (not (arrayp a)) (error "not an array - ~a" a))
  91   (if (or (> 0 k) (>= k (array-total-size a))) (error "index out of range"))
  92
  93   (let ((face 1)
  94         (rank (array-rank a))
  95         (dim (array-dimensions a)))
  96     (declare (fixnum face rank))
  97     (let ((cdim (make-next-element dim)))
  98       (dotimes (i rank)
  99         (declare (fixnum i))
 100         (setf face (* face (get-next-element cdim i)))))
 101     (let ((cdim (make-next-element dim)))
 102       (dotimes (i rank)
 103         (setf face (/ face (get-next-element cdim i)))
 104         (setf (aref result i) (floor (/ k face)))
 105         (setf k (rem k face))))))
 106
 107 (defun translate-index (i result x perm indices oldindices ilist)
 108   "Args: (i result x perm indices oldindices ilist).
 109 Translate row major index in original array to row major index in new
 110 array. Use indices vectors and ilist for temporary storage."
 111   (declare (fixnum i))
 112   (let ((rank (array-rank x)))
 113     (declare (fixnum rank))
 114     (indices-from-rowmajor x i oldindices)
 115     (permute-indices oldindices indices perm nil)
 116     (do ((next ilist (rest next))
 117          (k 0 (+ k 1)))
 118         ((not (and (< k rank) (consp next))))
 119       (setf (first next) (aref indices k)))
 120     (apply #'array-row-major-index result ilist)))
 121
 122 (defun permute-array (x perm)
 123   "Args: (a p)
 124 Returns a copy of the array A permuted according to the permutation P."
 125   (if (not (arrayp x)) (error "not an array - ~a" x))
 126   (check-sequence perm)
 127   (if (/= (length perm) (array-rank x))
 128     (error "bad permutation sequence - ~a" perm))
 129   (let* ((perm (coerce perm 'vector))
 130          (rank (array-rank x))
 131          (dim (make-array rank))
 132          (olddim (coerce (array-dimensions x) 'vector)))
 133     (declare (fixnum rank))
 134     ;; construct new dimension vector
 135     (permute-indices olddim dim perm t)
 136     ;; make result array and the index vectors and lists */
 137     (let* ((result (make-array (coerce dim 'list)))
 138           (indices (make-array rank))
 139           (oldindices (make-array rank))
 140           (ilist (make-list rank))
 141           (data (compound-data-seq x))
 142           (result_data (compound-data-seq result))
 143           (n (length data)))
 144       (declare (fixnum n))
 145       (dotimes (i rank)
 146         (declare (fixnum i))
 147         (setf (aref oldindices i) (list nil)))
 148       ;; fill in the result
 149       (if (/= n (length result_data)) (error "bad data"))
 150       (dotimes (i n result)
 151         (declare (fixnum i))
 152         (let ((k (translate-index i result x perm indices oldindices ilist)))
 153           (declare (fixnum k))
 154           (setf (aref result_data k) (aref data i)))))))
 155
 156 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 157 ;;;;
 158 ;;;;              SUM, PROD, COUNT-ELEMENTS, and MEAN Functions
 159 ;;;;
 160 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 161
 162 (defun sum-1 (x)
 163   (if (numberp x)
 164     x
 165     (let ((seq (compound-data-seq x))
 166           (sum 0))
 167       (if (consp seq)
 168         (dolist (x seq sum)
 169           (setf sum (+ sum (if (numberp x) x (sum-1 x)))))
 170         (let ((n (length seq)))
 171           (declare (fixnum n))
 172           (dotimes (i n sum)
 173             (declare (fixnum i))
 174             (let ((x (aref seq i)))
 175               (setf sum (+ sum (if (numberp x) x (sum-1 x)))))))))))
 176
 177 (defun sum (&rest args)
 178 "Args: (&rest number-data)
 179 Returns the sum of all the elements of its arguments. Returns 0 if there
 180 are no arguments. Vector reducing."
 181   (if args
 182     (sum-1 (if (rest args) args (first args)))
 183     0))
 184
 185 (defun prod-1 (x)
 186   (if (numberp x)
 187     x
 188     (let ((seq (compound-data-seq x))
 189           (prod 1))
 190       (if (consp seq)
 191         (dolist (x seq prod)
 192           (setf prod (* prod (if (numberp x) x (prod-1 x)))))
 193         (let ((n (length seq)))
 194           (declare (fixnum n))
 195           (dotimes (i n prod)
 196             (declare (fixnum i))
 197             (let ((x (aref seq i)))
 198               (setf prod (* prod (if (numberp x) x (prod-1 x)))))))))))
 199
 200 (defun prod (&rest args)
 201 "Args: (&rest number-data)
 202 Returns the product of all the elements of its arguments. Returns 1 if there
 203 are no arguments. Vector reducing."
 204   (if args
 205     (prod-1 (if (rest args) args (first args)))
 206     1))
 207
 208 (defun count-elements (x)
 209 "Args: (number &rest more-numbers)
 210 Returns the number of its arguments. Vector reducing"
 211   (if (compound-data-p x)
 212     (let ((seq (compound-data-seq x))
 213           (count 0))
 214       (if (consp seq)
 215         (dolist (x seq count)
 216           (incf count (if (compound-data-p x) (count-elements x) 1)))
 217         (let ((n (length seq)))
 218           (declare (fixnum n))
 219           (dotimes (i n count)
 220             (declare (fixnum i))
 221             (let ((x (aref seq i)))
 222               (incf count (if (compound-data-p x) (count-elements x) 1)))))))
 223     1))
 224
 225
 226 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 227 ;;;;
 228 ;;;;                    IF-ELSE Functions
 229 ;;;;
 230 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 231
 232 (defun if-else (a x y)
 233 "Args: (first x y)
 234 Takes simple or compound data items FIRST, X and Y and returns result of
 235 elementswise selecting from X if FIRST is not NIL and from Y otherwise."
 236   (flet ((base-if-else (a x y) (if a x y)))
 237     (recursive-map-elements #'base-if-else #'if-else a x y)))