(ns code.models.seheult
  (:use (code.utils utils))
  (:use (incanter core stats)))

(def SEHEULT [])

(defn seheult 
  ([x y sd] (seheult x y sd 1.5181932))
  ([x y sd mu] (seheult x y sd mu 0.0036737733))
  ([x y sd mu tau]
     (let [l      (/ sd tau)
	   ll     (* l l)
	   u      (/ (- x y) (* sd (Math/sqrt 2)))
	   z      (* 0.5 (+ x y))
	   v      (/ (- z mu) (* tau (Math/pow (+ 1 (* 0.5 ll)) 0.5)))
	   factor (* (/ (+ 1 ll) (* l (Math/pow (+ 2 ll) 0.5)))
		     (Math/exp (* (/ -1 (* 2 (+ 1 ll)))
				  (- (* u u) (* v v)))))]
       factor)))

(defn random-seheult
  ([] (random-seheult 1 1.4 1.6 1.4 1.6 1 1 1.5181932 15.181932 0.0036737733 0.036737733))
  ([samples] 
     (random-seheult samples 1.4 1.6 1.4 1.6 1 1 1.5181932 15.181932 0.0036737733 0.036737733))
  ([min-x max-x min-y max-y min-sd max-sd min-mu max-mu min-tau max-tau] 
     (random-seheult 1 min-x max-x min-y max-y min-sd max-sd min-mu max-mu min-tau max-tau))
  ([samples min-x max-x min-y max-y min-sd max-sd min-mu max-mu min-tau max-tau]
     (dotimes [stdev1 (+ (- max-sd min-sd) 1)]
       (dotimes [one samples]
	 (let [x      (nr min-x max-x)
	       y      (nr min-y max-y)
	       sd     (* (+ stdev1 1) (Math/pow 10 -5))
	       mu     (nr min-mu max-mu)
	       tau    (nr min-tau max-tau)
	       lr     (seheult x y sd mu tau)
	       new-lr (if (>= lr 10) 1 0)]
	   (def SEHEULT (conj SEHEULT [x y sd mu tau new-lr])))))
	   
       (let [result SEHEULT]
	 (def SEHEULT [])
	       result)))