/* 

  Icy road in Picat.

  From Hugin's icy_roads.net
  """
  This example shows d-separation in action. The risks of Holmes and Watson crashing are 
  only dependent the state of icy is not known.
 
  If it is known, that the roads are icy, they both have a large risk of crashing. Likewise, 
  if the roads are not icy they both have a small risk of crashing.

  But if the state of icy is not known, and Holmes crashes, the risk of Watson crashing 
  goes up, since the crash of Holmes indicates that the roads may be icy.
  """


  This is a port of my Gamble model gamble_icy_road.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.

main => go.

/*

  obs = no_observation
  var : holmes
  Probabilities:
  true: 0.5925000000000000
  false: 0.4075000000000000
  mean = [true = 0.5925,false = 0.4075]

  var : icy
  Probabilities:
  true: 0.7003000000000000
  false: 0.2997000000000000
  mean = [true = 0.7003,false = 0.2997]

  var : watson
  Probabilities:
  true: 0.5878000000000000
  false: 0.4122000000000000
  mean = [true = 0.5878,false = 0.4122]


  obs = (watson == true)
  var : holmes
  Probabilities:
  true: 0.7680491551459293
  false: 0.2319508448540707
  mean = [true = 0.768049,false = 0.231951]

  var : icy
  Probabilities:
  true: 0.9472606246799795
  false: 0.0527393753200205
  mean = [true = 0.947261,false = 0.0527394]

  var : watson
  Probabilities:
  true: 1.0000000000000000
  mean = [true = 1.0]


  obs = not watson
  var : holmes
  Probabilities:
  false: 0.6561333665090818
  true: 0.3438666334909181
  mean = [false = 0.656133,true = 0.343867]

  var : icy
  Probabilities:
  false: 0.6486688230903209
  true: 0.3513311769096790
  mean = [false = 0.648669,true = 0.351331]

  var : watson
  Probabilities:
  false: 1.0000000000000000
  mean = [false = 1.0]


  obs = (holmes == true)
  var : holmes
  Probabilities:
  true: 1.0000000000000000
  mean = [true = 1.0]

  var : icy
  Probabilities:
  true: 0.9520094963540784
  false: 0.0479905036459217
  mean = [true = 0.952009,false = 0.0479905]

  var : watson
  Probabilities:
  true: 0.7641173478039681
  false: 0.2358826521960319
  mean = [true = 0.764117,false = 0.235883]


  obs = not holmes
  var : holmes
  Probabilities:
  false: 1.0000000000000000
  mean = [false = 1.0]

  var : icy
  Probabilities:
  false: 0.6649188514357054
  true: 0.3350811485642946
  mean = [false = 0.664919,true = 0.335081]

  var : watson
  Probabilities:
  false: 0.6479400749063671
  true: 0.3520599250936330
  mean = [false = 0.64794,true = 0.35206]


  obs = (holmes == true,watson == true)
  var : holmes
  Probabilities:
  true: 1.0000000000000000
  mean = [true = 1.0]

  var : icy
  Probabilities:
  true: 0.9928170594837261
  false: 0.0071829405162738
  mean = [true = 0.992817,false = 0.00718294]

  var : watson
  Probabilities:
  true: 1.0000000000000000
  mean = [true = 1.0]


  obs = (not holmes,watson == true)
  var : holmes
  Probabilities:
  false: 1.0000000000000000
  mean = [false = 1.0]

  var : icy
  Probabilities:
  true: 0.8017179670722978
  false: 0.1982820329277022
  mean = [true = 0.801718,false = 0.198282]

  var : watson
  Probabilities:
  true: 1.0000000000000000
  mean = [true = 1.0]


  obs = (holmes == true,not watson)
  var : holmes
  Probabilities:
  true: 1.0000000000000000
  mean = [true = 1.0]

  var : icy
  Probabilities:
  true: 0.8052223006351447
  false: 0.1947776993648553
  mean = [true = 0.805222,false = 0.194778]

  var : watson
  Probabilities:
  false: 1.0000000000000000
  mean = [false = 1.0]


  obs = (not wolmes,not watson)
  var : holmes
  Probabilities:
  false: 1.0000000000000000
  mean = [false = 1.0]

  var : icy
  Probabilities:
  false: 0.8942632170978627
  true: 0.1057367829021372
  mean = [false = 0.894263,true = 0.105737]

  var : watson
  Probabilities:
  false: 1.0000000000000000
  mean = [false = 1.0]


  obs = (icy == true)
  var : holmes
  Probabilities:
  true: 0.8052354531227771
  false: 0.1947645468772229
  mean = [true = 0.805235,false = 0.194765]

  var : icy
  Probabilities:
  true: 1.0000000000000000
  mean = [true = 1.0]

  var : watson
  Probabilities:
  true: 0.7965571205007824
  false: 0.2034428794992175
  mean = [true = 0.796557,false = 0.203443]


  obs = not icy
  var : holmes
  Probabilities:
  false: 0.8895309313392250
  true: 0.1104690686607750
  mean = [false = 0.889531,true = 0.110469]

  var : icy
  Probabilities:
  false: 1.0000000000000000
  mean = [false = 1.0]

  var : watson
  Probabilities:
  false: 0.9038069340584637
  true: 0.0961930659415364
  mean = [false = 0.903807,true = 0.0961931]



*/
go ?=>
  member(Obs,0..10),
  reset_store,
  println(obs=case(Obs,[[0,no_observation],
                        [1,watson==true],
                        [2,not watson],
                        [3,holmes==true],
                        [4,not holmes],
                        [5,(holmes==true,watson==true)],
                        [6,(not holmes, watson==true)],
                        [7,(holmes==true, not watson)],
                        [8,(not wolmes, not watson)],
                        [9,icy==true],
                        [10,not icy]]
                      )),
  run_model(10_000,$model(Obs),[show_probs_trunc,mean]),
  nl,
  fail,
  nl.
go => true.

model(Obs) =>
  Icy = flip(0.7),
  % This node gives the risk of Watson crashing in his car, given the state of the roads.
  Watson = condt(Icy, flip(0.8), flip(0.1)),
  
  % This node gives the risk of Holmes crashing in his car, given the state of the roads.
  Holmes = condt(Icy, flip(0.8), flip(0.1)),
  if Obs > 0 then
    observe(case(Obs,[[1,Watson==true],
                      [2,not Watson],
                      [3,Holmes==true],
                      [4,not Holmes],
                      [5,(Holmes==true,Watson==true)],
                      [6,(not Holmes, Watson==true)],
                      [7,(Holmes==true, not Watson)],
                      [8,(not Holmes, not Watson)],
                      [9,Icy==true],
                      [10,not Icy]
                     ])),
  end,

  if Obs > 0, observed_ok then
    add("watson",Watson),
    add("holmes",Holmes),
    add("icy",Icy)
  elseif Obs == 0 then
    add("watson",Watson),
    add("holmes",Holmes),
    add("icy",Icy)
 
  end.