/* 

  Papers under 9 hours in Picat.

  From Pascal Bercker: Five-Minute Probability Problem — Netica vs. Python
  https://medium.com/@pbercker/five-minute-probability-problem-with-netica-d05299e454d5
  """
  I saw this simple little problem here 
  [https://www.youtube.com/watch?v=slbZ-SLpIgg ]
  using a Monte Carlo Simulation with Python. We will redo this with a 
  simple Bayesian Belief Network (Netica) [https://www.norsys.com/netica.html].

  Someone is tasked with writing two reports under a deadline. One report 
  could be anywhere from 1 to 5 hours, and the other anywhere from 2 to 6 
  hours. What is the probability you can get both done under 9 hours? 
  We want a probability distribution over the estimated time to completion. 
  """


  Cf my Gamble model gamble_papers_under_9_hours.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 : p
  Probabilities:
  true: 0.8830000000000000
  false: 0.1170000000000000
  mean = 0.883

  var : total
  Probabilities:
  7: 0.2061000000000000
  8: 0.1623000000000000
  6: 0.1612000000000000
  9: 0.1220000000000000
  5: 0.1125000000000000
  4: 0.0790000000000000
  10: 0.0786000000000000
  3: 0.0399000000000000
  11: 0.0384000000000000
  mean = 7.0129

  var : paper1
  Probabilities:
  4: 0.2038000000000000
  1: 0.2001000000000000
  2: 0.1992000000000000
  3: 0.1988000000000000
  5: 0.1981000000000000
  mean = 3.0006

  var : paper2
  Probabilities:
  6: 0.2037000000000000
  4: 0.2028000000000000
  2: 0.1990000000000000
  5: 0.1987000000000000
  3: 0.1958000000000000
  mean = 4.0123

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  Paper1 = 1 + random_integer(5),
  Paper2 = 2 + random_integer(5),  
  Total = Paper1 + Paper2,

  P = check(Total <= 9),

  add("p",P),
  add("total",Total),
  add("paper1",Paper1),
  add("paper2",Paper2).


/*
  Uniform distribution

  var : p
  Probabilities:
  1: 0.8774999999999999
  0: 0.1225000000000000
  mean = 0.8775
  HPD intervals:
  HPD interval (0.84): 1.00000000000000..1.00000000000000

  var : total
  Probabilities (truncated):
  10.89590526367347: 0.0001000000000000
  10.891206426495316: 0.0001000000000000
  10.871546188309578: 0.0001000000000000
  10.864835081559063: 0.0001000000000000
  .........
  3.099633069755339: 0.0001000000000000
  3.096429059326849: 0.0001000000000000
  3.093263773291029: 0.0001000000000000
  3.071430258486155: 0.0001000000000000
  mean = 6.98381
  HPD intervals:
  HPD interval (0.84): 4.48812390374398..9.29459089706400

  var : paper1
  Probabilities (truncated):
  4.999537274241233: 0.0001000000000000
  4.999150104820332: 0.0001000000000000
  4.998777760192183: 0.0001000000000000
  4.998306766151593: 0.0001000000000000
  .........
  1.001745915041187: 0.0001000000000000
  1.001607691869888: 0.0001000000000000
  1.00068890862199: 0.0001000000000000
  1.000133164227071: 0.0001000000000000
  mean = 2.99225
  HPD intervals:
  HPD interval (0.84): 1.03981102446085..4.38655558944985

  var : paper2
  Probabilities (truncated):
  5.999771052971376: 0.0001000000000000
  5.999410826712572: 0.0001000000000000
  5.999295657500296: 0.0001000000000000
  5.99927509762313: 0.0001000000000000
  .........
  2.001710511744819: 0.0001000000000000
  2.000289166346327: 0.0001000000000000
  2.000282784924043: 0.0001000000000000
  2.000050166621828: 0.0001000000000000
  mean = 3.99156
  HPD intervals:
  HPD interval (0.84): 2.55355379011182..5.90217264178310


*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[show_probs_trunc,mean,show_hpd_intervals, hpd_intervals=[0.84]]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  Paper1 = uniform(1,5),
  Paper2 = uniform(2,6),
  Total = Paper1 + Paper2,

  P = check1(Total <= 9),

  add("p",P),
  add("total",Total),
  add("paper1",Paper1),
  add("paper2",Paper2).