/* 

  Robot localization in Picat.

  From the SPPL model robot-localization.pynb
  The example just shows the model.

  For this model we observe that x is in 1..3


  Cf my Gamble model gamble_robot_localization.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.

main => go.

/*

  var : param
  Probabilities (truncated):
  19: 0.0386740331491713
  17: 0.0381215469613260
  18: 0.0370165745856354
  21: 0.0353591160220995
  .........
  55: 0.0005524861878453
  54: 0.0005524861878453
  50: 0.0005524861878453
  49: 0.0005524861878453
  mean = 21.8011
  HPD intervals:
  HPD interval (0.94): 0.00000000000000..40.00000000000000

  var : which
  Probabilities (truncated):
  0.949388316808915: 0.0005524861878453
  0.948285119583963: 0.0005524861878453
  0.947755493665: 0.0005524861878453
  0.946736392540269: 0.0005524861878453
  .........
  0.004249067979049: 0.0005524861878453
  0.002236513887642: 0.0005524861878453
  0.002079274040684: 0.0005524861878453
  0.000348746776743: 0.0005524861878453
  mean = 0.472915
  HPD intervals:
  HPD interval (0.94): 0.00444186432494..0.89108719438831

  var : x
  Probabilities (truncated):
  2.999902465250844: 0.0005524861878453
  2.99809129147233: 0.0005524861878453
  2.997736403931988: 0.0005524861878453
  2.996987481275492: 0.0005524861878453
  .........
  1.001699483406017: 0.0005524861878453
  1.001150748111563: 0.0005524861878453
  1.000511236470416: 0.0005524861878453
  1.000146812049666: 0.0005524861878453
  mean = 1.99052
  HPD intervals:
  HPD interval (0.94): 1.01429061254397..2.88310600613847

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,show_hpd_intervals,hpd_intervals=[0.94]]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
  
model() =>
  Param = random_integer(101),
  Which = uniform(0,1),

  X = cases([ [Which < 0.9,normal_dist(Param/10,1)],
              [Which < 0.95,uniform(0,10)],
              [true,10.0]]),

  observe(X >= 1.0,X <= 3.0),
  
 
  if observed_ok then
    add("param",Param),
    add("which",Which),
    add("x",X),    
  end.