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