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 ;;; sequence predicate SEQUENCEP
 218 ;;;
 219
 220 (defun sequencep (x)
 221 "Args: (x)
 222 Returns NIL unless X is a list or vector."
 223   (or (listp x) (vectorp x)))
 224
 225 ;;;
 226 ;;; ISEQ - generate a sequence of consecutive integers from a to b
 227 ;;;
 228
 229 (defun iseq (a &optional b)
 230 "Args: (n &optional m)
 231 With one argumant returns a list of consecutive integers from 0 to N - 1.
 232 With two returns a list of consecutive integers from N to M.
 233 Examples: (iseq 4) returns (0 1 2 3)
 234           (iseq 3 7)  returns (3 4 5 6 7)
 235           (iseq 3 -3) returns (3 2 1 0 -1 -2 -3)"
 236   (if b
 237       (let ((n (+ 1 (abs (- b a))))
 238             (x nil))
 239         (dotimes (i n x)
 240                  (setq x (cons (if (< a b) (- b i) (+ b i)) x))))
 241       (cond
 242        ((= 0 a) nil)
 243        ((< a 0) (iseq (+ a 1) 0))
 244        ((< 0 a) (iseq 0 (- a 1))))))
 245
 246 ;;;;
 247 ;;;; WHICH function
 248 ;;;;
 249
 250 (defun which (x)
 251 "Args: (x)
 252 Returns a list of the indices where elements of sequence X are not NIL."
 253   (let ((x (list (compound-data-seq x)))
 254         (result nil)
 255         (tail nil))
 256     (flet ((add-result (x)
 257              (if result (setf (rest tail) (list x)) (setf result (list x)))
 258              (setf tail (if tail (rest tail) result)))
 259            (get-next-element (seq-list i)
 260              (cond ((consp (first seq-list))
 261                     (let ((elem (first (first seq-list))))
 262                       (setf (first seq-list) (rest (first seq-list)))
 263                       elem))
 264                    (t (aref (first seq-list) i)))))
 265           (let ((n (length (first x))))
 266             (dotimes (i n result)
 267                      (if (get-next-element x i) (add-result i)))))))
 268
 269 ;;;;
 270 ;;;; REPEAT function
 271 ;;;;
 272
 273 (defun repeat (a b)
 274 "Args: (vals times)
 275 Repeats VALS. If TIMES is a number and VALS is a non-null, non-array atom,
 276 a list of length TIMES with all elements eq to VALS is returned. If VALS
 277 is a list and TIMES is a number then VALS is appended TIMES times. If
 278 TIMES is a list of numbers then VALS must be a list of equal length and
 279 the simpler version of repeat is mapped down the two lists.
 280 Examples: (repeat 2 5)                 returns (2 2 2 2 2)
 281           (repeat '(1 2) 3)            returns (1 2 1 2 1 2)
 282           (repeat '(4 5 6) '(1 2 3))   returns (4 5 5 6 6 6)
 283           (repeat '((4) (5 6)) '(2 3)) returns (4 4 5 6 5 6 5 6)"
 284   (cond ((compound-data-p b)
 285          (let* ((reps (coerce (compound-data-seq (map-elements #'repeat a b))
 286                               'list))
 287                 (result (first reps))
 288                 (tail (last (first reps))))
 289            (dolist (next (rest reps) result)
 290                    (when next
 291                          (setf (rest tail) next)
 292                          (setf tail (last next))))))
 293         (t (let* ((a (if (compound-data-p a)
 294                          (coerce (compound-data-seq a) 'list)
 295                          (list a)))
 296                   (result nil))
 297              (dotimes (i b result)
 298                       (let ((next (copy-list a)))
 299                         (if result (setf (rest (last next)) result))
 300                         (setf result next)))))))
 301
 302 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 303 ;;;;
 304 ;;;;               Subset Selection and Mutation Functions
 305 ;;;;
 306 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 307
 308 ;;;; is x an ordered sequence of nonnegative positive integers?
 309 (defun ordered-nneg-seq(x)
 310   (if (sequencep x)
 311       (let ((n (length x))
 312             (cx (make-next-element x))
 313             (m 0))
 314         (dotimes (i n t)
 315           (let ((elem (check-nonneg-fixnum (get-next-element cx i))))
 316             (if (> m elem) (return nil) (setf m elem)))))))
 317
 318 ;;;; select or set the subsequence corresponding to the specified indices
 319 (defun sequence-select(x indices &optional (values nil set-values))
 320   (let ((rlen 0)
 321         (dlen 0)
 322         (vlen 0)
 323         (data nil)
 324         (result nil))
 325     (declare (fixnum rlen dlen vlen))
 326
 327     ;; Check the input data
 328     (check-sequence x)
 329     (check-sequence indices)
 330     (if set-values (check-sequence values))
 331
 332     ;; Find the data sizes
 333     (setf data (if (ordered-nneg-seq indices) x (coerce x 'vector)))
 334     (setf dlen (length data))
 335     (setf rlen (length indices))
 336     (when set-values
 337       (setf vlen (length values))
 338       (if (/= vlen rlen) (error "value and index sequences do not match")))
 339
 340     ;; set up the result/value sequence
 341     (setf result
 342           (if set-values
 343               values
 344               (make-sequence (if (listp x) 'list 'vector) rlen)))
 345
 346     ;; get or set the sequence elements
 347     (if set-values
 348       (do ((nextx x)
 349            (cr (make-next-element result))
 350            (ci (make-next-element indices))
 351            (i 0 (+ i 1))
 352            (j 0)
 353            (index 0))
 354           ((>= i rlen))
 355         (declare (fixnum i j index))
 356         (setf index (get-next-element ci i))
 357         (if (<= dlen index) (error "index out of range - ~a" index))
 358         (let ((elem (get-next-element cr i)))
 359           (cond
 360            ((listp x)
 361             (when (> j index)
 362               (setf j 0)
 363               (setf nextx x))
 364             (do ()
 365                 ((not (and (< j index) (consp nextx))))
 366               (incf j 1)
 367               (setf nextx (rest nextx)))
 368             (setf (first nextx) elem))
 369            (t (setf (aref x index) elem)))))
 370       (do ((nextx data)
 371            (cr (make-next-element result))
 372            (ci (make-next-element indices))
 373            (i 0 (+ i 1))
 374            (j 0)
 375            (index 0)
 376            (elem nil))
 377           ((>= i rlen))
 378         (declare (fixnum i j index))
 379         (setf index (get-next-element ci i))
 380         (if (<= dlen index) (error "index out of range - ~a" index))
 381         (cond
 382          ((listp data) ;; indices must be ordered
 383           (do ()
 384               ((not (and (< j index) (consp nextx))))
 385             (incf j 1)
 386             (setf nextx (rest nextx)))
 387           (setf elem (first nextx)))
 388          (t (setf elem (aref data index))))
 389         (set-next-element cr i elem)))
 390
 391     result))
 392
 393 (defun old-rowmajor-index (index indices dim olddim)
 394   "translate row major index in resulting subarray to row major index
 395 in the original array."
 396   (declare (fixnum index))
 397   (let ((rank (length dim))
 398         (face 1)
 399         (oldface 1)
 400         (oldindex 0))
 401     (declare (fixnum rank face oldface))
 402
 403     (dotimes (i rank)
 404       (declare (fixnum i))
 405       (setf face (* face (aref dim i)))
 406       (setf oldface (* oldface (aref olddim i))))
 407
 408     (dotimes (i rank)
 409       (declare (fixnum i))
 410       (setf face (/ face (aref dim i)))
 411       (setf oldface (/ oldface (aref olddim i)))
 412       (incf oldindex
 413             (* oldface (aref (aref indices i) (floor (/ index face))))) ;;*** is this floor really needed???
 414       (setf index (rem index face)))
 415     oldindex))
 416
 417 (defun subarray-select (a indexlist &optional (values nil set_values))
 418   "extract or set subarray for the indices from a displaced array."
 419   (let ((indices nil)
 420         (index)
 421         (dim)
 422         (vdim)
 423         (data)
 424         (result_data)
 425         (olddim)
 426         (result)
 427         (rank 0)
 428         (n 0)
 429         (k 0))
 430     (declare (fixnum rank n))
 431
 432     (if (or (sequencep a) (not (arrayp a))) (error "not an array - ~a" a))
 433     (if (not (listp indexlist))  (error "bad index list - ~a" indices))
 434     (if (/= (length indexlist)  (array-rank a))
 435         (error "wrong number of indices"))
 436
 437     (setf indices (coerce indexlist 'vector))
 438
 439     (setf olddim (coerce (array-dimensions a) 'vector))
 440
 441     ;; compute the result dimension vector and fix up the indices
 442     (setf rank (array-rank a))
 443     (setf dim (make-array rank))
 444     (dotimes (i rank)
 445       (declare (fixnum i))
 446       (setf index (aref indices i))
 447       (setf n (aref olddim i))
 448       (setf index (if (fixnump index) (vector index) (coerce index 'vector)))
 449       (setf k (length index))
 450       (dotimes (j k)
 451         (declare (fixnum j))
 452         (if (<= n (check-nonneg-fixnum (aref index j)))
 453           (error "index out of bounds - ~a" (aref index j)))
 454         (setf (aref indices i) index))
 455       (setf (aref dim i) (length index)))
 456
 457     ;; set up the result or check the values
 458     (let ((dim-list (coerce dim 'list)))
 459       (cond
 460        (set_values
 461         (cond
 462          ((compound-data-p values)
 463           (if (or (not (arrayp values)) (/= rank (array-rank values)))
 464             (error "bad values array - ~a" values))
 465           (setf vdim (coerce (array-dimensions values) 'vector))
 466           (dotimes (i rank)
 467             (declare (fixnum i))
 468             (if (/= (aref vdim i) (aref dim i))
 469               (error "bad value array dimensions - ~a" values)))
 470           (setf result values))
 471          (t (setf result (make-array dim-list :initial-element values)))))
 472        (t (setf result (make-array dim-list)))))
 473
 474     ;; compute the result or set the values
 475     (setf data (compound-data-seq a))
 476     (setf result_data (compound-data-seq result))
 477     (setf n (length result_data))
 478     (dotimes (i n)
 479       (declare (fixnum i))
 480       (setf k (old-rowmajor-index i indices dim olddim))
 481       (if (or (> 0 k) (>= k (length data))) (error "index out of range"))
 482       (if set_values
 483         (setf (aref data k) (aref result_data i))
 484         (setf (aref result_data i) (aref data k))))
 485
 486     result))
 487
 488 ;;;;
 489 ;;;; SELECT function
 490 ;;;;
 491
 492 (defun select (x &rest args)
 493 "Args: (a &rest indices)
 494 A can be a list or an array. If A is a list and INDICES is a single number
 495 then the appropriate element of A is returned. If  is a list and INDICES is
 496 a list of numbers then the sublist of the corresponding elements is returned.
 497 If A in an array then the number of INDICES must match the ARRAY-RANK of A.
 498 If each index is a number then the appropriate array element is returned.
 499 Otherwise the INDICES must all be lists of numbers and the corresponding
 500 submatrix of A is returned. SELECT can be used in setf."
 501   (cond
 502    ((every #'fixnump args)
 503     (if (listp x) (nth (first args) x) (apply #'aref x args)))
 504    ((sequencep x) (sequence-select x (first args)))
 505    (t (subarray-select x args))))
 506
 507
 508 ;; Built in SET-SELECT (SETF method for SELECT)
 509 (defun set-select (x &rest args)
 510   (let ((indices (butlast args))
 511         (values (first (last args))))
 512     (cond
 513      ((sequencep x)
 514       (if (not (consp indices)) (error "bad indices - ~a" indices))
 515       (let* ((indices (first indices))
 516              (i-list (if (fixnump indices) (list indices) indices))
 517              (v-list (if (fixnump indices) (list values) values)))
 518         (sequence-select x i-list v-list)))
 519      ((arrayp x)
 520       (subarray-select x indices values))
 521      (t (error "bad argument type - ~a" x)))
 522     values))
 523
 524 (defsetf select set-select)
 525
 526
 527 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 528 ;;;;
 529 ;;;;                    Array Permutation Functions
 530 ;;;;
 531 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 532
 533 (defun permute-indices (x y perm check)
 534   "Args: (x y perm check).
 535 permute x into y using perm; all should be vectors; If check is TRUE
 536 the routine will check to make sure no indices are reused, but x
 537 will be destroyed."
 538   (let ((rank (length x)))
 539     (declare (fixnum rank))
 540     (dotimes (i rank)
 541       (declare (fixnum i))
 542       (let ((k (aref perm i)))
 543         (if (not (fixnump k)) (error "bad permutation sequence - ~a" perm))
 544         (if (or (< k 0) (>= k rank))
 545             (error "bad permutation sequence - ~a" perm))
 546         (setf (aref y i) (aref x k))
 547         ;; to insure dimensions are not re-used
 548         (if check (setf (aref x k) NIL))))))
 549
 550 (defun indices-from-rowmajor (a k result)
 551   "Args: (a k result).
 552 Compute indices in a from rowmajor index k, put in vector result."
 553   (declare (fixnum k))
 554
 555   (if (not (arrayp a)) (error "not an array - ~a" a))
 556   (if (or (> 0 k) (>= k (array-total-size a))) (error "index out of range"))
 557
 558   (let ((face 1)
 559         (rank (array-rank a))
 560         (dim (array-dimensions a)))
 561     (declare (fixnum face rank))
 562     (let ((cdim (make-next-element dim)))
 563       (dotimes (i rank)
 564         (declare (fixnum i))
 565         (setf face (* face (get-next-element cdim i)))))
 566     (let ((cdim (make-next-element dim)))
 567       (dotimes (i rank)
 568         (setf face (/ face (get-next-element cdim i)))
 569         (setf (aref result i) (floor (/ k face)))
 570         (setf k (rem k face))))))
 571
 572 (defun translate-index (i result x perm indices oldindices ilist)
 573   "Args: (i result x perm indices oldindices ilist).
 574 Translate row major index in original array to row major index in new
 575 array. Use indices vectors and ilist for temporary storage."
 576   (declare (fixnum i))
 577   (let ((rank (array-rank x)))
 578     (declare (fixnum rank))
 579     (indices-from-rowmajor x i oldindices)
 580     (permute-indices oldindices indices perm nil)
 581     (do ((next ilist (rest next))
 582          (k 0 (+ k 1)))
 583         ((not (and (< k rank) (consp next))))
 584       (setf (first next) (aref indices k)))
 585     (apply #'array-row-major-index result ilist)))
 586
 587 (defun permute-array (x perm)
 588   "Args: (a p)
 589 Returns a copy of the array A permuted according to the permutation P."
 590   (if (not (arrayp x)) (error "not an array - ~a" x))
 591   (check-sequence perm)
 592   (if (/= (length perm) (array-rank x))
 593     (error "bad permutation sequence - ~a" perm))
 594   (let* ((perm (coerce perm 'vector))
 595          (rank (array-rank x))
 596          (dim (make-array rank))
 597          (olddim (coerce (array-dimensions x) 'vector)))
 598     (declare (fixnum rank))
 599     ;; construct new dimension vector
 600     (permute-indices olddim dim perm t)
 601     ;; make result array and the index vectors and lists */
 602     (let* ((result (make-array (coerce dim 'list)))
 603           (indices (make-array rank))
 604           (oldindices (make-array rank))
 605           (ilist (make-list rank))
 606           (data (compound-data-seq x))
 607           (result_data (compound-data-seq result))
 608           (n (length data)))
 609       (declare (fixnum n))
 610       (dotimes (i rank)
 611         (declare (fixnum i))
 612         (setf (aref oldindices i) (list nil)))
 613       ;; fill in the result
 614       (if (/= n (length result_data)) (error "bad data"))
 615       (dotimes (i n result)
 616         (declare (fixnum i))
 617         (let ((k (translate-index i result x perm indices oldindices ilist)))
 618           (declare (fixnum k))
 619           (setf (aref result_data k) (aref data i)))))))
 620
 621 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 622 ;;;;
 623 ;;;;              SUM, PROD, COUNT-ELEMENTS, and MEAN Functions
 624 ;;;;
 625 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 626
 627 (defun sum-1 (x)
 628   (if (numberp x)
 629     x
 630     (let ((seq (compound-data-seq x))
 631           (sum 0))
 632       (if (consp seq)
 633         (dolist (x seq sum)
 634           (setf sum (+ sum (if (numberp x) x (sum-1 x)))))
 635         (let ((n (length seq)))
 636           (declare (fixnum n))
 637           (dotimes (i n sum)
 638             (declare (fixnum i))
 639             (let ((x (aref seq i)))
 640               (setf sum (+ sum (if (numberp x) x (sum-1 x)))))))))))
 641
 642 (defun sum (&rest args)
 643 "Args: (&rest number-data)
 644 Returns the sum of all the elements of its arguments. Returns 0 if there
 645 are no arguments. Vector reducing."
 646   (if args
 647     (sum-1 (if (rest args) args (first args)))
 648     0))
 649
 650 (defun prod-1 (x)
 651   (if (numberp x)
 652     x
 653     (let ((seq (compound-data-seq x))
 654           (prod 1))
 655       (if (consp seq)
 656         (dolist (x seq prod)
 657           (setf prod (* prod (if (numberp x) x (prod-1 x)))))
 658         (let ((n (length seq)))
 659           (declare (fixnum n))
 660           (dotimes (i n prod)
 661             (declare (fixnum i))
 662             (let ((x (aref seq i)))
 663               (setf prod (* prod (if (numberp x) x (prod-1 x)))))))))))
 664
 665 (defun prod (&rest args)
 666 "Args: (&rest number-data)
 667 Returns the product of all the elements of its arguments. Returns 1 if there
 668 are no arguments. Vector reducing."
 669   (if args
 670     (prod-1 (if (rest args) args (first args)))
 671     1))
 672
 673 (defun count-elements (x)
 674 "Args: (number &rest more-numbers)
 675 Returns the number of its arguments. Vector reducing"
 676   (if (compound-data-p x)
 677     (let ((seq (compound-data-seq x))
 678           (count 0))
 679       (if (consp seq)
 680         (dolist (x seq count)
 681           (incf count (if (compound-data-p x) (count-elements x) 1)))
 682         (let ((n (length seq)))
 683           (declare (fixnum n))
 684           (dotimes (i n count)
 685             (declare (fixnum i))
 686             (let ((x (aref seq i)))
 687               (incf count (if (compound-data-p x) (count-elements x) 1)))))))
 688     1))
 689
 690 (defun mean (x)
 691 "Args: (x)
 692 Returns the mean of the elements x. Vector reducing."
 693   (let ((mean 0.0)
 694         (count 0.0))
 695     (labels ((add-to-mean (x)
 696               (let ((count+1 (+ count 1.0)))
 697                 (setf mean (+ (* (/ count count+1) mean) (* (/ count+1) x)))
 698                 (setf count count+1)))
 699              (find-mean (x)
 700                (if (numberp x)
 701                  (add-to-mean x)
 702                  (let ((seq (compound-data-seq x)))
 703                    (if (consp seq)
 704                      (dolist (x seq)
 705                        (if (numberp x) (add-to-mean x) (find-mean x)))
 706                      (let ((n (length seq)))
 707                        (dotimes (i n)
 708                          (declare (fixnum i))
 709                          (let ((x (aref seq i)))
 710                            (if (numberp x)
 711                                (add-to-mean x)
 712                                (find-mean x))))))))))
 713       (find-mean x)
 714       mean)))
 715
 716 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 717 ;;;;
 718 ;;;;                         Sorting Functions
 719 ;;;;
 720 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 721
 722 (defun sort-data (x)
 723 "Args: (sequence)
 724 Returns a sequence with the numbers or strings in the sequence X in order."
 725   (flet ((less (x y) (if (numberp x) (< x y) (string-lessp x y))))
 726     (stable-sort (copy-seq (compound-data-seq x)) #'less)))
 727
 728 (defun order (x)
 729 "Args (x)
 730 Returns a sequence of the indices of elements in the sequence of numbers
 731 or strings X in order."
 732   (let* ((seq (compound-data-seq x))
 733          (type (if (consp seq) 'list 'vector))
 734          (i -1))
 735     (flet ((entry (x) (setf i (+ i 1)) (list x i))
 736            (less (a b)
 737                  (let ((x (first a))
 738                        (y (first b)))
 739                    (if (numberp x) (< x y) (string-lessp x y)))))
 740       (let ((sorted-seq (stable-sort (map type #'entry seq) #'less)))
 741         (map type #'second sorted-seq)))))
 742
 743 ;; this isn't destructive -- do we document destructive only, or any
 744 ;; variant?
 745 (defun rank (x)
 746 "Args (x)
 747 Returns a sequence with the elements of the list or array of numbers or
 748 strings X replaced by their ranks."
 749   (let ((ranked-seq (order (order x))))
 750     (make-compound-data (compound-data-shape x) ranked-seq)))
 751
 752 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 753 ;;;;
 754 ;;;;                    IF-ELSE and SAMPLE Functions
 755 ;;;;
 756 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 757
 758 (defun if-else (a x y)
 759 "Args: (first x y)
 760 Takes simple or compound data items FIRST, X and Y and returns result of
 761 elementswise selecting from X if FIRST is not NIL and from Y otherwise."
 762   (flet ((base-if-else (a x y) (if a x y)))
 763     (recursive-map-elements #'base-if-else #'if-else a x y)))
 764
 765 (defun sample (x ssize &optional replace)
 766 "Args: (x n &optional (replace nil))
 767 Returns a list of a random sample of size N from sequence X drawn with or
 768 without replacement."
 769   (check-sequence x)
 770   (let ((n (length x))
 771         (x (if (consp x) (coerce x 'vector) (copy-vector x)))
 772         (result nil))
 773     (if (< 0 n)
 774         (dotimes (i ssize result)
 775                  (let ((j (if replace (random n) (+ i (random (- n i))))))
 776                    (setf result (cons (aref x j) result))
 777                    (unless replace     ;; swap elements i and j
 778                            (let ((temp (aref x i)))
 779                              (setf (aref x i) (aref x j))
 780                              (setf (aref x j) temp))))))))
 781