/* 

  Election in Picat.

  This is a port of the SPPL model election.pynb
  """
  from sppl.sym_util import binspace
  n = 4000
  param ~= randint(low=250, high=350)
  switch (param) cases (b in range(250, 350)):
    p ~= beta(a=277, b=b)
  switch (p) cases (x in binspace(0, 1, 20)):
    votes ~= binom(n=n, p=(x.left + x.right)/2)
  win ~= votes > 0.5*n
  """

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

/*

  * using binomial_dist/2
    var : param
    Probabilities (truncated):
    280: 0.0127000000000000
    334: 0.0124000000000000
    324: 0.0123000000000000
    307: 0.0121000000000000
    .........
    262: 0.0083000000000000
    314: 0.0082000000000000
    296: 0.0075000000000000
    267: 0.0073000000000000
    mean = 299.9221
    HPD intervals:
    HPD interval (0.94): 256.00000000000000..349.00000000000000

    var : p
    Probabilities (truncated):
    0.594754475851686: 0.0001000000000000
    0.588974544203882: 0.0001000000000000
    0.581122016286476: 0.0001000000000000
    0.579975564495532: 0.0001000000000000
    .........
    0.392598747078684: 0.0001000000000000
    0.39060755065802: 0.0001000000000000
    0.388060588168771: 0.0001000000000000
    0.386976901816159: 0.0001000000000000
    mean = 0.48141
    HPD intervals:
    HPD interval (0.94): 0.42630517751817..0.54228415898071

    var : votes
    Probabilities (truncated):
    1953: 0.0044000000000000
    1857: 0.0044000000000000
    1914: 0.0041000000000000
    1876: 0.0039000000000000
    .........
    1566: 0.0001000000000000
    1561: 0.0001000000000000
    1546: 0.0001000000000000
    1532: 0.0001000000000000
    mean = 1925.14
    HPD intervals:
    HPD interval (0.94): 1691.00000000000000..2168.00000000000000

    var : win
    Probabilities:
    false: 0.7114000000000000
    true: 0.2886000000000000
    mean = [false = 0.7114,true = 0.2886]
    HPD intervals:
    show_hpd_intervals: data is not numeric

*/
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() =>
  N = 4000,
  Param = 250 + random_integer(100),
  P = beta_dist(277,Param),
  % Votes = binomial_dist(N,P),
  % Votes = binomial_dist_smart(N,P), % faster than binomial_dist/2
  Votes = binomial_dist_btpe(N,P), % faster than binomial_dist/2
  
  Win = check(Votes > N*0.5),
  
  add("param",Param),
  add("p",P),
  add("votes",Votes),
  add("win",Win).