external/clem/src/statistics.lisp

   1 ;;; statistics.lisp
   2 ;;; macros, functions and methods for matrix element access
   3 ;;;
   4 ;;; Copyright (c) 2004-2006 Cyrus Harmon (ch-lisp@bobobeach.com)
   5 ;;; All rights reserved.
   6 ;;;
   7 ;;; Redistribution and use in source and binary forms, with or without
   8 ;;; modification, are permitted provided that the following conditions
   9 ;;; are met:
  10 ;;;
  11 ;;;   * Redistributions of source code must retain the above copyright
  12 ;;;     notice, this list of conditions and the following disclaimer.
  13 ;;;
  14 ;;;   * Redistributions in binary form must reproduce the above
  15 ;;;     copyright notice, this list of conditions and the following
  16 ;;;     disclaimer in the documentation and/or other materials
  17 ;;;     provided with the distribution.
  18 ;;;
  19 ;;; THIS SOFTWARE IS PROVIDED BY THE AUTHOR 'AS IS' AND ANY EXPRESSED
  20 ;;; OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  21 ;;; WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  22 ;;; ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
  23 ;;; DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  24 ;;; DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
  25 ;;; GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  26 ;;; INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
  27 ;;; WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  28 ;;; NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  29 ;;; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  30 ;;;
  31
  32 (in-package :clem)
  33
  34 ;;; slow fallback methods
  35
  36 (defmethod sum ((m matrix))
  37   (destructuring-bind (mr mc) (dim m)
  38     (sum-range m 0 (- mr 1) 0 (- mc 1))))
  39
  40 (defmethod sum-cols ((m matrix) &key (matrix-class (class-of m)))
  41   (let ((mr (rows m)) (mc (cols m)))
  42     (let ((n (make-instance matrix-class :rows 1 :cols mc)))
  43       (dotimes (i mr)
  44         (dotimes (j mc)
  45           (incf (mref n 0 j) (mref m i j))))
  46       n)))
  47
  48 (defmethod sum-rows ((m matrix) &key (matrix-class (class-of m)))
  49   (let ((mr (rows m)) (mc (cols m)))
  50     (let ((n (make-instance matrix-class :rows mr :cols 1)))
  51       (dotimes (i mr)
  52         (dotimes (j mc)
  53           (incf (mref n i 0) (mref m i j))))
  54       n)))
  55
  56 (defmethod sum-square-range ((m matrix) (startr fixnum) (endr fixnum) (startc fixnum) (endc fixnum))
  57   (declare (dynamic-extent startr endr startc endc)
  58            (fixnum startr endr startc endc))
  59   (let ((acc 0))
  60     (map-range m startr endr startc endc
  61                #'(lambda (v i j)
  62                    (declare (ignore i j))
  63                    (incf acc (* v v))))
  64     acc))
  65
  66 (defmethod sum-square ((m matrix))
  67   (destructuring-bind (mr mc) (dim m)
  68     (sum-square-range m 0 (- mr 1) 0 (- mc 1))))
  69
  70 (defmethod sum-range ((m matrix) (startr fixnum) (endr fixnum) (startc fixnum) (endc fixnum))
  71   (declare (dynamic-extent startr endr startc endc)
  72            (fixnum startr endr startc endc))
  73   (let ((acc 0))
  74     (map-range m startr endr startc endc
  75                #'(lambda (v i j)
  76                    (declare (ignore i j))
  77                    (incf acc v)))
  78     acc))
  79
  80 (defmethod mean-range ((m matrix) startr endr startc endc)
  81   (double-float-divide (sum-range m startr endr startc endc)
  82                 (count-range startr endr startc endc)))
  83
  84 (defmethod mean ((m matrix))
  85   (destructuring-bind (mr mc) (dim m)
  86     (mean-range m 0 (- mr 1) 0 (- mc 1))))
  87
  88 (defmethod variance-range ((m matrix) startr endr startc endc)
  89   (declare (dynamic-extent startr endr startc endc)
  90            (fixnum startr endr startc endc))
  91   (let ((mu (mean-range m startr endr startc endc)))
  92     (let ((musq (* mu mu)))
  93       (let ((ssr (sum-square-range m startr endr startc endc)))
  94         (let ((cr (count-range startr endr startc endc)))
  95           (declare (fixnum cr))
  96           (- (double-float-divide ssr cr)
  97              musq))))))
  98
  99 (defmethod variance ((m matrix))
 100   (destructuring-bind (mr mc) (dim m)
 101     (variance-range m 0 (- mr 1) 0 (- mc 1))))
 102
 103 (defmethod sample-variance-range ((m matrix) startr endr startc endc)
 104   (let* ((acc 0)
 105          (mu (mean-range m startr endr startc endc))
 106          (musq (* mu mu)))
 107     (map-range m startr endr startc endc
 108                #'(lambda (v i j)
 109                    (declare (ignore i j))
 110                    (incf acc (- (* v v) musq))))
 111     (double-float-divide acc (1- (count-range startr endr startc endc)))))
 112
 113 (defmethod sample-variance ((m matrix))
 114   (destructuring-bind (mr mc) (dim m)
 115     (sample-variance-range m 0 (- mr 1) 0 (- mc 1))))
 116
 117 ;;; faster type-specific methods
 118
 119 (defmacro %%sum-range (m startr endr startc endc element-type accumulator-type)
 120   `(let ((acc (coerce 0 ',accumulator-type))
 121          (a (matrix-vals ,m)))
 122      (declare (type ,accumulator-type acc)
 123               (type (simple-array ,element-type *) a))
 124      (do ((i ,startr (1+ i)))
 125          ((> i ,endr))
 126        (declare (dynamic-extent i) (type fixnum i))
 127        (do ((j ,startc (1+ j)))
 128            ((> j ,endc))
 129          (declare (dynamic-extent j) (type fixnum j))
 130          (setf acc (+ acc (aref a i j)))))
 131      acc))
 132
 133 (macrolet
 134     ((frob-sum-range (matrix-type accumulator-type)
 135        (let ((element-type (element-type (find-class matrix-type))))
 136          `(defmethod sum-range ((m ,matrix-type)
 137                                 (startr fixnum) (endr fixnum)
 138                                 (startc fixnum) (endc fixnum))
 139             (%%sum-range m startr endr startc endc
 140                          ,element-type ,accumulator-type)))))
 141   (frob-sum-range double-float-matrix double-float)
 142   (frob-sum-range single-float-matrix single-float)
 143
 144   (frob-sum-range ub8-matrix (unsigned-byte 32))
 145   (frob-sum-range ub16-matrix (unsigned-byte 32))
 146   (frob-sum-range ub32-matrix (unsigned-byte 32))
 147
 148   (frob-sum-range sb8-matrix (signed-byte 32))
 149   (frob-sum-range sb16-matrix (signed-byte 32))
 150   (frob-sum-range sb32-matrix (signed-byte 32))
 151
 152   (frob-sum-range fixnum-matrix (signed-byte 32))
 153   (frob-sum-range bit-matrix (signed-byte 32)))
 154
 155
 156 (macrolet
 157     ((frob-sum-square-range (matrix-type accumulator-type)
 158        (let ((element-type (element-type (find-class matrix-type))))
 159          `(defmethod sum-square-range ((m ,matrix-type) (startr fixnum) (endr fixnum) (startc fixnum) (endc fixnum))
 160             (let ((acc (coerce 0 ',accumulator-type))
 161                   (a (matrix-vals m)))
 162               (declare (type ,accumulator-type acc)
 163                        (type (simple-array ,element-type *) a))
 164               (do ((i startr (1+ i)))
 165                   ((> i endr))
 166                 (declare (dynamic-extent i) (type fixnum i))
 167                 (do ((j startc (1+ j)))
 168                     ((> j endc))
 169                   (declare (dynamic-extent j) (type fixnum j))
 170                   (incf acc (* (aref a i j) (aref a i j)))))
 171               acc)))))
 172
 173   (frob-sum-square-range double-float-matrix double-float)
 174   (frob-sum-square-range single-float-matrix single-float)
 175
 176   (frob-sum-square-range ub8-matrix (unsigned-byte 32))
 177   (frob-sum-square-range ub16-matrix (unsigned-byte 32))
 178   (frob-sum-square-range ub32-matrix (unsigned-byte 32))
 179
 180   (frob-sum-square-range sb8-matrix (signed-byte 32))
 181   (frob-sum-square-range sb16-matrix (signed-byte 32))
 182   (frob-sum-square-range sb32-matrix (signed-byte 32))
 183
 184   (frob-sum-square-range fixnum-matrix (signed-byte 32))
 185   (frob-sum-square-range bit-matrix (signed-byte 32)))
 186