;;; PLT Scheme Simulation Collection
;;; statistics.ss
;;; Copyright (c) 2005-2008 M. Douglas Williams
;;;
;;; This library is free software; you can redistribute it and/or
;;; modify it under the terms of the GNU Lesser General Public
;;; License as published by the Free Software Foundation; either
;;; version 2.1 of the License, or (at your option) any later version.
;;;
;;; This library is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
;;; Lesser General Public License for more details.
;;;
;;; You should have received a copy of the GNU Lesser General Public
;;; License along with this library; if not, write to the Free
;;; Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
;;; 02111-1307 USA.
;;;
;;; Version  Date      Description
;;; 1.0.0    07/17/05  Cleaned up initialization code.  Added time parameter
;;;                    to the constructor.  (Doug Williams)
;;; 1.0.1    05/08/06  Removed time-last-synchronized field.  This is maintained
;;;                    at the variable level.  (Doug Williams)
;;; 1.0.2    11/20/07  Added vector statistics.  (Doug Williams)
;;; 3.0.0    06/24/08  Updated for V4.0.  (Doug Williams)

;;; statistics structure
;;;
;;;  #  slot                    description
;;;  -  ----                    -----------
;;;  0  time-dependant?         #t if the statistics are time dependant
;;;  1  minimum                 minimum value
;;;  2  maximum                 maximum value
;;;  3  n                       number of values
;;;  4  sum                     sum of the values
;;;  5  sum-of-squares          sum of the squares of the values
;;;  6  dimension               0 for scalar, > 0 for vector
(define-values (struct:statistics
                statistics-constructor
                statistics?
                statistics-field-ref
                set-statistics-field!)
  (make-struct-type 'statistics #f 7 0))

;;; statistics structure, time-dependant? slot
;;; statistics-time-dependant?: statistics -> boolean
;;; set-statistics-time-dependant?!: statistics x boolean -> void
(define statistics-time-dependant?
  (make-struct-field-accessor statistics-field-ref 0 'time-dependant?))
(define set-statistics-time-dependant?!
  (make-struct-field-mutator set-statistics-field! 0 'time-dependant?))

;;; statistics structure, minimum slot
;;; statistics-minimum: statistics -> real
;;; set-statistics-minimum!: statistics x real -> void
(define statistics-minimum
  (make-struct-field-accessor statistics-field-ref 1 'minimum))
(define set-statistics-minimum!
  (make-struct-field-mutator set-statistics-field! 1 'minimum))

;;; statistics structure, maximum slot
;;; statistics-maximum: statistics -> real
;;; set-statistics-maximum!: statistics x real -> void
(define statistics-maximum
  (make-struct-field-accessor statistics-field-ref 2 'maximum))
(define set-statistics-maximum!
  (make-struct-field-mutator set-statistics-field! 2 'maximum))

;;; statistics structure, n slot
;;; statistics-n: statistics -> real
;;; set-statistics-n!: statistics x real -> void
(define statistics-n
  (make-struct-field-accessor statistics-field-ref 3 'n))
(define set-statistics-n!
  (make-struct-field-mutator set-statistics-field! 3 'n))

;;; statistics structure, sum slot
;;; statistics-sum: statistics -> real
;;; set-statistics-sum!: statistics x real -> void
(define statistics-sum
  (make-struct-field-accessor statistics-field-ref 4 'sum))
(define set-statistics-sum!
  (make-struct-field-mutator set-statistics-field! 4 'sum))

;;; statistics structure, sum-of-squares slot
;;; statistics-sum-of-squares: statistics -> real
;;; set-statistics-sum-of-squares!: statistics x real -> void
(define statistics-sum-of-squares
  (make-struct-field-accessor statistics-field-ref 5 'sum-of-squares))
(define set-statistics-sum-of-squares!
  (make-struct-field-mutator set-statistics-field! 5 'sum-of-squares))

;;; statistics structure, dimension slot
;;; statistics-dimension: statistics -> natural
;;; set-statistics-dimension!: statistics x integer -> void
(define statistics-dimension
  (make-struct-field-accessor statistics-field-ref 6 'dimension))
(define set-statistics-dimension!
  (make-struct-field-mutator set-statistics-field! 6 'dimension))

;;; make-statistics: boolean x real -> statistics
(define (make-statistics time-dependant? time)
  (statistics-constructor time-dependant? +inf.0 -inf.0 0 0.0 0.0 0))

;;; make-vector-statistics: boolean x real x positive -> statistics
(define (make-vector-statistics time-dependant? time dimension)
  (statistics-constructor time-dependant?
                          (make-vector dimension +inf.0)
                          (make-vector dimension -inf.0)
                          0
                          (make-vector dimension 0.0)
                          (make-vector dimension 0.0)
                          dimension))

;;; statistics-accumulate!: statistics x real x real -> void
;;; Accumulate the statistics for the value for the specified time
;;; (actually, duration).
(define (statistics-accumulate! statistics value time)
  (set-statistics-n!
   statistics (+ (statistics-n statistics) time))
  (if (= (statistics-dimension statistics) 0)
      (begin
        (when (< value (statistics-minimum statistics))
          (set-statistics-minimum! statistics value))
        (when (> value (statistics-maximum statistics))
          (set-statistics-maximum! statistics value))
        (let ((weighted-value (* value time))
              (weighted-value-squared (* value value time)))
          (set-statistics-sum!
           statistics (+ (statistics-sum statistics) weighted-value))
          (set-statistics-sum-of-squares! 
           statistics (+ (statistics-sum-of-squares statistics)
                         weighted-value-squared))))
      (do ((i 0 (+ i 1)))
          ((> i (statistics-dimension statistics)) (void))
        (when (< (vector-ref value i) (vector-ref (statistics-minimum statistics) i))
          (vector-set! (statistics-minimum statistics) i (vector-ref value i)))
        (when (> (vector-ref value i) (vector-ref (statistics-maximum statistics) i))
          (vector-set! (statistics-maximum statistics) i (vector-ref value i)))
        (let ((weighted-value (* (vector-ref value i) time))
              (weighted-value-squared (* (vector-ref value i)
                                       (vector-ref value i)
                                       time)))
          (vector-set! (statistics-sum statistics) i
                       (+ (vector-ref (statistics-sum statistics) i)
                          weighted-value))
          (vector-set! (statistics-sum-of-squares statistics) i
                       (+ (vector-ref (statistics-sum-of-squares statistics) i)
                          weighted-value-squared))))))

;;; statistics-tally!: statistics x real -> void
;;; Tally the statistics for the value.
(define (statistics-tally! statistics value)
  (statistics-accumulate! statistics value 1))

;;; Pseudo slots

;;; statistics-mean: statistics -> real
(define (statistics-mean statistics)
  (if (= (statistics-dimension statistics) 0)
      (/ (statistics-sum statistics)
         (statistics-n statistics))
      (let ((n (statistics-n statistics))
            (sum (statistics-sum statistics))
            (mean (make-vector (statistics-dimension statistics))))
        (do ((i 0 (+ i 1)))
            ((> i (statistics-dimension statistics)) (mean))
          (vector-set! mean i
                       (/ (vector-ref sum i) n))))))
                               
;(define (statistics-mean-square statistics)
;  (/ (statistics-sum-of-squares statistics)
;     (statistics-n statistics)))
(define (statistics-mean-square statistics)
  (if (= (statistics-dimension statistics) 0)
      (/ (statistics-sum-of-squares statistics)
         (statistics-n statistics))
      (let ((n (statistics-n statistics))
            (sum-of-squares (statistics-sum-of-squares statistics))
            (mean-square (make-vector (statistics-dimension statistics))))
        (do ((i 0 (+ i 1)))
            ((> i (statistics-dimension statistics)) (mean-square))
          (vector-set! mean-square i
                       (/ (vector-ref sum-of-squares i) n))))))

;(define (statistics-variance statistics)
;  (- (statistics-mean-square statistics)
;     (* (statistics-mean statistics)
;        (statistics-mean statistics))))
(define (statistics-variance statistics)
  (if (= (statistics-dimension statistics) 0)
      (- (statistics-mean-square statistics)
         (* (statistics-mean statistics)
            (statistics-mean statistics)))
      (let ((n (statistics-n statistics))
            (mean (statistics-mean statistics))
            (mean-square (statistics-mean-square statistics))
            (variance (make-vector (statistics-dimension statistics))))
        (do ((i 0 (+ i 1)))
            ((> i (statistics-dimension statistics)) (variance))
          (vector-set! variance i
                       (- (vector-ref mean-square i)
                          (* (vector-ref mean i)
                             (vector-ref mean i))))))))

;(define (statistics-standard-deviation statistics)
;  (sqrt (statistics-variance statistics)))
(define (statistics-standard-deviation statistics)
  (if (= (statistics-dimension statistics) 0)
      (sqrt (statistics-variance statistics))
      (let ((variance (statistics-variance statistics))
            (standard-deviation (make-vector (statistics-dimension statistics))))
        (do ((i 0 (+ i 1)))
            ((> i (statistics-dimension statistics)) (standard-deviation))
          (vector-set! standard-deviation i
                       (sqrt (vector-ref variance i)))))))