;;;;Peter Santiago, Started July 8th 2008
;;;;compiled under:
;;;;slime 1:20070927-3
;;;;SBCL  1:1.0.11.0-1

(load "utests.lisp") ;loads unit-testing file

;;;;this check is for funsies and isn't actually relevant to the whole thing
(deftest friedman-test()
  "runs friedman on data with known results.
   Note: The demsar paper uses the wrong ranks, so there
   is a slight difference"
  (check (= 299/75 (friedman '((0.763 0.768 0.771 0.798)
                          (0.599 0.591 0.590 0.569)
                          (0.954 0.971 0.968 0.967)
                          (0.628 0.661 0.654 0.657)
                          (0.882 0.888 0.886 0.898)
                          (0.936 0.931 0.916 0.931)
                          (0.661 0.668 0.609 0.685)
                          (0.583 0.583 0.563 0.625)
                          (0.775 0.838 0.866 0.875)
                          (1.000 1.000 1.000 1.000)
                          (0.940 0.962 0.965 0.962)
                          (0.619 0.666 0.614 0.669)
                          (0.972 0.981 0.975 0.975) ;demsar pg 14 is wrong here
                          (0.957 0.978 0.946 0.970))))))

(defun friedman (data)
  "each row of data represents a specific dataset,
   while the column represents the test run on that data."
  (labels ((as-ranks (l r) (mapcar #'(lambda (x) (gethash x r)) l)))
    (let* ((templist 0) ;this is needed as sort is destructive
           (data.rows (list-length data))
           (data.cols (list-length (first data)))
           (e.avg (/ (loop for item from 1 to data.cols
                           summing item) data.cols)) ;expected average
           (ranklist '()) ;holds averaged ranks
           (f.stat 0)
           (f.f 0))

      (dotimes (i data.rows)
        (setf templist (copy-list (nth i data)))
        (setf (nth i data) (as-ranks templist
                                     (rank (sort (nth i data) #'>)))))
      ;;;prettifully print the ranks
      ;(dotimes (i (list-length data))
       ;  (format t "~a~%" (nth i data)))

      ;;this loop is destructive to data!
      (dotimes (i data.cols)
        (setf ranklist (append ranklist  `(,(/ (loop for item in
                                                     (mapcar #'car data)
                                                     summing item) data.rows))))
        (setf data (mapcar #'cdr data)))

      (print ranklist)
      ;;determine the friedman statistic 
      (setf f.stat (* (/ (* 12 data.rows) (* data.cols (+ data.cols 1)))
                      (- (loop for item in ranklist
                               summing (expt item 2))
                         (/ (* data.cols (expt (+ data.cols 1) 2)) 4))))

      ;;iman and davenport's "better" f.stat
      (setf f.f (/ (* (- data.rows 1) f.stat)
                   (- (* data.rows (- data.cols 1)) f.stat)))
      f.f
  )))

(defun rank (l &optional (ranks (make-hash-table)) (n 0))
  "Returns a hash of the ranks in a sorted list.  All numbers in a run of
   repeated entries get the average rank of that run."
  (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)
          (cond ((= now 0)
                 (setf repeats 1 sum 0 now 0 n 0)))
	  (setf (gethash now ranks) (/ sum repeats))
	  (rank l ranks n)))))