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