/* 

  Galaxy in Picat.

  From BLOG example/galaxy.blog
  (run the BLOG model with -s blog.sample.MHSampler to get some result)

  Cf my Gamble model gamble_galaxy.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 : g0
  Probabilities (truncated):
  49910.028818952865549: 0.0010000000000000
  49878.552174139091221: 0.0010000000000000
  49859.54298165604996: 0.0010000000000000
  49859.415304315938556: 0.0010000000000000
  .........
  5217.632071682081005: 0.0010000000000000
  5208.685060594550123: 0.0010000000000000
  5188.906597992827301: 0.0010000000000000
  5029.069655588395108: 0.0010000000000000
  mean = 27254.6
  HPD intervals:
  HPD interval (0.84): 5188.90659799282730..41989.03480637307803

  var : g5
  Probabilities (truncated):
  49981.893517534197599: 0.0010000000000000
  49972.033423824250349: 0.0010000000000000
  49945.601804622267082: 0.0010000000000000
  49912.624820094846655: 0.0010000000000000
  .........
  5356.110457962430701: 0.0010000000000000
  5221.120934105068955: 0.0010000000000000
  5149.671328323740454: 0.0010000000000000
  5083.753548135866367: 0.0010000000000000
  mean = 27962.9
  HPD intervals:
  HPD interval (0.84): 11182.59433479169093..48016.35358855890081

  var : g9
  Probabilities (truncated):
  49971.956980867289531: 0.0010000000000000
  49965.80317614870728: 0.0010000000000000
  49922.16214766826306: 0.0010000000000000
  49917.380229997157585: 0.0010000000000000
  .........
  5130.676783216501462: 0.0010000000000000
  5112.94518137022169: 0.0010000000000000
  5075.148225331282447: 0.0010000000000000
  5057.621749610463667: 0.0010000000000000
  mean = 27307.2
  HPD intervals:
  HPD interval (0.84): 7338.30782693731999..45464.40999278072559

  var : g14
  Probabilities (truncated):
  49980.215816749361693: 0.0010000000000000
  49969.334094770878437: 0.0010000000000000
  49936.263211973135185: 0.0010000000000000
  49913.208347239160503: 0.0010000000000000
  .........
  5125.835081620996789: 0.0010000000000000
  5069.266455280252558: 0.0010000000000000
  5048.386044822813346: 0.0010000000000000
  5002.956548707074035: 0.0010000000000000
  mean = 27518.4
  HPD intervals:
  HPD interval (0.84): 12553.39877566015275..49936.26321197313518

  var : g15
  Probabilities (truncated):
  49976.983503427814867: 0.0010000000000000
  49945.655490712102619: 0.0010000000000000
  49901.452371804720315: 0.0010000000000000
  49820.125761823786888: 0.0010000000000000
  .........
  5116.479033658503795: 0.0010000000000000
  5067.682338444367815: 0.0010000000000000
  5059.669762877593712: 0.0010000000000000
  5021.426093308919008: 0.0010000000000000
  mean = 26630.8
  HPD intervals:
  HPD interval (0.84): 5021.42609330891901..42124.92333125552250

  var : g31
  Probabilities (truncated):
  49992.235968817134562: 0.0010000000000000
  49897.040242746952572: 0.0010000000000000
  49889.568886668217601: 0.0010000000000000
  49838.292300625842472: 0.0010000000000000
  .........
  5171.342571345782744: 0.0010000000000000
  5142.736015907831643: 0.0010000000000000
  5095.738433346030433: 0.0010000000000000
  5025.676263042574647: 0.0010000000000000
  mean = 27689.7
  HPD intervals:
  HPD interval (0.84): 10474.49243277054848..47863.57156832868350

  var : g40
  Probabilities (truncated):
  49729.938118127145572: 0.0010000000000000
  49728.64028519421845: 0.0010000000000000
  49692.491902826579462: 0.0010000000000000
  49689.477225620983518: 0.0010000000000000
  .........
  5250.544326962225568: 0.0010000000000000
  5231.108369879008933: 0.0010000000000000
  5082.996767052913128: 0.0010000000000000
  5019.572812607312699: 0.0010000000000000
  mean = 27523.5
  HPD intervals:
  HPD interval (0.84): 7904.09796307985562..44657.52319184016233

  var : g70
  Probabilities (truncated):
  49952.242849372443743: 0.0010000000000000
  49872.700390346675704: 0.0010000000000000
  49847.56520010883105: 0.0010000000000000
  49822.188443421473494: 0.0010000000000000
  .........
  5127.491597145559354: 0.0010000000000000
  5046.615786872159333: 0.0010000000000000
  5046.354732497760779: 0.0010000000000000
  5028.394639970918433: 0.0010000000000000
  mean = 27592.9
  HPD intervals:
  HPD interval (0.84): 7175.41526173027887..44231.52983152844536

  var : g80
  Probabilities (truncated):
  49998.698814305804262: 0.0010000000000000
  49861.559374193457188: 0.0010000000000000
  49852.674996877402009: 0.0010000000000000
  49818.569892001600238: 0.0010000000000000
  .........
  5115.412604117492265: 0.0010000000000000
  5097.028827339890086: 0.0010000000000000
  5065.814473231236661: 0.0010000000000000
  5063.57407433147273: 0.0010000000000000
  mean = 27244.3
  HPD intervals:
  HPD interval (0.84): 11241.03943875107871..48647.71564660953300

  var : num cluster
  Probabilities (truncated):
  7: 0.1710000000000000
  6: 0.1430000000000000
  8: 0.1310000000000000
  9: 0.1200000000000000
  .........
  2: 0.0080000000000000
  14: 0.0060000000000000
  15: 0.0040000000000000
  18: 0.0010000000000000
  mean = 7.404
  HPD intervals:
  HPD interval (0.84): 3.00000000000000..10.00000000000000

  var : orig cluster 0
  Probabilities (truncated):
  2: 0.1710000000000000
  3: 0.1470000000000000
  1: 0.1470000000000000
  4: 0.1370000000000000
  .........
  11: 0.0120000000000000
  12: 0.0030000000000000
  14: 0.0020000000000000
  17: 0.0010000000000000
  mean = 4.19
  HPD intervals:
  HPD interval (0.84): 1.00000000000000..7.00000000000000

*/
go ?=>
  Data = [9172,9350,9483,9558,9775,10227,10406,16084,16170,18419,18552,18600,18927,19052,
          19070,19330,19343,19349,19440,19473,19529,19541,19547,19663,19846,19856,19863,
          19914,19918,19973,19989,20166,20175,20179,20196,20215,20221,20415,20629,20795,
          20821,20846,20875,20986,21137,21492,21701,21814,21921,21960,22185,22209,22242,
          22249,22314,22374,22495,22746,22747,22888,22914,23206,23241,23263,23484,23538,
          23542,23666,23706,23711,24129,24285,24289,24366,24717,24990,25633,26960,26995,
          32065,32789,34279],
  
  reset_store,
  run_model(30_000,$model(Data),[show_probs_trunc,mean,
                           % show_percentiles,show_histogram,
                           show_hpd_intervals,hpd_intervals=[0.84],
                           min_accepted_samples=1000,show_accepted_samples=false
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  

orig_cluster(G,NumCluster) = random_integer1(NumCluster).
clus_velocity(C) = uniform(5000,50000).

model(Data) =>
  NumCluster = 1 + poisson_dist(10.0),
  % ClusVelocity = uniform_n(5000,50000,NumCluster),

  % OrigCluster = random_integer1_n(NumCluster,NumCluster),
  Y = [ clus_velocity(orig_cluster(G,NumCluster)): G in 1..NumCluster],
  
  observe_abc(Data,Y,1/2),
  
  if observed_ok then
    foreach(G in [0,5,9,14,15,31,40,70,80])
      add("g"++G.to_string,clus_velocity(orig_cluster(G,NumCluster))),
    end,
    add("num cluster",NumCluster),
    add("orig cluster 0",orig_cluster(0,NumCluster)),
  end.