/* 

  Probability of Alien Extinction in Picat.

  From MindYourDecisions
  """
  A single alien lands on Earth. Every day after that, each alien on Earth
  undergoes a transformation, which could be any of the four equally likely
  events:
  (a) the alien dies
  (b) the alien does nothing
  (c) the alien replicates itself (2 aliens total)
  (d) the alien replicates itself twice (3 alients total)

  What is the probability the alien race eventually dies out
  and goes extinct?
  """

  (Via Pascal Bercker)

  Cf my Gamble model gamble_alien_extinction_puzzle.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 : a
  Probabilities:
  0: 0.2503000000000000
  3: 0.0948000000000000
  2: 0.0743000000000000
  4: 0.0655000000000000
  5: 0.0621000000000000
  1: 0.0610000000000000
  6: 0.0509000000000000
  7: 0.0441000000000000
  8: 0.0408000000000000
  9: 0.0348000000000000
  10: 0.0307000000000000
  11: 0.0271000000000000
  12: 0.0233000000000000
  13: 0.0197000000000000
  14: 0.0162000000000000
  15: 0.0144000000000000
  16: 0.0125000000000000
  17: 0.0098000000000000
  18: 0.0093000000000000
  20: 0.0067000000000000
  19: 0.0067000000000000
  21: 0.0061000000000000
  22: 0.0049000000000000
  23: 0.0047000000000000
  24: 0.0040000000000000
  25: 0.0038000000000000
  26: 0.0030000000000000
  27: 0.0025000000000000
  28: 0.0023000000000000
  29: 0.0018000000000000
  34: 0.0015000000000000
  32: 0.0015000000000000
  33: 0.0013000000000000
  31: 0.0011000000000000
  30: 0.0011000000000000
  36: 0.0008000000000000
  35: 0.0008000000000000
  39: 0.0007000000000000
  37: 0.0006000000000000
  46: 0.0003000000000000
  45: 0.0003000000000000
  38: 0.0003000000000000
  50: 0.0002000000000000
  48: 0.0002000000000000
  44: 0.0002000000000000
  43: 0.0002000000000000
  41: 0.0002000000000000
  40: 0.0002000000000000
  58: 0.0001000000000000
  52: 0.0001000000000000
  49: 0.0001000000000000
  42: 0.0001000000000000
  mean = 5.9235

  var : p
  Probabilities:
  false: 0.7497000000000000
  true: 0.2503000000000000
  mean = 0.2503

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs,mean]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  

prob() = uniform_draw(0..3).

f(NumAliens,Steps,Limit) = Res =>
  if NumAliens == 0 ; Steps >= Limit then
    Res = NumAliens
  else
    A = prob(),
    if A == 0 then
      Res = 0
    else
      Res = A + sum([f(1,Steps+1,Limit) : _ in 1..NumAliens])
    end
  end.

model() =>
  Limit = 100,
  A = f(1,0,Limit),
  P = check(A == 0),
  
  add("a",A),
  add("p",P).