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