;;;;Peter Santiago, Started June 16th 2008
;;;;compiled under:
;;;;slime 1:20070927-3
;;;;SBCL  1:1.0.11.0-1
;;;;Wilcoxon algorithm and information shamelessly stolen from Richard Lowry:
;;;;http://faculty.vassar.edu/lowry/ch12a.html

;;;To use the wilcoxon signed-rank test, the data must have the following properties:
;;; 1. the scale of measurement for XA and XB has the properties of an equal-interval scale
;;; 2. the differences between the paired values of XA and XB have been randomly drawn from the source population;
;;; 3. the source population from which these differences have been drawn can be reasonably supposed to have a normal distribution. 

(defun wilcoxon (datarray)
  "Datarray is assumed to be a nx2 matrix, where column 0 represents measurements of x_a and column 1 represents measurements of x_b"
  (let ((expanded (adjust-array datarray (list (array-dimension datarray 0) 6)))
        (w 0)
        (n 0)
        (temp nil))


    (dotimes (i (array-dimension expanded 0))
      (setf (aref expanded i 2) (- (aref expanded i 0) (aref expanded i 1)))
      (setf (aref expanded i 3) (abs (aref expanded i 2))))

    (setf temp (rank (quicksort-array expanded 3 0) 3))
    (setf expanded (first temp))
    (setf n (second temp))
    ;(format t "n=~a~%" n)
    ;(format t "~a~%" expanded)
    
    (dotimes (i (array-dimension expanded 0))
      (setf w (+ w (aref expanded i 5))))

    ;(setf n (aref expanded (- (array-dimension expanded 0) 1) 5))
    (/ (- w .5) (sqrt (/ (* n (+ n 1) (+ (* 2 n) 1)) 6)))))

(defun rank (datarray column)
  "The specified column of datarray is ranked.The ranking is stored in the subsequent column. The signed ranking is stored in the final column.  Also, the code looks really, really bad, mostly because I bandaided a few bugs. I could maybe make it better, but 'if it ain't broke....'"
  (let ((temprawval 0)
        (rankiterator 1)
        (numsame 0)
        (tempsum 0))
    
    (dotimes (i (array-dimension datarray 0))
      ;(format t "Pass: ~a~%" i)
      (cond ((= (aref datarray i column) 0) ;special case when value is zero
             ;(format t "Initial loop! Setting ~a to 0~%" i)
             (setf (aref datarray i (+ column 1)) 0))

            ((= i (- (array-dimension datarray 0) 1)) ;special case when you're out of values
             (cond ((= (aref datarray i column) temprawval)
                    (setf tempsum (+ tempsum rankiterator))
                    (dotimes (k (+ numsame 1))
                      ;(format t "Final Loop! Setting ~a to ~a~%" (+ (- i numsame) k) (/ tempsum numsame)) 
                      (setf (aref datarray (+ (- i numsame) k)  (+ column 1)) (/ tempsum numsame))
                      (cond ((< (aref datarray (+ (- i numsame) k) (- column 1)) 0)
                             (setf (aref datarray (+ (- i numsame) k)  (+ column 2)) (* (/ tempsum numsame) -1)))
                            ((>= (aref datarray (+ (- i numsame) k) (- column 1)) 0)
                             (setf (aref datarray (+ (- i numsame) k)  (+ column 2)) (/ tempsum numsame))))))
                   
                   ((/= (aref datarray i column) temprawval)
                    (dotimes (k numsame)
                      ;(format t "Inner Loop! Setting ~a to ~a~%" (+ (- i numsame) k) (/ tempsum numsame)) 
                      (setf (aref datarray (+ (- i numsame) k)  (+ column 1)) (/ tempsum numsame))
                      (cond ((< (aref datarray (+ (- i numsame) k) (- column 1)) 0)
                             (setf (aref datarray (+ (- i numsame) k)  (+ column 2)) (* (/ tempsum numsame) -1)))
                            ((>= (aref datarray (+ (- i numsame) k) (- column 1)) 0)
                             (setf (aref datarray (+ (- i numsame) k)  (+ column 2)) (/ tempsum numsame)))))
                    ;(format t "Setting final value ~a to ~a~%" i rankiterator)
                    (setf (aref datarray i (+ column 1)) rankiterator)
                    (cond ((< (aref datarray i (- column 1)) 0)
                           (setf (aref datarray i  (+ column 2)) (* rankiterator -1)))
                          ((>= (aref datarray i (- column 1)) 0)
                           (setf (aref datarray i (+ column 2)) rankiterator))))))                    

            ((/= (aref datarray i column) 0) ;normal case
             (cond ((/= (aref datarray i column) temprawval)
                    (dotimes (k numsame)
                      ;(format t "Inner Loop! Setting ~a to ~a~%" (+ (- i numsame) k) (/ tempsum numsame)) 
                      (setf (aref datarray (+ (- i numsame) k)  (+ column 1)) (/ tempsum numsame))
                      (cond ((< (aref datarray (+ (- i numsame) k) (- column 1)) 0)
                             (setf (aref datarray (+ (- i numsame) k)  (+ column 2)) (* (/ tempsum numsame) -1)))
                            ((>= (aref datarray (+ (- i numsame) k) (- column 1)) 0)
                             (setf (aref datarray (+ (- i numsame) k)  (+ column 2)) (/ tempsum numsame)))))
                    (setf temprawval 0)
                    (setf numsame 0)
                    (setf tempsum 0)))
             ;;if i moved the below stuff somewhere above, it might make the code a lot nicer looking,
             ;;but the code is moderately close to optimal running-time regardless, so someone in the future(possibly myself)
             ;;is welcome to play with it for aesthetics
             (setf temprawval (aref datarray i column))
             (setf numsame (+ numsame 1))
             (setf tempsum (+ tempsum rankiterator))
             (setf rankiterator (+ rankiterator 1)))))
    (list datarray rankiterator)))
      
                        
