(load "evett-5.lisp")
(load "asterisk.lisp")

(defparameter *x-freqs*
  '((1.46833  4)
    (1.476159 0)
    (1.4839879 18)
    (1.4918169 0)
    (1.4996458 0)
    (1.5074748 114)
    (1.5153037 0)
    (1.5231327 4106)
    (1.5309616 27)
    (1.5387906 0)))

(defparameter *y-freqs*
  '((1.46833  4)
    (1.476159  0)
    (1.4839879  18)
    (1.4918169  0)
    (1.4996458  0)
    (1.5074748  114)
    (1.5153037  0)
    (1.5231327  4106)
    (1.5309616  27)
    (1.5387906  0)))

(defparameter *sd-freqs*
  '((1  4024)
    (81/10  0)
    (76/5  235)
    (223/10  0)
    (147/5  7)
    (73/2  0)
    (218/5  0)
    (507/10  2)
    (289/5  0)
    (649/10  1)))

(defmodel sh
	     ((:scale x 1.4 1.6 .001)
	      (:scale y 1.4 1.6 .001)
	      (:scale sd .00001 .00025 .00001)
	      (:scale mu 1.5181932 15.181932 1.5181932)
	      (:scale tau 0.0036737733 0.036737733 0.0036737733))


    ((format t "Wut~%")
     (let* ((l (/ sd tau))
         (ll (* l l))
         (u (/ (- x y) (* sd (sqrt 2))))
         (z (* .5 (+ x y)))
         (v (/ (- z mu) (* tau (expt (+ 1 (* .5 ll)) .5))))
         (factor (* (/ (+ 1 ll) (* l (expt (+ 2 ll) .5)))
                    (exp (* (/ -1 (* 2 (+ 1 ll)))
                            (- (* u u) (* v v)))))))
    factor))
                          

    ((format t "Start!~%")
     (let ((s (demo-seheult)) str (xf 0))
       (labels ((between (left right val) (and (<= left val) (<= val right))))

         (dotimes (i (length *x-freqs*) xf)
           (when (<= (value x) (car (nth i *x-freqs*)))
               (unless (< (1- i) 0)
                 (configure x :label (format nil "x (~A)" (second (nth (1- i) *x-freqs*))))
                 (return))))

         (format t "Frequency of x=~A is ~A~%" (value x) xf)
         
         (dolist (NANER (gethash 'x s))
           (when (between (first NANER) (second NANER) (value x))
             (configure x :troughcolor (intern (format nil "gray~A" (round (* 100 (first (last NANER)))))))
             (format t "Value: ~A~%" (round (* 100 (first (last NANER)))))))

         (dotimes (i (length *y-freqs*) xf)
           (when (<= (value y) (car (nth i *y-freqs*)))
               (unless (< (1- i) 0)
                 (configure y :label (format nil "y (~A)" (second (nth (1- i) *y-freqs*))))
                 (return))))
         
         (dolist (NANER (gethash 'y s))
           (if (between (first NANER) (second NANER) (value y))
               (configure y :troughcolor (intern (format nil "gray~A" (round (* 100 (first (last NANER)))))))))

         (dotimes (i (length *sd-freqs*) xf)
           (when (<= (value sd) (car (nth i *sd-freqs*)))
             (unless (< (1- i) 0)
               (configure sd :label (format nil "x (~A)" (second (nth (1- i) *sd-freqs*))))
               (return))))
         
         (dolist (NANER (gethash 'sd s))
           (if (between (first NANER) (second NANER) (value sd))
               (configure sd :troughcolor (intern (format nil "gray~A" (round (* 100 (first (last NANER)))))))))
         
         (dolist (NANER (gethash 'mu s))
           (if (between (first NANER) (second NANER) (value mu))
               (configure mu :troughcolor (intern (format nil "gray~A" (round (* 100 (first (last NANER)))))))))
         
         (dolist (NANER (gethash 'tau s))
           (if (between (first NANER) (second NANER) (value tau))
               (configure tau :troughcolor (intern (format nil "gray~A" (round (* 100 (first (last NANER)))))))))

         (format t "END~%"))))
  
  (format t "Woot~%")
  (new-random-seheult 1.4 1.6 1 1 mu tau)
  (let (tempval)
    (handler-case
        (progn
          (setf tempval (sh-model x y sd mu tau)))
      
      (floating-point-overflow ()
        (setf tempval 10)))
    (if (> tempval 10)
        (setf tempval 10))
    (with-open-file (out "randomdata.dat"
                         :direction :output
                         :if-exists :append)
      (format out "~%~%~a ~a ~a" x y tempval))
    tempval)
  (format-gnuplot "set title 'Seheult Data'")
  (format-gnuplot "unset key")
  (set-labels "x-ri" "y-ri" "LR")
  (new-generate-thumbs 30 40 40 30)
  (format-gnuplot "splot 'randomdata.dat' index 0, 'randomdata.dat' index 1, 'seheultless10.dat' index 0, 'seheultless10.dat' index 1, 'seheultless10.dat' index 2, 'seheultless10.dat' index 3, 'seheultless10.dat' index 4, 'seheultless10.dat' index 5, 'seheultless10.dat' index 6, 'seheultless10.dat' index 7, 'seheultless10.dat' index 8, 'seheultless10.dat' index 9, 'seheultless10.dat' index 10, 'seheultless10.dat' index 11, 'seheultless10.dat' index 12, 'seheultless10.dat' index 13, 'seheultless10.dat' index 14, 'seheultless10.dat' index 15, 'seheultless10.dat' index 16, 'seheultless10.dat' index 17, 'seheultless10.dat' index 18, 'seheultless10.dat' index 19, 'seheultless10.dat' index 20, 'seheultless10.dat' index 21, 'seheultless10.dat' index 22, 'seheultless10.dat' index 23, 'seheultless10.dat' index 24, 'seheultless10.dat' index 25, 'seheultless10.dat' index 26, 'seheultless10.dat' index 27, 'seheultless10.dat' index 28, 'seheultless10.dat' index 29"))
  

(defmacro new-generate-thumbs (viewx viewy dx dy)
  `(progn
     (format-gnuplot "set terminal gif")
     (format-gnuplot "set size .6,.6")
     (dotimes (i 3)
       (format-gnuplot (format nil "set view ~A,~A" (+ ,viewx (* (1+ i) ,dx)) (+ ,viewy (* (1+ i) ,dy))))
       (format-gnuplot (format nil "set output \"~A.gif\"" i))
       (format-gnuplot "splot 'randomdata.dat' index 0, 'randomdata.dat' index 1, 'seheultless10.dat' index 0, 'seheultless10.dat' index 1, 'seheultless10.dat' index 2, 'seheultless10.dat' index 3, 'seheultless10.dat' index 4, 'seheultless10.dat' index 5, 'seheultless10.dat' index 6, 'seheultless10.dat' index 7, 'seheultless10.dat' index 8, 'seheultless10.dat' index 9, 'seheultless10.dat' index 10, 'seheultless10.dat' index 11, 'seheultless10.dat' index 12, 'seheultless10.dat' index 13, 'seheultless10.dat' index 14, 'seheultless10.dat' index 15, 'seheultless10.dat' index 16, 'seheultless10.dat' index 17, 'seheultless10.dat' index 18, 'seheultless10.dat' index 19, 'seheultless10.dat' index 20, 'seheultless10.dat' index 21, 'seheultless10.dat' index 22, 'seheultless10.dat' index 23, 'seheultless10.dat' index 24, 'seheultless10.dat' index 25, 'seheultless10.dat' index 26, 'seheultless10.dat' index 27, 'seheultless10.dat' index 28, 'seheultless10.dat' index 29"))
     (format-gnuplot "set terminal wxt")
     (format-gnuplot "set size 1,1")
     (format-gnuplot "set output")))

(defun new-random-seheult (min max min.sd max.sd
                           &optional
                           (mu 1.5181932) (tau 0.0036737733) (timespersd 5000) (filename "randomdata.dat"))
  (let ((x 0)
        (y 0)
        (tempval 0))
    (with-open-file (out filename
                         :direction :output
                         :if-exists :supersede)
      (format out "### splot '~a' index 0" filename)
      (dotimes (i (- max.sd min.sd))
        (format out ", '~a' index ~a" filename (+ i 1)))
      (dotimes (i (+ (- max.sd min.sd) 1))
                                        ;(format out "~%~%###SD: ~a~%" (+ i 1))
        (format out "~%~%###SD: ~a~%" (+ i min.sd))
        (dotimes (j timespersd)
          (setf x (+ (random (- max min)) min))
          (setf y (+ (random (- max min)) min))
          (handler-case
              (progn
                (setf tempval
                      (seheult x y (* (+ i 1)
                                      (expt 10 -5)) mu tau)))
            (floating-point-overflow ()
              (setf tempval 10)))
          (if (> tempval 10)
              (setf tempval 10))
          (format out "~a ~a ~a~%" x y tempval))))
    'done))

(defun seheult (x y sd &optional (mu 1.5181932)
                (tau 0.0036737733))
  (let* ((l (/ sd tau))
         (ll (* l l))
         (u (/ (- x y) (* sd (sqrt 2))))
         (z (* .5 (+ x y)))
         (v (/ (- z mu) (* tau (expt (+ 1 (* .5 ll)) .5))))
         (factor (* (/ (+ 1 ll) (* l (expt (+ 2 ll) .5)))
                    (exp (* (/ -1 (* 2 (+ 1 ll)))
                            (- (* u u) (* v v)))))))
    factor))

(defun handle-seheult (x y sd &optional  (mu 1.5181932)
                       (tau 0.0036737733) (maxallowed 10) (filename "sliderdata.dat"))
  (let (tempval)
    (handler-case
        (progn
          (setf tempval (seheult x y sd mu tau)))
      
      (floating-point-overflow ()
        (setf tempval maxallowed)))
    (if (> tempval maxallowed)
        (setf tempval maxallowed))
    (with-open-file (out filename
                         :direction :output
                         :if-exists :supersede)
      (format out "~a ~a ~a" x y tempval))
    tempval))