;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; This file is part of "NOVA":  NOVA = search + COCOMO tools
; Copyright, 2008, Tim Menzies tim@menzies.us
;
; NOVA is free software: you can redistribute it and/or modify
; it under the terms of the GNU General Public License as published by
; the Free Software Foundation, either version 3 of the License, or
; (at your option) any later version.
;
; NOVA 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 General Public License for more details.
; You should have received a copy of the GNU General Public License
; a long with NOVA.  If not, see <http://www.gnu.org/licenses/>.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(in-package :wvu-lib.tricks)

(defun sum (l)
  "Computes the sum of a list of numbers."
  (let ((sum 0))
    (dolist (x l sum)
      (incf sum x))))

(defun mean (l)
  "Computes the mean of a list of numbers."
  (let ((n   0)
	(sum 0))
    (dolist (x l (/ sum n))
      (incf n)
      (incf sum x))))

(defun stdev-from-list (l)
  "Returns the standard deviation of a list of numbers."
  (let ((n      0)
	(sum    0)
	(sum-sq 0))
    (dolist (x l (stdev n sum sum-sq))
      (incf n)
      (incf sum x)
      (incf sum-sq (* x x)))))

(defun stdev (n sum sum-sq)
  "Computes standard deviation."
  (sqrt (/ (- sum-sq (/ (* sum sum) n)) (- n 1))))

