lsbasics.lsp

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