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