/* 

  Bullets of fate in Picat.

  https://brainstellar.com/puzzles/probability/1008
  """
  In a room stand N armed and angry people. At each chime of a clock, everyone 
  simultaneously spins around and shoots a random other person. The persons who 
  get shot, fall dead and the survivors spin and shoot again at the next chime. 
  Eventually, either everyone is dead or there is a single survivor.

  As N grows, what is the limiting probability that there will be a survivor?

  Warning: I could not solve it. The solution only shows what I tried.
  """

  The limiting probability is around 0.5.

  Cf
  https://colab.research.google.com/gist/varun-seth/c57ed32726bf7492daa9dea375398c09/puzzle_1008_rdd_exact.ipynb

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

/*

  n = 3
  var : p
  Probabilities:
  false: 0.5647000000000000
  true: 0.4353000000000000
  mean = [false = 0.5647,true = 0.4353]


  n = 4
  var : p
  Probabilities:
  false: 0.5192000000000000
  true: 0.4808000000000000
  mean = [false = 0.5192,true = 0.4808]


  n = 5
  var : p
  Probabilities:
  true: 0.5612000000000000
  false: 0.4388000000000000
  mean = [true = 0.5612,false = 0.4388]


  n = 6
  var : p
  Probabilities:
  true: 0.5288000000000000
  false: 0.4712000000000000
  mean = [true = 0.5288,false = 0.4712]


  n = 7
  var : p
  Probabilities:
  false: 0.5179000000000000
  true: 0.4821000000000000
  mean = [false = 0.5179,true = 0.4821]


  n = 10
  var : p
  Probabilities:
  false: 0.5098000000000000
  true: 0.4902000000000000
  mean = [false = 0.5098,true = 0.4902]


  n = 50
  var : p
  Probabilities:
  false: 0.5270000000000000
  true: 0.4730000000000000
  mean = [false = 0.527,true = 0.473]


  n = 100
  var : p
  Probabilities:
  true: 0.5213000000000000
  false: 0.4787000000000000
  mean = [true = 0.5213,false = 0.4787]

*/
go ?=>
  member(N,[3,4,5,6,7,10,50,100]),
  println(n=N),
  reset_store,
  run_model(10_000,$model(N),[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.
  

f(A) = Res =>
  if A.len <= 1 then
    Res = A
  else
    % Don't forget to sort the sets when using ordsets
    Targets = [uniform_draw(A.delete(P)) : P in A].sort,
    Res = f(subtract(A,Targets).sort)
  end.
  
model(N) =>
  People = draw_without_replacement(1..N),  
  A = f(People),
  Survivor = cond(A.len > 0,A.first,0), % 0: No survivor
  P = check(Survivor != 0),
  % add("survivor",Survivor),  
  add("p",P).