/* 

  Clinical trial (R2) in Picat.

  From R2
    with the datafiles 
  - clinical_trial_control.csv  (size 1000)
  - clinical_trial_treatment.csv (size 1000)

  length :1000 mean(control): 0.513 mean(treatment): 0.51

  Output from the R2 model ClinicalTrial.cs (for isEffective)
  """
  Mean: 0.113
  Variance: 0.100331
  Number of accepted samples = 318
  0:00:14:048.362.900
  """

  Cf my Gamble model gamble_clinical_trial_r2.rkt

  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

import ppl_distributions, ppl_utils.
import util.
% import ordset.

main => go.

/*

  var : is effective
  mean = [false = 0.896,true = 0.104]

  var : prob control
  mean = 0.509796

  var : prob treatment
  mean = 0.519695

  var : prob all
  mean = 0.505956


*/
go ?=>
  Control = [0,0,1,1,1,0,1,1,0,0,0,0,0,0,1,0,1,0,0,1,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,0,0,1,1,0,0,1,0,0,1,0,0,1,1,1,0,0,1,0,1,0,1,1,0,1,0,0,1,1,1,0,0,0,1,0,1,1,0,0,1,1,1,0,0,0,1,1,0,0,0,1,1,0,0,0,1,0,1,1,0,0,0,1,0,1,1,1,0,0,1,1,0,0,1,0,1,1,1,1,1,1,1,1,0,0,1,0,0,0,1,0,1,1,0,1,0,0,1,0,1,1,1,1,0,1,1,0,1,0,0,0,1,0,0,1,0,1,0,0,0,1,1,1,1,1,1,1,0,1,1,1,1,0,0,1,1,0,0,1,0,0,0,0,0,1,1,1,0,1,1,1,0,0,0,1,0,1,0,1,1,0,0,1,1,1,0,0,0,0,1,0,0,0,0,1,1,1,0,0,1,0,0,1,0,1,1,0,1,0,1,0,1,1,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,1,0,0,1,0,1,0,0,0,1,1,0,1,0,1,0,0,0,1,0,0,1,0,0,0,1,1,1,1,1,1,0,1,0,0,0,0,0,1,1,1,1,1,1,1,0,0,1,1,0,1,0,1,1,1,1,1,0,1,0,0,1,0,1,1,1,0,1,0,0,1,0,1,1,0,0,1,1,1,1,1,1,1,1,0,1,1,0,0,0,1,0,1,0,1,0,1,0,1,0,1,1,1,1,0,1,0,1,0,1,0,1,0,1,1,0,1,1,1,1,0,0,0,1,1,0,0,0,0,0,1,0,0,1,1,1,1,0,0,0,0,1,0,1,0,1,0,0,0,0,0,1,1,1,1,1,0,0,1,0,0,0,1,0,0,1,1,1,1,1,1,1,1,0,0,0,1,1,1,1,0,1,0,1,1,0,1,0,1,1,1,1,0,0,1,0,1,0,1,0,0,1,0,1,1,1,0,0,0,0,0,0,1,1,1,1,1,0,0,0,1,0,1,1,0,0,0,0,0,0,0,1,0,1,0,1,1,0,1,0,1,1,1,0,1,0,1,1,0,0,0,1,1,0,0,0,1,0,0,1,0,1,1,0,0,0,0,1,0,0,1,1,0,0,0,0,1,1,1,1,0,1,1,1,1,1,0,1,0,0,1,0,0,0,1,0,1,0,0,1,1,1,0,0,0,1,1,1,0,1,1,0,0,1,0,1,1,1,1,1,0,1,1,0,0,1,1,1,0,1,1,0,0,0,1,0,1,1,1,1,1,0,0,1,1,1,1,1,1,0,1,1,0,1,1,0,0,0,0,1,0,1,0,1,0,0,1,0,1,0,0,0,0,0,1,1,0,1,0,0,1,0,0,0,0,0,1,1,1,1,1,0,1,1,0,0,0,1,1,1,0,0,1,0,0,1,0,0,0,1,0,1,1,1,1,0,0,0,0,1,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0,0,0,1,1,1,1,1,1,1,0,0,1,1,1,0,1,1,0,1,0,1,1,1,1,1,0,1,0,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0,1,1,1,1,1,1,1,0,0,1,1,0,1,0,1,0,0,0,0,0,1,1,0,0,1,0,1,1,0,1,1,1,0,1,1,0,1,1,1,1,1,0,1,1,1,0,1,0,1,0,0,1,0,1,0,0,0,1,1,1,0,0,1,1,1,0,1,0,1,0,0,1,1,0,0,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,0,1,1,1,1,0,1,0,1,0,1,0,0,1,1,0,1,0,1,1,0,0,0,0,1,1,1,0,0,1,0,1,1,1,0,1,1,0,0,1,0,0,0,0,1,0,1,1,0,0,0,0,0,1,1,1,0,0,0,1,0,1,1,0,1,1,1,1,0,1,0,1,1,0,0,0,1,0,1,1,0,0,1,0,0,1,1,1,0,1,0,1,0,0,0,1,0,0,0,0,0,0,0,0,1,1,1,0,1,1,1,0,1,1,1,0,1,1,0,0,0,1,1,1,0,1,0,0,0,1,0,1,1,0,0,0,1,1,1,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,1,0,1,1,1,0,1,1,1,1,0,1,0,0,1,1,0,0,1,0,1,0,0,1,0,1,1,0,1,0,1],
  Treatment = [0,1,0,0,1,1,1,1,1,1,0,0,1,1,0,1,0,0,1,0,1,0,1,0,0,0,0,1,1,0,0,0,1,1,1,0,1,1,1,1,1,1,1,1,0,1,0,1,0,0,1,0,1,1,1,1,1,1,0,0,1,0,1,0,0,0,1,0,1,0,1,0,0,1,1,0,1,0,0,0,1,0,0,1,1,1,1,0,0,1,1,0,1,1,0,0,1,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,1,0,1,1,0,1,1,0,1,1,0,0,1,0,1,1,1,0,0,0,1,0,1,1,0,0,0,1,1,1,1,0,0,0,1,0,0,1,1,1,0,1,1,0,0,1,1,0,0,1,0,1,0,1,0,1,1,1,0,0,1,1,1,0,0,1,1,0,0,0,0,1,1,0,1,1,1,0,1,1,0,1,0,1,0,1,0,0,1,0,0,1,1,1,1,1,0,1,1,0,1,1,1,0,1,0,0,1,0,0,0,0,1,0,0,1,1,1,1,1,1,1,1,0,1,1,0,1,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,0,1,0,1,0,0,0,1,0,1,0,0,1,0,0,0,0,1,0,0,0,0,1,1,1,1,0,1,1,0,0,1,1,1,0,0,0,0,0,1,0,1,1,1,1,1,0,0,1,1,0,1,0,1,1,1,1,1,0,1,1,1,0,1,1,1,1,0,1,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,0,1,1,1,1,1,0,0,1,0,0,0,0,1,1,0,1,1,1,1,0,1,1,1,0,1,0,0,1,1,1,1,0,0,1,0,0,1,1,0,0,0,0,1,1,0,1,1,0,1,0,0,0,0,0,1,1,1,0,0,0,1,0,1,1,1,1,1,0,1,1,0,1,0,0,0,0,0,0,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,1,1,1,0,0,1,1,0,0,1,1,1,1,0,0,0,0,0,0,1,1,0,0,1,1,1,0,1,0,1,0,0,1,1,1,1,0,0,1,1,0,0,0,1,1,0,1,1,0,1,0,1,1,1,1,0,0,0,0,1,0,1,1,0,0,0,0,0,1,1,0,1,0,0,0,0,0,1,0,0,0,1,1,1,0,1,1,1,1,0,1,1,1,0,0,1,0,1,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,1,1,1,1,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,0,1,0,1,1,0,0,1,0,1,1,1,0,0,1,0,1,1,1,0,1,0,1,1,0,1,1,1,0,0,0,1,0,1,0,0,1,0,0,0,1,0,1,0,1,0,0,1,1,1,1,0,0,1,1,0,0,1,0,1,1,1,0,0,1,0,0,0,0,0,0,0,1,1,0,1,0,0,0,1,1,0,0,0,0,1,0,1,0,1,1,1,0,0,1,0,1,0,0,0,0,1,1,0,1,0,1,1,1,0,0,1,0,0,1,1,1,1,0,1,0,1,0,0,1,0,0,0,1,1,1,1,0,0,1,1,1,1,1,0,0,0,1,1,1,1,1,1,0,1,1,1,0,1,0,1,0,1,1,0,1,0,1,0,0,1,1,0,0,0,0,1,1,0,1,1,0,1,0,1,1,0,1,0,1,1,0,0,1,0,1,1,0,1,0,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0,0,1,1,1,0,0,0,0,0,0,0,1,0,1,0,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,1,0,1,1,1,1,1,0,1,1,0,1,1,1,1,0,1,0,0,0,0,0,1,0,1,1,0,1,1,0,1,1,1,1,1,1,0,0,0,0,0,1,1,0,0,1,1,1,1,0,1,0,1,1,1,1,0,0,1,1,0,0,1,1,0,1,1,1,0,1,0,0,0,0,1,0,0,1,1,1,1,1,0,0,0,1,0,0,0,1,1,1,1,0,0,0,0,0,1,1,0,0,0,1,0,1,1,0,1,0,1,1,0,1,1,1,0,1,0,1,0,1,0,0,0,0,0,1,0,1,0,1,1,0,1,1,0,0,1,0,1,1,0,1,0,1,0,1,1,1,0,0,0,0,0,0,1,0,0,1,0,1,0,0],
  reset_store,
  println(control=[mean=Control.mean,stdev=Control.stdev]),
  println(treatment=[mean=Treatment.mean,stdev=Treatment.stdev]),  
  run_model(10_000,$model(Control,Treatment),[% show_probs_trunc,
                                              mean,
                           % show_percentiles,show_histogram,
                           % show_hpd_intervals,hpd_intervals=[0.84],
                           min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model(Control,Treatment) =>
  N = Control.len,

  IsEffective   = flip(0.5),
  ProbControl   = beta_dist(1,1),
  ProbTreatment = beta_dist(1,1),
  ProbAll       = beta_dist(1,1),

  C = [ cond(IsEffective==true,bern(ProbControl),bern(ProbAll)) : I in 1..N],
  T = [ cond(IsEffective==true,bern(ProbTreatment),bern(ProbAll)) : I in 1..N],

  % Scale the acceptable error interval
  observe_abc(Control[1..N],C,1/14),
  observe_abc(Treatment[1..N],T,1/14),

  if observed_ok then
    add("is effective",IsEffective),
    add("prob control",ProbControl),
    add("prob treatment",ProbTreatment),
    add("prob all",ProbAll),    
  end.