/* 

  Bugs book, example 2.5.1 in Picat.
  "The how many trick"

  This is a changepoint problem.

  Cf my Gamble model gamble_bugs_book_2_5_1.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 : number
  Probabilities (truncated):
  47284.736856080489815: 0.0005000000000000
  46595.351136583449261: 0.0005000000000000
  46441.004669077148719: 0.0005000000000000
  45338.305399431577825: 0.0005000000000000
  45273.217524456194951: 0.0005000000000000
  44620.553284690548026: 0.0005000000000000
  44361.440293705833028: 0.0005000000000000
  43676.795038188676699: 0.0005000000000000
  43294.890294664161047: 0.0005000000000000
  43053.555481928997324: 0.0005000000000000
  .........
  23387.399921357173298: 0.0005000000000000
  23371.699983153106587: 0.0005000000000000
  23276.480927529071778: 0.0005000000000000
  22819.370284743032244: 0.0005000000000000
  22667.459072135425231: 0.0005000000000000
  22369.607613393789507: 0.0005000000000000
  22310.937586571453721: 0.0005000000000000
  21883.351418376965739: 0.0005000000000000
  21102.857352604180051: 0.0005000000000000
  19852.395885078745778: 0.0005000000000000
  mean = 32135.8
  HPD intervals:
  HPD interval (0.84): 27097.00170842656371..37205.45559959030652

  var : ix
  Probabilities:
  12: 0.2375000000000000
  13: 0.1980000000000000
  11: 0.1805000000000000
  14: 0.1190000000000000
  10: 0.1135000000000000
  9: 0.0565000000000000
  15: 0.0505000000000000
  16: 0.0210000000000000
  8: 0.0150000000000000
  17: 0.0030000000000000
  18: 0.0020000000000000
  7: 0.0020000000000000
  6: 0.0010000000000000
  19: 0.0005000000000000
  mean = 12.049
  HPD intervals:
  HPD interval (0.84): 10.00000000000000..14.00000000000000

  var : x
  Probabilities:
  [1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0]: 0.2375000000000000
  [1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0]: 0.1980000000000000
  [1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0]: 0.1805000000000000
  [1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0]: 0.1190000000000000
  [1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0]: 0.1135000000000000
  [1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0]: 0.0565000000000000
  [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0]: 0.0505000000000000
  [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0]: 0.0210000000000000
  [1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0]: 0.0150000000000000
  [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0]: 0.0030000000000000
  [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0]: 0.0020000000000000
  [1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0]: 0.0020000000000000
  [1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0]: 0.0010000000000000
  [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0]: 0.0005000000000000

*/
go ?=>
  reset_store,
  run_model(2_000,$model,[show_probs_trunc,truncate_size=10,mean,show_hpd_intervals,hpd_intervals=[0.84]]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
cum(I,Y) = Res =>
  if I == 1 then
    Res = Y[1] + 1
  else
    Res = cum(I-1,Y) + Y[I] + (I-1)
  end.
  
model() =>
  N = 20,
  Y = gamma_dist_n(3.0,1/0.04, N),
  Cum = [cum(I,Y) : I in 1..N],
  CumStep = [cond(C > 1000.0, C*(I-1),0) : I in 1..N, C = Cum[I] ],
  
  % A change point
  X = [cond(Cum[I] < 1001.0,1,0) : I in 1..N],
  Number = CumStep.max,

  Ix = [ cond(X[I]==1, I,0) : I in 1..N].max,
 
  add("number",Number),
  add("ix",Ix),
  add("x",X).