#lang racket/base
;;; Racket Simulation Collection
;;; statistics.rkt
;;; Copyright (c) 2004-2010 M. Douglas Williams
;;;
;;; This file is part of the Racket Simulation Collection.
;;;
;;; The Racket Simulation Collection 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 3 of the License,
;;; or (at your option) any later version.
;;;
;;; The Racket Simulation Collection 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 the Racket Simulation Collection.  If not, see
;;; <http://www.gnu.org/licenses/>.
;;;
;;; -----------------------------------------------------------------------------
;;;
;;; Version  Date      Description
;;; 1.0.0    07/17/05  Cleaned up initialization code. Added time parameter to
;;;                    the constructor. (MDW)
;;; 1.0.1    05/08/06  Removed time-last-synchronized field. This is maintained
;;;                    at the variable level. (MDW)
;;; 1.0.2    11/20/07  Added vector statistics. (MDW)
;;; 3.0.0    06/24/08  Updated for V4.0. (MDW)
;;; 3.0.1    11/27/08  Converted to a module. (MDW)
;;; 4.0.0    08/15/10  Converted to Racket. (MDW)

(struct statistics (time-dependant?
                    minimum
                    maximum
                    n
                    sum
                    sum-of-squares
                    dimension)
  #:mutable)

;;; (make-statistics time-dependant? time) -> statistics?
;;;   time-dependant? : boolean?
;;;   time : (>=/c 0.0)
(define (make-statistics time-dependant? time)
  (statistics time-dependant? +inf.0 -inf.0 0 0.0 0.0 0))

;;; (make-vector-statistics time-dependant? time dimension) -> statistics
;;;   time-dependant? : boolean?
;;;   time : (>=/c 0.0)
;;;   dimension : exact-positive-integer?
(define (make-vector-statistics time-dependant? time dimension)
  (statistics 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 value time) -> void?
;;;   statistics : statistics?
;;;   value : real?
;;;   time : (>=/c 0.0)
;;; 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 value) -> void?
;;;   statistics : statistics?
;;;   value : real?
;;; Tally the statistics for the value.
(define (statistics-tally! statistics value)
  (statistics-accumulate! statistics value 1))

;;; Pseudo slots

;;; (statistics-mean statistics) -> real?
;;;   statistics : statistics?
(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))))))
    
;;; (statistics-mean-square statistics) -? real?
;;;   statistics : 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))))))

;;; (statistics-variance statistics) -> real?
;;;   statistics : 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))))))))

;;; (statistics-standard-deviation statistics) -> real?
;;;   statistics : 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)))))))

;;; Module Contracts

(provide (all-defined-out))