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   (:nicknames #:ls-basics)
  17   (:use ;;:common-lisp
  18         :lisp-stat-object-system
  19         :lisp-stat-fastmap
  20         :lisp-stat-float
  21         :lisp-stat-macros)
  22   ;;(:shadow (package-shadowing-symbols 'lisp-stat-object-system))
  23   ;;(:import-from packagename symbols)
  24   ;;(:shadowing-import-from packagename symbols)
  25   (:export
  26    ;; lsbasics.lisp
  27    sequencep copy-vector copy-array iseq which repeat select
  28    permute-array sum prod count-elements mean if-else sample sort-data
  29    order rank
  30    ;; matrices.lisp
  31    matrixp num-rows num-cols matmult identity-matrix diagonal row-list
  32    column-list inner-product outer-product cross-product transpose
  33    bind-columns bind-rows
  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 eigen
  39    ;; in linalg.lisp, possibly not supported by matlisp
  40    spline kernel-dens kernel-smooth
  41    ;; lispstat-macros
  42    make-rv-function make-rv-function-1
  43    ;; dists
  44    log-gamma uniform-rand normal-cdf normal-quant normal-dens
  45    normal-rand bivnorm-cdf cauchy-cdf cauchy-quant cauchy-dens
  46    cauchy-rand gamma-cdf gamma-quant gamma-dens gamma-rand
  47    chisq-cdf chisq-quant chisq-dens chisq-rand beta-cdf beta-quant
  48    beta-dens beta-rand t-cdf t-quant t-dens t-rand f-cdf f-quant
  49    f-dens f-rand poisson-cdf poisson-quant poisson-pmf poisson-rand
  50    binomial-cdf binomial-quant binomial-pmf binomial-rand
  51    ;;
  52    ))
  53
  54
  55 ;; (defpackage :lisp-stat-basics
  56 ;;   (:nicknames :ls-basics)
  57 ;;   (:use ;; :common-lisp
  58 ;;      :lisp-stat-object-system
  59 ;;         :lisp-stat-fastmap
  60 ;;      :lisp-stat-macros)
  61 ;;   ;;(:shadowing-import-from (package-shadowing-symbols #:lisp-stat-object-system))
  62 ;;   (:export
  63
  64 ;;    ;; lsbasics.lsp
  65 ;;    sequencep copy-vector copy-array iseq which repeat select
  66 ;;    permute-array sum prod count-elements mean if-else
  67 ;;    sample sort-data order rank
  68
  69 ;;    ;; kclpatch.lsp
  70 ;;    ;; #+  kcl (export '(function-lambda-expression realp fixnump))
  71
  72 ;;    ;; compound.lsp
  73
  74 ;;    compound-data-p map-elements compound-data-seq
  75 ;;    compound-data-length element-seq compound-data-proto
  76
  77 ;;    ;; dists.lsp
  78 ;;    log-gamma uniform-rand normal-cdf normal-quant normal-dens
  79 ;;    normal-rand bivnorm-cdf cauchy-cdf cauchy-quant cauchy-dens
  80 ;;    cauchy-rand gamma-cdf gamma-quant gamma-dens gamma-rand
  81 ;;    chisq-cdf chisq-quant chisq-dens chisq-rand beta-cdf beta-quant
  82 ;;    beta-dens beta-rand t-cdf t-quant t-dens t-rand f-cdf f-quant
  83 ;;    f-dens f-rand poisson-cdf poisson-quant poisson-pmf poisson-rand
  84 ;;    binomial-cdf binomial-quant binomial-pmf binomial-rand
  85
  86 ;;    ;; linalg.lsp
  87
  88 ;;    chol-decomp lu-decomp lu-solve determinant inverse sv-decomp
  89 ;;    qr-decomp rcondest make-rotation spline kernel-dens kernel-smooth
  90 ;;    fft make-sweep-matrix sweep-operator ax+y numgrad numhess
  91 ;;    split-list eigen
  92
  93 ;;    ;; matrices.lsp
  94 ;;    matrixp num-rows num-cols matmult identity-matrix diagonal
  95 ;;    row-list column-list inner-product outer-product cross-product
  96 ;;    transpose bind-columns bind-rows
  97
  98 ;;    ;; lsfloat.lsp
  99
 100 ;;    +stat-float-typing+ +stat-cfloat-typing+ +stat-float-template+
 101 ;;    machine-epsilon
 102
 103 ;;    ;; mclglue.lsp
 104 ;;    ;; #+:mcl
 105 ;;    ;; (import '(ccl:def-logical-directory ccl:ff-load ccl:deffcfun ccl:defccallable))
 106
 107 ;;    ))
 108
 109 ;; ;;(in-package #:lisp-stat-basics)
 110
 111
 112
 113 (in-package #:lisp-stat-basics)
 114
 115 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 116 ;;;;
 117 ;;;;                      Type Checking Functions
 118 ;;;;
 119 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 120
 121 (defun fixnump (x)
 122 "Args: (x)
 123 Returns T if X is a fixnum; NIL otherwise."
 124   (declare (inline typep))
 125   (typep x 'fixnum))
 126
 127 (defun check-nonneg-fixnum (x)
 128   (if (and (fixnump x) (<= 0 x)) x (error "not a non-negative fixnum")))
 129
 130 (defun check-one-fixnum (x)
 131   (if (not (fixnump x)) (error "not a fixnum - ~a" x)))
 132
 133 (defun check-one-real (a)
 134   (if (not (or (rationalp a) (floatp a))) (error "not a real number ~s" a)))
 135
 136 (defun check-one-number (a)
 137   (if (not (numberp a)) (error "not a number ~s" a)))
 138
 139 (defun check-sequence (a)
 140   (if (not (or (vectorp a) (consp a))) (error "not a sequence - ~s" a)))
 141
 142 (defun check-matrix (a)
 143   (if (not (and (arrayp a) (= (array-rank a) 2)))
 144       (error "not a matrix - ~s" a)))
 145
 146 (defun check-square-matrix (a)
 147   (check-matrix a)
 148   (let ((m (array-dimension a 0))
 149         (n (array-dimension a 1)))
 150     (if (/= n m) (error "not a square matrix - ~s" a))))
 151
 152
 153 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 154 ;;;;
 155 ;;;;                       Sequence Element Access
 156 ;;;;
 157 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 158
 159 (defun get-next-element (x i)
 160   (let ((seq (first x)))
 161     (if (consp seq)
 162         (let ((elem (first seq)))
 163           (setf (first x) (rest seq))
 164           elem)
 165         (aref seq i))))
 166
 167 (defun set-next-element (x i v)
 168   (let ((seq (first x)))
 169     (cond ((consp seq)
 170            (setf (first seq) v)
 171            (setf (first x) (rest seq)))
 172           (t (setf (aref seq i) v)))))
 173
 174 (defun make-next-element (x) (list x))
 175
 176 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 177 ;;;;
 178 ;;;;          Array to Row-Major Data Vector Conversion Functions
 179 ;;;;
 180 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 181
 182 (defun array-data-vector (a)
 183 "Args: (a)
 184 Displaces array A to a vector"
 185   (make-array (array-total-size a)
 186               :displaced-to a
 187               :element-type (array-element-type a)))
 188
 189 (defun vector-to-array (v dims)
 190 "Args: (v dims)
 191 Displaces vector V to array with dimensions DIMS"
 192   (make-array dims
 193               :displaced-to v
 194               :element-type (array-element-type v)))
 195
 196 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 197 ;;;;
 198 ;;;;                           Copying Functions
 199 ;;;;
 200 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 201
 202 ;;;
 203 ;;; COPY-VECTOR function
 204 ;;;
 205
 206 (defun copy-vector (x)
 207 "Args: (x)
 208 Returns a copy of the vector X"
 209   (copy-seq x))
 210
 211 ;;;
 212 ;;; COPY-ARRAY function
 213 ;;;
 214
 215 (defun copy-array (a)
 216 "Args: (a)
 217 Returns a copy of the array A"
 218   (vector-to-array (copy-seq (array-data-vector a))
 219                    (array-dimensions a)))
 220
 221 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 222 ;;;;
 223 ;;;;                         Sequence Functions
 224 ;;;;
 225 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 226
 227 ;;;
 228 ;;; sequence predicate SEQUENCEP
 229 ;;;
 230
 231 (defun sequencep (x)
 232 "Args: (x)
 233 Returns NIL unless X is a list or vector."
 234   (or (listp x) (vectorp x)))
 235
 236 ;;;
 237 ;;; ISEQ - generate a sequence of consecutive integers from a to b
 238 ;;;
 239
 240 (defun iseq (a &optional b)
 241 "Args: (n &optional m)
 242 With one argumant returns a list of consecutive integers from 0 to N - 1.
 243 With two returns a list of consecutive integers from N to M.
 244 Examples: (iseq 4) returns (0 1 2 3)
 245           (iseq 3 7)  returns (3 4 5 6 7)
 246           (iseq 3 -3) returns (3 2 1 0 -1 -2 -3)"
 247   (if b
 248       (let ((n (+ 1 (abs (- b a))))
 249             (x nil))
 250         (dotimes (i n x)
 251                  (setq x (cons (if (< a b) (- b i) (+ b i)) x))))
 252       (cond
 253        ((= 0 a) nil)
 254        ((< a 0) (iseq (+ a 1) 0))
 255        ((< 0 a) (iseq 0 (- a 1))))))
 256
 257 ;;;;
 258 ;;;; WHICH function
 259 ;;;;
 260
 261 (defun which (x)
 262 "Args: (x)
 263 Returns a list of the indices where elements of sequence X are not NIL."
 264   (let ((x (list (compound-data-seq x)))
 265         (result nil)
 266         (tail nil))
 267     (flet ((add-result (x)
 268              (if result (setf (rest tail) (list x)) (setf result (list x)))
 269              (setf tail (if tail (rest tail) result)))
 270            (get-next-element (seq-list i)
 271              (cond ((consp (first seq-list))
 272                     (let ((elem (first (first seq-list))))
 273                       (setf (first seq-list) (rest (first seq-list)))
 274                       elem))
 275                    (t (aref (first seq-list) i)))))
 276           (let ((n (length (first x))))
 277             (dotimes (i n result)
 278                      (if (get-next-element x i) (add-result i)))))))
 279
 280 ;;;;
 281 ;;;; REPEAT function
 282 ;;;;
 283
 284 (defun repeat (a b)
 285 "Args: (vals times)
 286 Repeats VALS. If TIMES is a number and VALS is a non-null, non-array atom,
 287 a list of length TIMES with all elements eq to VALS is returned. If VALS
 288 is a list and TIMES is a number then VALS is appended TIMES times. If
 289 TIMES is a list of numbers then VALS must be a list of equal length and
 290 the simpler version of repeat is mapped down the two lists.
 291 Examples: (repeat 2 5)                 returns (2 2 2 2 2)
 292           (repeat '(1 2) 3)            returns (1 2 1 2 1 2)
 293           (repeat '(4 5 6) '(1 2 3))   returns (4 5 5 6 6 6)
 294           (repeat '((4) (5 6)) '(2 3)) returns (4 4 5 6 5 6 5 6)"
 295   (cond ((compound-data-p b)
 296          (let* ((reps (coerce (compound-data-seq (map-elements #'repeat a b))
 297                               'list))
 298                 (result (first reps))
 299                 (tail (last (first reps))))
 300            (dolist (next (rest reps) result)
 301                    (when next
 302                          (setf (rest tail) next)
 303                          (setf tail (last next))))))
 304         (t (let* ((a (if (compound-data-p a)
 305                          (coerce (compound-data-seq a) 'list)
 306                          (list a)))
 307                   (result nil))
 308              (dotimes (i b result)
 309                       (let ((next (copy-list a)))
 310                         (if result (setf (rest (last next)) result))
 311                         (setf result next)))))))
 312
 313 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 314 ;;;;
 315 ;;;;               Subset Selection and Mutation Functions
 316 ;;;;
 317 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 318
 319 ;;;; is x an ordered sequence of nonnegative positive integers?
 320 (defun ordered-nneg-seq(x)
 321   (if (sequencep x)
 322       (let ((n (length x))
 323             (cx (make-next-element x))
 324             (m 0))
 325         (dotimes (i n t)
 326           (let ((elem (check-nonneg-fixnum (get-next-element cx i))))
 327             (if (> m elem) (return nil) (setf m elem)))))))
 328
 329 ;;;; select or set the subsequence corresponding to the specified indices
 330 (defun sequence-select(x indices &optional (values nil set-values))
 331   (let ((rlen 0)
 332         (dlen 0)
 333         (vlen 0)
 334         (data nil)
 335         (result nil))
 336     (declare (fixnum rlen dlen vlen))
 337
 338     ;; Check the input data
 339     (check-sequence x)
 340     (check-sequence indices)
 341     (if set-values (check-sequence values))
 342
 343     ;; Find the data sizes
 344     (setf data (if (ordered-nneg-seq indices) x (coerce x 'vector)))
 345     (setf dlen (length data))
 346     (setf rlen (length indices))
 347     (when set-values
 348       (setf vlen (length values))
 349       (if (/= vlen rlen) (error "value and index sequences do not match")))
 350
 351     ;; set up the result/value sequence
 352     (setf result
 353           (if set-values
 354               values
 355               (make-sequence (if (listp x) 'list 'vector) rlen)))
 356
 357     ;; get or set the sequence elements
 358     (if set-values
 359       (do ((nextx x)
 360            (cr (make-next-element result))
 361            (ci (make-next-element indices))
 362            (i 0 (+ i 1))
 363            (j 0)
 364            (index 0))
 365           ((>= i rlen))
 366         (declare (fixnum i j index))
 367         (setf index (get-next-element ci i))
 368         (if (<= dlen index) (error "index out of range - ~a" index))
 369         (let ((elem (get-next-element cr i)))
 370           (cond
 371            ((listp x)
 372             (when (> j index)
 373               (setf j 0)
 374               (setf nextx x))
 375             (do ()
 376                 ((not (and (< j index) (consp nextx))))
 377               (incf j 1)
 378               (setf nextx (rest nextx)))
 379             (setf (first nextx) elem))
 380            (t (setf (aref x index) elem)))))
 381       (do ((nextx data)
 382            (cr (make-next-element result))
 383            (ci (make-next-element indices))
 384            (i 0 (+ i 1))
 385            (j 0)
 386            (index 0)
 387            (elem nil))
 388           ((>= i rlen))
 389         (declare (fixnum i j index))
 390         (setf index (get-next-element ci i))
 391         (if (<= dlen index) (error "index out of range - ~a" index))
 392         (cond
 393          ((listp data) ;; indices must be ordered
 394           (do ()
 395               ((not (and (< j index) (consp nextx))))
 396             (incf j 1)
 397             (setf nextx (rest nextx)))
 398           (setf elem (first nextx)))
 399          (t (setf elem (aref data index))))
 400         (set-next-element cr i elem)))
 401
 402     result))
 403
 404 (defun old-rowmajor-index (index indices dim olddim)
 405   "translate row major index in resulting subarray to row major index
 406 in the original array."
 407   (declare (fixnum index))
 408   (let ((rank (length dim))
 409         (face 1)
 410         (oldface 1)
 411         (oldindex 0))
 412     (declare (fixnum rank face oldface))
 413
 414     (dotimes (i rank)
 415       (declare (fixnum i))
 416       (setf face (* face (aref dim i)))
 417       (setf oldface (* oldface (aref olddim i))))
 418
 419     (dotimes (i rank)
 420       (declare (fixnum i))
 421       (setf face (/ face (aref dim i)))
 422       (setf oldface (/ oldface (aref olddim i)))
 423       (incf oldindex
 424             (* oldface (aref (aref indices i) (floor (/ index face))))) ;;*** is this floor really needed???
 425       (setf index (rem index face)))
 426     oldindex))
 427
 428 (defun subarray-select (a indexlist &optional (values nil set_values))
 429   "extract or set subarray for the indices from a displaced array."
 430   (let ((indices nil)
 431         (index)
 432         (dim)
 433         (vdim)
 434         (data)
 435         (result_data)
 436         (olddim)
 437         (result)
 438         (rank 0)
 439         (n 0)
 440         (k 0))
 441     (declare (fixnum rank n))
 442
 443     (if (or (sequencep a) (not (arrayp a))) (error "not an array - ~a" a))
 444     (if (not (listp indexlist))  (error "bad index list - ~a" indices))
 445     (if (/= (length indexlist)  (array-rank a))
 446         (error "wrong number of indices"))
 447
 448     (setf indices (coerce indexlist 'vector))
 449
 450     (setf olddim (coerce (array-dimensions a) 'vector))
 451
 452     ;; compute the result dimension vector and fix up the indices
 453     (setf rank (array-rank a))
 454     (setf dim (make-array rank))
 455     (dotimes (i rank)
 456       (declare (fixnum i))
 457       (setf index (aref indices i))
 458       (setf n (aref olddim i))
 459       (setf index (if (fixnump index) (vector index) (coerce index 'vector)))
 460       (setf k (length index))
 461       (dotimes (j k)
 462         (declare (fixnum j))
 463         (if (<= n (check-nonneg-fixnum (aref index j)))
 464           (error "index out of bounds - ~a" (aref index j)))
 465         (setf (aref indices i) index))
 466       (setf (aref dim i) (length index)))
 467
 468     ;; set up the result or check the values
 469     (let ((dim-list (coerce dim 'list)))
 470       (cond
 471        (set_values
 472         (cond
 473          ((compound-data-p values)
 474           (if (or (not (arrayp values)) (/= rank (array-rank values)))
 475             (error "bad values array - ~a" values))
 476           (setf vdim (coerce (array-dimensions values) 'vector))
 477           (dotimes (i rank)
 478             (declare (fixnum i))
 479             (if (/= (aref vdim i) (aref dim i))
 480               (error "bad value array dimensions - ~a" values)))
 481           (setf result values))
 482          (t (setf result (make-array dim-list :initial-element values)))))
 483        (t (setf result (make-array dim-list)))))
 484
 485     ;; compute the result or set the values
 486     (setf data (compound-data-seq a))
 487     (setf result_data (compound-data-seq result))
 488     (setf n (length result_data))
 489     (dotimes (i n)
 490       (declare (fixnum i))
 491       (setf k (old-rowmajor-index i indices dim olddim))
 492       (if (or (> 0 k) (>= k (length data))) (error "index out of range"))
 493       (if set_values
 494         (setf (aref data k) (aref result_data i))
 495         (setf (aref result_data i) (aref data k))))
 496
 497     result))
 498
 499 ;;;;
 500 ;;;; SELECT function
 501 ;;;;
 502
 503 (defun select (x &rest args)
 504 "Args: (a &rest indices)
 505 A can be a list or an array. If A is a list and INDICES is a single number
 506 then the appropriate element of A is returned. If  is a list and INDICES is
 507 a list of numbers then the sublist of the corresponding elements is returned.
 508 If A in an array then the number of INDICES must match the ARRAY-RANK of A.
 509 If each index is a number then the appropriate array element is returned.
 510 Otherwise the INDICES must all be lists of numbers and the corresponding
 511 submatrix of A is returned. SELECT can be used in setf."
 512   (cond
 513    ((every #'fixnump args)
 514     (if (listp x) (nth (first args) x) (apply #'aref x args)))
 515    ((sequencep x) (sequence-select x (first args)))
 516    (t (subarray-select x args))))
 517
 518
 519 ;; Built in SET-SELECT (SETF method for SELECT)
 520 (defun set-select (x &rest args)
 521   (let ((indices (butlast args))
 522         (values (first (last args))))
 523     (cond
 524      ((sequencep x)
 525       (if (not (consp indices)) (error "bad indices - ~a" indices))
 526       (let* ((indices (first indices))
 527              (i-list (if (fixnump indices) (list indices) indices))
 528              (v-list (if (fixnump indices) (list values) values)))
 529         (sequence-select x i-list v-list)))
 530      ((arrayp x)
 531       (subarray-select x indices values))
 532      (t (error "bad argument type - ~a" x)))
 533     values))
 534
 535 (defsetf select set-select)
 536
 537
 538 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 539 ;;;;
 540 ;;;;                    Array Permutation Functions
 541 ;;;;
 542 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 543
 544 (defun permute-indices (x y perm check)
 545   "Args: (x y perm check).
 546 permute x into y using perm; all should be vectors; If check is TRUE
 547 the routine will check to make sure no indices are reused, but x
 548 will be destroyed."
 549   (let ((rank (length x)))
 550     (declare (fixnum rank))
 551     (dotimes (i rank)
 552       (declare (fixnum i))
 553       (let ((k (aref perm i)))
 554         (if (not (fixnump k)) (error "bad permutation sequence - ~a" perm))
 555         (if (or (< k 0) (>= k rank))
 556             (error "bad permutation sequence - ~a" perm))
 557         (setf (aref y i) (aref x k))
 558         ;; to insure dimensions are not re-used
 559         (if check (setf (aref x k) NIL))))))
 560
 561 (defun indices-from-rowmajor (a k result)
 562   "Args: (a k result).
 563 Compute indices in a from rowmajor index k, put in vector result."
 564   (declare (fixnum k))
 565
 566   (if (not (arrayp a)) (error "not an array - ~a" a))
 567   (if (or (> 0 k) (>= k (array-total-size a))) (error "index out of range"))
 568
 569   (let ((face 1)
 570         (rank (array-rank a))
 571         (dim (array-dimensions a)))
 572     (declare (fixnum face rank))
 573     (let ((cdim (make-next-element dim)))
 574       (dotimes (i rank)
 575         (declare (fixnum i))
 576         (setf face (* face (get-next-element cdim i)))))
 577     (let ((cdim (make-next-element dim)))
 578       (dotimes (i rank)
 579         (setf face (/ face (get-next-element cdim i)))
 580         (setf (aref result i) (floor (/ k face)))
 581         (setf k (rem k face))))))
 582
 583 (defun translate-index (i result x perm indices oldindices ilist)
 584   "Args: (i result x perm indices oldindices ilist).
 585 Translate row major index in original array to row major index in new
 586 array. Use indices vectors and ilist for temporary storage."
 587   (declare (fixnum i))
 588   (let ((rank (array-rank x)))
 589     (declare (fixnum rank))
 590     (indices-from-rowmajor x i oldindices)
 591     (permute-indices oldindices indices perm nil)
 592     (do ((next ilist (rest next))
 593          (k 0 (+ k 1)))
 594         ((not (and (< k rank) (consp next))))
 595       (setf (first next) (aref indices k)))
 596     (apply #'array-row-major-index result ilist)))
 597
 598 (defun permute-array (x perm)
 599   "Args: (a p)
 600 Returns a copy of the array A permuted according to the permutation P."
 601   (if (not (arrayp x)) (error "not an array - ~a" x))
 602   (check-sequence perm)
 603   (if (/= (length perm) (array-rank x))
 604     (error "bad permutation sequence - ~a" perm))
 605   (let* ((perm (coerce perm 'vector))
 606          (rank (array-rank x))
 607          (dim (make-array rank))
 608          (olddim (coerce (array-dimensions x) 'vector)))
 609     (declare (fixnum rank))
 610     ;; construct new dimension vector
 611     (permute-indices olddim dim perm t)
 612     ;; make result array and the index vectors and lists */
 613     (let* ((result (make-array (coerce dim 'list)))
 614           (indices (make-array rank))
 615           (oldindices (make-array rank))
 616           (ilist (make-list rank))
 617           (data (compound-data-seq x))
 618           (result_data (compound-data-seq result))
 619           (n (length data)))
 620       (declare (fixnum n))
 621       (dotimes (i rank)
 622         (declare (fixnum i))
 623         (setf (aref oldindices i) (list nil)))
 624       ;; fill in the result
 625       (if (/= n (length result_data)) (error "bad data"))
 626       (dotimes (i n result)
 627         (declare (fixnum i))
 628         (let ((k (translate-index i result x perm indices oldindices ilist)))
 629           (declare (fixnum k))
 630           (setf (aref result_data k) (aref data i)))))))
 631
 632 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 633 ;;;;
 634 ;;;;              SUM, PROD, COUNT-ELEMENTS, and MEAN Functions
 635 ;;;;
 636 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 637
 638 (defun sum-1 (x)
 639   (if (numberp x)
 640     x
 641     (let ((seq (compound-data-seq x))
 642           (sum 0))
 643       (if (consp seq)
 644         (dolist (x seq sum)
 645           (setf sum (+ sum (if (numberp x) x (sum-1 x)))))
 646         (let ((n (length seq)))
 647           (declare (fixnum n))
 648           (dotimes (i n sum)
 649             (declare (fixnum i))
 650             (let ((x (aref seq i)))
 651               (setf sum (+ sum (if (numberp x) x (sum-1 x)))))))))))
 652
 653 (defun sum (&rest args)
 654 "Args: (&rest number-data)
 655 Returns the sum of all the elements of its arguments. Returns 0 if there
 656 are no arguments. Vector reducing."
 657   (if args
 658     (sum-1 (if (rest args) args (first args)))
 659     0))
 660
 661 (defun prod-1 (x)
 662   (if (numberp x)
 663     x
 664     (let ((seq (compound-data-seq x))
 665           (prod 1))
 666       (if (consp seq)
 667         (dolist (x seq prod)
 668           (setf prod (* prod (if (numberp x) x (prod-1 x)))))
 669         (let ((n (length seq)))
 670           (declare (fixnum n))
 671           (dotimes (i n prod)
 672             (declare (fixnum i))
 673             (let ((x (aref seq i)))
 674               (setf prod (* prod (if (numberp x) x (prod-1 x)))))))))))
 675
 676 (defun prod (&rest args)
 677 "Args: (&rest number-data)
 678 Returns the product of all the elements of its arguments. Returns 1 if there
 679 are no arguments. Vector reducing."
 680   (if args
 681     (prod-1 (if (rest args) args (first args)))
 682     1))
 683
 684 (defun count-elements (x)
 685 "Args: (number &rest more-numbers)
 686 Returns the number of its arguments. Vector reducing"
 687   (if (compound-data-p x)
 688     (let ((seq (compound-data-seq x))
 689           (count 0))
 690       (if (consp seq)
 691         (dolist (x seq count)
 692           (incf count (if (compound-data-p x) (count-elements x) 1)))
 693         (let ((n (length seq)))
 694           (declare (fixnum n))
 695           (dotimes (i n count)
 696             (declare (fixnum i))
 697             (let ((x (aref seq i)))
 698               (incf count (if (compound-data-p x) (count-elements x) 1)))))))
 699     1))
 700
 701 (defun mean (x)
 702 "Args: (x)
 703 Returns the mean of the elements x. Vector reducing."
 704   (let ((mean 0.0)
 705         (count 0.0))
 706     (labels ((add-to-mean (x)
 707               (let ((count+1 (+ count 1.0)))
 708                 (setf mean (+ (* (/ count count+1) mean) (* (/ count+1) x)))
 709                 (setf count count+1)))
 710              (find-mean (x)
 711                (if (numberp x)
 712                  (add-to-mean x)
 713                  (let ((seq (compound-data-seq x)))
 714                    (if (consp seq)
 715                      (dolist (x seq)
 716                        (if (numberp x) (add-to-mean x) (find-mean x)))
 717                      (let ((n (length seq)))
 718                        (dotimes (i n)
 719                          (declare (fixnum i))
 720                          (let ((x (aref seq i)))
 721                            (if (numberp x)
 722                                (add-to-mean x)
 723                                (find-mean x))))))))))
 724       (find-mean x)
 725       mean)))
 726
 727 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 728 ;;;;
 729 ;;;;                         Sorting Functions
 730 ;;;;
 731 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 732
 733 (defun sort-data (x)
 734 "Args: (sequence)
 735 Returns a sequence with the numbers or strings in the sequence X in order."
 736   (flet ((less (x y) (if (numberp x) (< x y) (string-lessp x y))))
 737     (stable-sort (copy-seq (compound-data-seq x)) #'less)))
 738
 739 (defun order (x)
 740 "Args (x)
 741 Returns a sequence of the indices of elements in the sequence of numbers
 742 or strings X in order."
 743   (let* ((seq (compound-data-seq x))
 744          (type (if (consp seq) 'list 'vector))
 745          (i -1))
 746     (flet ((entry (x) (setf i (+ i 1)) (list x i))
 747            (less (a b)
 748                  (let ((x (first a))
 749                        (y (first b)))
 750                    (if (numberp x) (< x y) (string-lessp x y)))))
 751       (let ((sorted-seq (stable-sort (map type #'entry seq) #'less)))
 752         (map type #'second sorted-seq)))))
 753
 754 ;; this isn't destructive -- do we document destructive only, or any
 755 ;; variant?
 756 (defun rank (x)
 757 "Args (x)
 758 Returns a sequence with the elements of the list or array of numbers or
 759 strings X replaced by their ranks."
 760   (let ((ranked-seq (order (order x))))
 761     (make-compound-data (compound-data-shape x) ranked-seq)))
 762
 763 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 764 ;;;;
 765 ;;;;                    IF-ELSE and SAMPLE Functions
 766 ;;;;
 767 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 768
 769 (defun if-else (a x y)
 770 "Args: (first x y)
 771 Takes simple or compound data items FIRST, X and Y and returns result of
 772 elementswise selecting from X if FIRST is not NIL and from Y otherwise."
 773   (flet ((base-if-else (a x y) (if a x y)))
 774     (recursive-map-elements #'base-if-else #'if-else a x y)))
 775
 776 (defun sample (x ssize &optional replace)
 777 "Args: (x n &optional (replace nil))
 778 Returns a list of a random sample of size N from sequence X drawn with or
 779 without replacement."
 780   (check-sequence x)
 781   (let ((n (length x))
 782         (x (if (consp x) (coerce x 'vector) (copy-vector x)))
 783         (result nil))
 784     (if (< 0 n)
 785         (dotimes (i ssize result)
 786                  (let ((j (if replace (random n) (+ i (random (- n i))))))
 787                    (setf result (cons (aref x j) result))
 788                    (unless replace     ;; swap elements i and j
 789                            (let ((temp (aref x i)))
 790                              (setf (aref x i) (aref x j))
 791                              (setf (aref x j) temp))))))))
 792