(let [a-prior 5
      b-prior 5
      z 1
      N 10
      a-post (+ a-prior z)
      b-post (+ b-prior (- N z))]

(- (beta-cdf a-post b-post 0.6)
   (beta-cdf a-post b-post 0.4)))

(def result (atom
(with-default-bayadera
  (let [a 5 b 5
        z 1 N 10]
    (with-release
      [prior-dist (beta a b)
       prior-sampler (sampler prior-dist)
       prior-sample (dataset (sample! prior-sampler))
       prior-pdf (pdf prior-dist prior-sample)
       post (posterior (posterior-model (:binomial likelihoods)
                                        (:beta distributions)))
       post-dist (post (sv (op (binomial-lik-params N z) (beta-params a b))))
       post-sampler (time (doto (sampler post-dist) (mix!)))
       post-sample (dataset (sample! post-sampler))
       post-pdf (scal! (/ 1.0 (evidence post-dist prior-sample))
                       (pdf post-dist post-sample))]

      {:prior {:sample (native (row (p/data prior-sample) 0))
               :pdf (native prior-pdf)}
       :posterior {:sample (native (row (p/data post-sample) 0))
                   :pdf (native post-pdf)}})))))

(defn read-data [in-file]
  (loop [c 0 data (drop 1 (csv/read-csv in-file)) hw (transient [])]
    (if data
      (let [[_ h w] (first data)]
        (recur (inc c) (next data)
               (-> hw (conj! (double (read-string h)))
                   (conj! (double (read-string w))))))
      (op [c] (persistent! hw)))))

(def params-300 (sv (read-data (slurp (io/file "ht-wt-data-300.csv")))))

(def rlr-source (slurp (io/file "robust-linear-regression.h")))

(def rlr-prior
  (cl-distribution-model [(:gaussian source-library)
                          (:uniform source-library)
                          (:exponential source-library)
                          (:t source-library)
                          rlr-source]
                         :name "rlr" :mcmc-logpdf "rlr_mcmc_logpdf"
                         :params-size 7 :dimension 4))

(defn rlr-likelihood [n]
  (cl-likelihood-model rlr-source :name "rlr" :params-size n))

REAL rlr_loglik(__constant const REAL* params, REAL* x) {
    const REAL nu = x[0];
    const REAL b0 = x[1];
    const REAL b1 = x[2];
    const REAL sigma = x[3];
    const uint n = (uint)params[0];
    const bool valid = (0.0f < nu) && (0.0f < sigma);
    if (valid) {
        const REAL scale = t_log_scale(nu, sigma);
        REAL res = 0.0;
        for (uint i = 0; i < n; i = i+2) {
            res += t_log_unscaled(nu,
                                  b0 + b1 * params[i+1],
                                  sigma, params[i+2])
                + scale;
        }
        return res;
    }
    return NAN;
}

REAL rlr_logpdf(__constant const REAL* params, REAL* x) {
    return exponential_log(params[0], x[0] - 1)
        + gaussian_log(params[1], params[2], x[1])
        + gaussian_log(params[3], params[4], x[2])
        + uniform_log(params[5], params[6], x[3]);

}

(def result (atom
(with-default-bayadera
  (with-release
    [prior (distribution rlr-prior)
     prior-dist (prior (sv 10 -100 100 5 10 0.001 1000))

     post (posterior "rlr_300" (rlr-likelihood (dim params-300)) prior-dist)
     post-dist (post params-300)
     post-sampler (sampler post-dist {:limits (sge 2 4 [1 10 -400 100 0 20 0.01 100])})]

    (mix! post-sampler {:step 384})

    (histogram! post-sampler 1000)))))

(defn var-plot [] (plot2d (qa/current-applet) {:width 480 :height 240}))

(defn setup []
  (q/background 0)
  (q/image (show (render-histogram (var-plot) @result 1)) 0 0)
  (q/image (show (render-histogram (var-plot) @result 2)) 480 0)
  (q/image (show (render-histogram (var-plot) @result 3)) 0 260)
  (q/image (show (render-histogram (var-plot) @result 0)) 480 260)

  (q/save "robust-linear-regression.png"))

#_(q/defsketch diagrams :renderer :p2d :size :fullscreen
    :setup setup :middleware [pause-on-error])