(defun median (l)
  "Computes the median, spread, and quartiles of a list of numbers."
  (labels ((mean (x y) (/ (+ x y) 2)))
    (let* ((sorted    (sort (copy-list l) #'<))
	   (len       (length sorted))
	   (mid-pos   (floor (/ len 2)))
	   (mid-val   (nth mid-pos sorted))
	   (25percent (nth (floor (* len 0.25)) sorted))
	   (75percent (nth (floor (* len 0.75)) sorted))
	   (50percent (if (oddp len)
			  mid-val
			  (mean mid-val (nth (1- mid-pos) sorted)))))
      (values 50percent
	      (- 75percent 50percent)
	      75percent
	      25percent))))

(defun spread (l)
  (multiple-value-bind (median spread 75p 25p) (median l)
    (declare (ignore median 75p 25p))
    spread))

(defun quartiles (l)
  (let ((max (apply #'max l))
        (min (apply #'min l))
        (medvals (multiple-value-list (median l))))
    (values min
            (nth 3 medvals)
            (nth 0 medvals)
            (nth 2 medvals)
            max)))

;;;
;;; ttest stuff
;;;

(defun ttest-from-lists (one two &optional (conf 95))
  "Returns the ttest result from two lists."
  (let ((as  0)
	(asq 0)
	(an  0)
	(bs  0)
	(bsq 0)
	(bn  0))
    (dolist (a one)
      (incf an)
      (incf as a)
      (incf asq (* a a)))
    (dolist (b two)
      (incf bn)
      (incf bs b)
      (incf bsq (* b b)))
    (ttest as asq an bs bsq bn conf)))

(defun ttest (as asq an bs bsq bn &optional (conf 95))
  "Compares means of two independent samples a and b.
     as, bs   = sample sums
     asq, bsq = sample sums of squares
     an, bn   = sample sizes
     conf     = confidence level (95 or 99)"
  (labels
      ((less () (< (/ as an) (/ bs bn)))
       (same () (let* ((tcrit  (tcritical (+ an bn -2) conf))
		       (ssa    (- asq (/ (* as as) an)))
		       (ssb    (- bsq (/ (* bs bs) bn)))
		       (pooled (/ (+ ssa ssb) (+ bn  an -2)))
		       (sxasb  (sqrt (* pooled (+ (/ 1 an) (/ 1 bn)))))
		       (tvalue (abs (/ (- (/ bs bn) (/ as an)) sxasb))))
		  ; (oo tcrit tvalue)
		  (> tcrit tvalue))))
    (cond ((same)  0)    ; H0  : mean of a = mean of b
	  ((less) -1)    ; H1a : mean of a < mean of b
	  (t       1)))) ; H1b : mean of a > mean of b

(let ((ttable '((90 . ((  1 .  6.314  )
		       (  3 .  2.353 )
		       (  5 .  2.015 )
		       ( 10 .  1.812 )
		       ( 20 .  1.725 )
		       ( 80 .  1.67  )
		       (320 .  1.65  )))
                (95 . ((  1 . 12.70  )
		       (  3 .  3.1820)
		       (  5 .  2.5710)
		       ( 10 .  2.2280)
		       ( 20 .  2.0860)
		       ( 80 .  1.99  )
		       (320 .  1.97  )))
		(99 . ((  1 . 63.6570)
		       (  3 .  5.8410)
		       (  5 .  4.0320)
		       ( 10 .  3.1690)
		       ( 20 .  2.8450)
		       ( 80 .  2.64  )
		       (320 .  2.58  ))))))
  (defun tcritical (n conf)
    "Returns the ttest critical values.  Keeps those values as a set of lines,
     and intermediary values are interpolated between the points."
    (multiple-value-bind (points foundp) (geta conf ttable)
      (if foundp          
          (interpolates n points)
          (error "unsupported confidence ~a.  Valid choices: ~a" conf (mapcar #'car ttable))))))


;;; report stuff

(defstruct dist
  (min    most-positive-double-float)
  (max    most-negative-double-float)
  (bars   nil)
  (sum    0)
  (sumsq  0)
  (n      0)
  (sorted t)
  (fuzz   1)
  (fat    nil)
  (all    nil))

(defun dist-compare (d1 d2 &optional (conf 95))
  "Compares the means of two independent samples."
  (ttest (dist-sum d1) (dist-sumsq d1) (dist-n d1)
	 (dist-sum d2) (dist-sumsq d2) (dist-n d2)
	 conf))

(defun as-dist (l)
  "Return a dist with all the numbers in the given list."
  (if (eq 'dist (type-of l))
      l
      (dist-adds l)))

(defun dist-adds (l &optional (d (make-dist)))
  "Add the numbers in the given list to a dist."
  (dolist (x l d)
    (dist-add x d)))

(defun dist-mean (d)
  "Returns the mean of a dist."
  (/ (dist-sum d) (dist-n d)))

(defun dist-sd (d)
  "Returns the standard deviation of a dist."
  (stdev (dist-n d) (dist-sum d) (dist-sumsq d)))

(defun dist-add (num &optional (d (make-dist)))
  "Adds a number to a dist."
  (incf (dist-n     d))
  (incf (dist-sum   d) num)
  (incf (dist-sumsq d) (* num num))
  (if (dist-fat d)
      (push num (dist-all d)))
  (if (< num (dist-min d))
      (setf (dist-min d) num))
  (if (> num (dist-max d))
      (setf (dist-max d) num))
  (let* ((fuzz (dist-fuzz d))
	 (num1 (* fuzz (round (/ num fuzz))))
	 (val  (cdr (assoc num1 (dist-bars d)))))
    (if val
	(setf (cdr (assoc num1 (dist-bars d))) (1+ val))
	(push (cons num1 1) (dist-bars d)))
    (setf (dist-sorted d) nil))
  d)

(defun dist-sort (d)
  "Sorts a dist."
  (labels ((car< (a b) (< (car a ) (car b))))
    (unless (dist-sorted d)
      (setf (dist-bars d) (sort (dist-bars d) #'car<)
	    (dist-sorted d) t))
    d))

(defun distogram (d &key header (lwidth 5) (decimals 2) (rwidth 5)
		  (shrink 1) (on "-") (pad 4) (str t))
  "Prints a dist on stream str, shrinking the right hand side bars by shrink.
   header is some text to show on top.  The other variables control how each
   line is printed."
  (declare (ignore pad))
  (let ((d (as-dist d)))
    (unless (dist-sorted d)
      (dist-sort d))
    (let (;(sum         (dist-sum  d))
	  (bars        (dist-bars d))
	  ;(min         (dist-min  d))
	  ;(max         (dist-max  d))
	  (n           (dist-n    d))
	  (fuzz        (dist-fuzz d))
	  (percent-sum 0)
	  (fmt         (format nil " ~~~a,~af * ~~~ad = ~~3d% ~~a "
			       lwidth decimals rwidth)))
      (if header
	  (format str "~a~%" header))
      (dolist (bar bars t)
	(let* ((key     (car bar))
	       (value   (cdr bar))
	       (percent (floor (* 100 (/ value n))))
	       (stars   (round (/ value shrink)))
	       (halfp   (< percent-sum 50)))
	  (incf percent-sum percent)
	  (format str fmt
		  (+ key (/ fuzz 2))
		  value
		  percent
		  (if halfp "<" ">"))
	  (chars stars on)
	  (terpri str))))))

(defun demo-distogram ()
  "Creates a dist with random numbers and displays it."
  (let ((log (make-dist :fuzz 2)))
    (dotimes (i 1000)
      (dist-add (sqrt (random 100)) log))
    (distogram log :shrink 10)))


(defun quartile-chart (absolute-min min 25p med 75p max absolute-max &key (total-width 20) (stream t))
  (let ((range (- absolute-max absolute-min)))
    (labels ((calc-pos (val)
               (min (floor (* (/ val range)
                              total-width))
                    (1- total-width))))
      (let* ((amin-pos 0)
             (min-pos (calc-pos min))
             (25p-pos (calc-pos 25p))
             (med-pos (calc-pos med))
             (75p-pos (calc-pos 75p))
             (max-pos (calc-pos max))
             (amax-pos (1- total-width))
             (chart (make-array (list total-width) :element-type 'character :initial-element  #\  )))
        ;; first and fourth quartiles
        (loop for pos from min-pos to max-pos
           do (setf (aref chart pos ) #\-))
        ;; second and third quartiles
        (loop for pos from 25p-pos to 75p-pos
           do (setf (aref chart pos ) #\=))
        ;; set min and max
        (setf (aref chart min-pos) #\|)
        (setf (aref chart max-pos) #\|)
        ;; median
        (setf (aref chart med-pos) #\#)
        ;; absolute min and absolute max
        (setf (aref chart amin-pos) #\[)
        (setf (aref chart amax-pos) #\])
        ;; printx
        (format stream "~a" (coerce chart 'string))))))

;;; ranking stats stuff

(defun showh (h)
  "Prints hash table entries sorted by key."
  (let (l)
    (maphash #'(lambda (k v) (push (cons k v) l)) h)
    (mapc #'print (sort l #'(lambda (a b) (< (car a) (car b)))))))

(defun gen-rank-ht (l &optional (ranks (make-hash-table)) (n 0))
  (if (null l)
      ranks
      (let (repeats sum now)
        (labels ((walk () (incf n) (pop l))
                 (new  () (setf repeats 1) (setf sum n))
                 (same () (incf sum n) (incf repeats))
                 (spin () (when (eql now (car l))
                            (walk) (same) (spin))))
          (setf now (walk))
          (new)
          (spin)
          (setf (gethash now ranks) (/ sum repeats))
          (gen-rank-ht l ranks n)))))

(defun rank (list-of-numbers &optional (sort-fn #'>))
  "Returns a list of the elements of l ranked.  Ties are given their average rank"
  (let ((rank-ht (gen-rank-ht (sort (copy-list list-of-numbers) sort-fn))))
    (mapcar #'(lambda (x) (gethash x rank-ht)) list-of-numbers)))

(defun mann-whitney-demo ()
  "Performs a Mann-Whitney test."
  (mann-whitney '(4.6 4.7 4.9 5.1 5.2 5.5 5.8 6.1 6.5 6.5 7.2)
		'(5.2 5.3 5.4 5.6 6.2 6.3 6.8 7.7 8.0 8.1)))

(defun mann-whitney-demo-big (&optional (fudge 1))
  "Generates two lists of 10,000 random floats.  Multiplies the second list by
   fudge.  Performs a Mann-Whitney test on the two lists."
  (labels ((big (n s) (let (out) (dotimes (i n out)
				   (push (* s (random 100)) out)))))
    (let ((one (big 10000 1))
	  (two (big 10000 fudge)))
      (time (mann-whitney one two)))))

(defun mann-whitney (a b &optional (conf 95))
  "Performs a Mann-Whitney test as described in method 1 of
   http://faculty.vassar.edu/lowry/ch11a.html."
  (labels ((as-ranks (l r) (mapcar #'(lambda (x) (gethash x r)) l))
	   (sum      (l)   (let ((s 0)) (dolist (x l s) (incf s x))))
	   (median   (l n) (let ((sorted (sort l #'<))
				 (midv (floor (/ n 2))))
			     (if (oddp n)
				 (nth midv sorted)
				 (/ (+ (nth (1- midv) sorted)
				       (nth midv sorted)) 2)))))
    (let* ((all    (sort (copy-list (append a b)) #'<))
	   (ranks  (gen-rank-ht all))
	   (ranksa (as-ranks a ranks))
	   (ranksb (as-ranks b ranks))
	   (na     (length a))
	   (nb     (length b))
	   (n      (+ na nb))
	   (tcrit  (tcritical n conf))
	   (suma   (* 1.0 (sum ranksa)))
	   (ta     (/ (* na (+ n 1)) 2.0))
	   (sigma  (sqrt (/  (* na nb (+ n 1)) 12.0)))
	   (za     (/ (+ (- suma ta) 0.5) sigma)))
      ; (oo suma ta sigma tcrit)
      (cond  ((< (abs za) tcrit)                         0)
	     ((< (median ranksa na) (median ranksb nb)) -1)
	     (t                                          1)))))

(defun mann-whitney-compare (l1 l2 combine-fn compare-fn &optional (conf 95))
  "compare l1 and l2, if mann-whitney can't prove one is different,
compare by comparing result of calling combine-fn on l1 and l2 (cl1 and cl2).

l1 and l2; each are a sequence of numbers.
combine-fn should return a value usable by compare-fn
compare-fn with cl1 cl2 should return t if cl1 comes before cl2,
if cl2 should come before cl1 or they could go either way, return nil."
  (macrolet ((combine (l) `(funcall combine-fn ,l))
             (compare (cl1 cl2) `(funcall compare-fn ,cl1 ,cl2)))
    (declare (sequence l1 l2)
             (function combine-fn compare-fn)
             (number conf))
    (let ((mw-result (mann-whitney l1 l2 conf)))
      (if (= mw-result 0)
          mw-result
          (let ((cl1 (combine l1))
                (cl2 (combine l2)))
            (cond ((not (or (compare cl1 cl2)
                            (compare cl2 cl1)))
                   0)
                  ((compare cl1 cl2)
                   1)
                  ((compare cl2 cl1)
                   -1)))))))



;;FROM DEMSAR
(defun friedman-test (data &optional (confidence 95))
  "Performs a Friedman test as described in demsar06
   data is a list of lists, N lists of k elements where
   each inner list refers to a row."
  ;;rank the rows
  (let* ((N (length data))
         (k (length (first data)))
         (ranked-data (let (l)
                        (dolist (row data (nreverse l))
                          (push (rank row) l))))
         (ranked-col-means (mapcar #'mean (transpose ranked-data)))
         (friedman-stat (* (/ (* 12 N)
                              (* k (1+ k)))
                           (- (reduce #'+
                                      (mapcar #'(lambda (x) (expt x 2))
                                              ranked-col-means))
                              (/ (* k
                                    (expt (1+ k)
                                          2))
                                 4))))
         (fudged-friedman-stat (/ (* (1- N)
                                     friedman-stat)
                                  (- (* N
                                        (1- k))
                                     friedman-stat))))
    ;;return t if different
    (>= fudged-friedman-stat (f-value k N confidence))))

(defun nemenyi-critical-difference (data &optional (confidence 95))
  "Calculates nemenyi critical difference as described in demsar06
   data is a list of lists, N lists of k elements where
   each inner list refers to a row."
  (let* ((N (length data))
         (k (length (first data)))
         (ranked-data (let (l)
                        (dolist (row data (nreverse l))
                          (push (rank row) l))))
         (ranked-col-means (mapcar #'mean (transpose ranked-data)))
         (critical-difference (* (nemenyi-critical-value k confidence)
                                 (sqrt (/ (* k (1+ k))
                                          (* 6 N))))))
    (values critical-difference ranked-col-means)))