(defun quicksort-array (data key &optional (direction 1))
  "quicksort-array performs a quicksort on the given key (0 indexed) in each row, such that the matrix is eventually sorted by the given key. Direction is either 1 or 0: 1 is greatest to smallest, 0 is smallest to greatest."
  (let ((pivotvector (make-array (array-dimension data  1) :initial-element nil))
        (lesser (make-array (list (array-dimension data 0) (array-dimension data 1)) :initial-element nil))
        (greater (make-array (list (array-dimension data 0) (array-dimension data 1)) :initial-element nil))
        (rows.lesser 0)
        (rows.greater 0)
        (x (array-dimension data 0))
        (y (array-dimension data 1)))
    (if (<= x 1)
        (return-from quicksort-array data))
    ;;set pivot
    (dotimes (j y)
      (setf (aref pivotvector j) (aref data (- x 1) j)))
    (setf data (adjust-array data (list (- x 1)  y)))
    (setf x (array-dimension data 0))
    ;;iteration!
    (dotimes (i x)
      (cond ((<= (aref data i key) (aref pivotvector key))
             (dotimes (j y)
               (setf (aref lesser rows.lesser j) (aref data i j)))
             (setf rows.lesser (+ rows.lesser 1)))

            ((> (aref data i key) (aref pivotvector key))
             (dotimes (j y)
               (setf (aref greater rows.greater j) (aref data i j)))
             (setf rows.greater (+ rows.greater 1)))))
    ;;adjust matrices
    (setf lesser (adjust-array lesser (list rows.lesser y)))
    (setf greater (adjust-array greater (list rows.greater y)))
    ;;recursion
    (cond ((= direction 1)
           (array-concat (quicksort-array greater key direction) pivotvector (quicksort-array lesser key direction)))
          ((= direction 0)
           (array-concat (quicksort-array lesser key direction) pivotvector (quicksort-array greater key direction))))))


(defun array-concat (a1 vector a3)
  "array-append: returns a new array with all elements of a1 vector a3 in it. Note: the row size must be the same for all arrays"
  (let ((arr (make-array (list (+ (array-dimension a1 0) 1 (array-dimension a3 0)) (array-dimension a1 1)) :initial-element nil))
        (currow -1))
    (dotimes (i (array-dimension a1 0))
      (setf currow (+ currow 1)) 
      (dotimes (j (array-dimension a1 1))
        (setf (aref arr currow j) (aref a1 i j))))
    
    (setf currow (+ currow 1)) 
    (dotimes (j (array-dimension vector 0))
      (setf (aref arr currow j) (aref vector j)))
    
    (dotimes (i (array-dimension a3 0))
      (setf currow (+ currow 1)) 
      (dotimes (j (array-dimension a3 1))
        (setf (aref arr currow j) (aref a3 i j))))
    arr))

(defun save-item (filename item)
  "save-item: allows the variable item to be stored in the file at filename as a lisp object"
 (with-open-file (out filename
                  :direction :output
                  :if-exists :supersede)
    (with-standard-io-syntax
      (print item out)))) ;end save-item


(defun load-item (filename)
  "load-item: allows the file at filename to be parsed and stored"
  (let ((x 0))
  (with-open-file (in filename)
    (with-standard-io-syntax
      (setf x (read in))
      x)))) ;end load-item