/* 

  Birthday coincidence in Picat.

  From https://www.doc.ic.ac.uk/~mpd37/teaching/2014/ml_tutorials/2014-10-15-wood-anglican.pdf
  page 18f:
  """
  Birthday Coincidence
  Approximately, what’s the probability that in a room
  filled with 23 people at least one pair of people
  have the same birthday? 

  (assume birthday (mem (lambda (i) (uniform-discrete 1 366))))
  (assume N 23)
  (assume pair-equal
    (lambda (i j)
      (if (> i N)
        false
          (if (> j N)
            (pair-equal (+ i 1) (+ i 2))
            (if (= (birthday i) (birthday j))
              true
              (pair-equal i (+ j 1)))))))
  (predict (pair-equal 1 2)
  """

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

/*
  (Slighly edited)

  2 = 0.004
  3 = 0.006
  4 = 0.0185
  5 = 0.022
  6 = 0.0435
  7 = 0.053
  8 = 0.0695
  9 = 0.0995
  10 = 0.13
  11 = 0.133
  12 = 0.174
  13 = 0.2
  14 = 0.2165
  15 = 0.249
  16 = 0.2735
  17 = 0.304
  18 = 0.342
  19 = 0.3725
  20 = 0.413
  21 = 0.4445
  22 = 0.487
  23 = 0.5045
  24 = 0.532
  25 = 0.5615
  26 = 0.606
  27 = 0.6415
  28 = 0.645
  29 = 0.6855
  30 = 0.723

*/
go ?=>
  member(N,2..30),
  reset_store,
  run_model(2000,$model(N),[]),
  % The the probability of at least one common birthday
  println(N=get_store().get("p").get_probs().new_map().get(true)),
  fail,
  nl.
go => true.
  

pair_equal(I,J,N,Birthdays) = Res =>
  if I > N then
    Res = false
  else
    if J > N then
      Res = pair_equal(I + 1, 2 + I,N,Birthdays)
    else
      if Birthdays[I] == Birthdays[J] then
        Res = true
      else
        Res = pair_equal(I,J+1,N,Birthdays)
      end
    end
  end.

model(N) =>
  Birthdays = random_integer1_n(365,N),
  P = pair_equal(1,2,N,Birthdays),
  add("p",P).


/* 

  A neater, and more Picat-ish, model 

  2 = 0.0025
  3 = 0.009
  4 = 0.0215
  5 = 0.0265
  6 = 0.0445
  7 = 0.0565
  8 = 0.0655
  9 = 0.0895
  10 = 0.1205
  11 = 0.1385
  12 = 0.1485
  13 = 0.1945
  14 = 0.227
  15 = 0.2715
  16 = 0.2745
  17 = 0.3235
  18 = 0.3365
  19 = 0.362
  20 = 0.4265
  21 = 0.45
  22 = 0.4755
  23 = 0.528
  24 = 0.5445
  25 = 0.556
  26 = 0.5875
  27 = 0.6265
  28 = 0.651
  29 = 0.68
  30 = 0.703

*/ 
go2 ?=>
  member(N,2..30),
  reset_store,
  run_model(2000,$model(N),[]),
  println(N=get_store().get("p").get_probs().new_map().get(true)),
  fail,
  nl.
go2 => true.
  
model2(N) =>
  Birthdays = random_integer1_n(365,N),
  P = false,
  foreach(I in 1..N, J in I+1..N, break(B==true))
    if Birthdays[I] == Birthdays[J] then
      P := true
    end
  end,
  add("p",P).