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 ;;; Package Setup
  12
  13 (defpackage :lisp-stat-basics
  14     (:use :common-lisp
  15           :lisp-stat-object-system
  16           :lisp-stat-types
  17           :lisp-stat-fastmap
  18           :lisp-stat-float
  19           :lisp-stat-macros
  20           :lisp-stat-compound-data
  21           ;;:lisp-stat-matrix
  22           ;;:lisp-stat-linalg
  23           :lisp-stat-probability)
  24   (:shadowing-import-from :lisp-stat-object-system
  25                           slot-value call-method call-next-method)
  26   (:export
  27    permute-array sum prod count-elements mean if-else sample ))
  28
  29
  30
  31 (in-package #:lisp-stat-basics)
  32
  33
  34
  35
  36 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  37 ;;;;
  38 ;;;;                    Array Permutation Functions
  39 ;;;;
  40 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  41
  42 (defun permute-indices (x y perm check)
  43   "Args: (x y perm check).
  44 permute x into y using perm; all should be vectors; If check is TRUE
  45 the routine will check to make sure no indices are reused, but x
  46 will be destroyed."
  47   (let ((rank (length x)))
  48     (declare (fixnum rank))
  49     (dotimes (i rank)
  50       (declare (fixnum i))
  51       (let ((k (aref perm i)))
  52         (if (not (fixnump k)) (error "bad permutation sequence - ~a" perm))
  53         (if (or (< k 0) (>= k rank))
  54             (error "bad permutation sequence - ~a" perm))
  55         (setf (aref y i) (aref x k))
  56         ;; to insure dimensions are not re-used
  57         (if check (setf (aref x k) NIL))))))
  58
  59 (defun indices-from-rowmajor (a k result)
  60   "Args: (a k result).
  61 Compute indices in a from rowmajor index k, put in vector result."
  62   (declare (fixnum k))
  63
  64   (if (not (arrayp a)) (error "not an array - ~a" a))
  65   (if (or (> 0 k) (>= k (array-total-size a))) (error "index out of range"))
  66
  67   (let ((face 1)
  68         (rank (array-rank a))
  69         (dim (array-dimensions a)))
  70     (declare (fixnum face rank))
  71     (let ((cdim (make-next-element dim)))
  72       (dotimes (i rank)
  73         (declare (fixnum i))
  74         (setf face (* face (get-next-element cdim i)))))
  75     (let ((cdim (make-next-element dim)))
  76       (dotimes (i rank)
  77         (setf face (/ face (get-next-element cdim i)))
  78         (setf (aref result i) (floor (/ k face)))
  79         (setf k (rem k face))))))
  80
  81 (defun translate-index (i result x perm indices oldindices ilist)
  82   "Args: (i result x perm indices oldindices ilist).
  83 Translate row major index in original array to row major index in new
  84 array. Use indices vectors and ilist for temporary storage."
  85   (declare (fixnum i))
  86   (let ((rank (array-rank x)))
  87     (declare (fixnum rank))
  88     (indices-from-rowmajor x i oldindices)
  89     (permute-indices oldindices indices perm nil)
  90     (do ((next ilist (rest next))
  91          (k 0 (+ k 1)))
  92         ((not (and (< k rank) (consp next))))
  93       (setf (first next) (aref indices k)))
  94     (apply #'array-row-major-index result ilist)))
  95
  96 (defun permute-array (x perm)
  97   "Args: (a p)
  98 Returns a copy of the array A permuted according to the permutation P."
  99   (if (not (arrayp x)) (error "not an array - ~a" x))
 100   (check-sequence perm)
 101   (if (/= (length perm) (array-rank x))
 102     (error "bad permutation sequence - ~a" perm))
 103   (let* ((perm (coerce perm 'vector))
 104          (rank (array-rank x))
 105          (dim (make-array rank))
 106          (olddim (coerce (array-dimensions x) 'vector)))
 107     (declare (fixnum rank))
 108     ;; construct new dimension vector
 109     (permute-indices olddim dim perm t)
 110     ;; make result array and the index vectors and lists */
 111     (let* ((result (make-array (coerce dim 'list)))
 112           (indices (make-array rank))
 113           (oldindices (make-array rank))
 114           (ilist (make-list rank))
 115           (data (compound-data-seq x))
 116           (result_data (compound-data-seq result))
 117           (n (length data)))
 118       (declare (fixnum n))
 119       (dotimes (i rank)
 120         (declare (fixnum i))
 121         (setf (aref oldindices i) (list nil)))
 122       ;; fill in the result
 123       (if (/= n (length result_data)) (error "bad data"))
 124       (dotimes (i n result)
 125         (declare (fixnum i))
 126         (let ((k (translate-index i result x perm indices oldindices ilist)))
 127           (declare (fixnum k))
 128           (setf (aref result_data k) (aref data i)))))))
 129
 130 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 131 ;;;;
 132 ;;;;              SUM, PROD, COUNT-ELEMENTS, and MEAN Functions
 133 ;;;;
 134 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 135
 136 (defun sum-1 (x)
 137   (if (numberp x)
 138     x
 139     (let ((seq (compound-data-seq x))
 140           (sum 0))
 141       (if (consp seq)
 142         (dolist (x seq sum)
 143           (setf sum (+ sum (if (numberp x) x (sum-1 x)))))
 144         (let ((n (length seq)))
 145           (declare (fixnum n))
 146           (dotimes (i n sum)
 147             (declare (fixnum i))
 148             (let ((x (aref seq i)))
 149               (setf sum (+ sum (if (numberp x) x (sum-1 x)))))))))))
 150
 151 (defun sum (&rest args)
 152 "Args: (&rest number-data)
 153 Returns the sum of all the elements of its arguments. Returns 0 if there
 154 are no arguments. Vector reducing."
 155   (if args
 156     (sum-1 (if (rest args) args (first args)))
 157     0))
 158
 159 (defun prod-1 (x)
 160   (if (numberp x)
 161     x
 162     (let ((seq (compound-data-seq x))
 163           (prod 1))
 164       (if (consp seq)
 165         (dolist (x seq prod)
 166           (setf prod (* prod (if (numberp x) x (prod-1 x)))))
 167         (let ((n (length seq)))
 168           (declare (fixnum n))
 169           (dotimes (i n prod)
 170             (declare (fixnum i))
 171             (let ((x (aref seq i)))
 172               (setf prod (* prod (if (numberp x) x (prod-1 x)))))))))))
 173
 174 (defun prod (&rest args)
 175 "Args: (&rest number-data)
 176 Returns the product of all the elements of its arguments. Returns 1 if there
 177 are no arguments. Vector reducing."
 178   (if args
 179     (prod-1 (if (rest args) args (first args)))
 180     1))
 181
 182 (defun count-elements (x)
 183 "Args: (number &rest more-numbers)
 184 Returns the number of its arguments. Vector reducing"
 185   (if (compound-data-p x)
 186     (let ((seq (compound-data-seq x))
 187           (count 0))
 188       (if (consp seq)
 189         (dolist (x seq count)
 190           (incf count (if (compound-data-p x) (count-elements x) 1)))
 191         (let ((n (length seq)))
 192           (declare (fixnum n))
 193           (dotimes (i n count)
 194             (declare (fixnum i))
 195             (let ((x (aref seq i)))
 196               (incf count (if (compound-data-p x) (count-elements x) 1)))))))
 197     1))
 198
 199
 200 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 201 ;;;;
 202 ;;;;                    IF-ELSE Functions
 203 ;;;;
 204 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 205
 206 (defun if-else (a x y)
 207 "Args: (first x y)
 208 Takes simple or compound data items FIRST, X and Y and returns result of
 209 elementswise selecting from X if FIRST is not NIL and from Y otherwise."
 210   (flet ((base-if-else (a x y) (if a x y)))
 211     (recursive-map-elements #'base-if-else #'if-else a x y)))