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