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