/* 

  Always take the middle taxi in Picat.

  From 
  https://www.janestreet.com/static/pdfs/real-numbers/RN_PDF_02.pdf
  """
  In the busy city center of Metropolis, exactly half of the cars are yellow taxis. You are
  trying to catch a taxi with your friend, who is superstitious and will only ride in a taxi that
  is both directly behind another taxi and directly in front of another taxi. If a stoplight in
  Metropolis has six cars stopped at it, and each car is a taxi independently with probability
  1/2, what is the probability that there are three yellow taxis in a row at the traffic light?
  """

  Via Pascal Bercker: "Jane Street — Problem Of The Week #2 — Take the Taxi in the Middle"
  https://medium.com/@pbercker/jane-street-problem-of-the-week-2-take-the-taxi-in-the-middle-fb5cf988cb84

  I also wrote a Gamble model: gamble_always_take_the_middle_taxi.pi.
  which gave the exact probabilities:

  variable : s
  0: 11/16 (0.6875)
  1: 3/16 (0.1875)
  2: 5/64 (0.078125)
  3: 1/32 (0.03125)
  4: 1/64 (0.015625)
  mean: 1/2 (0.5)

  variable : p
  #f: 11/16 (0.6875)
  #t: 5/16 (0.3125)
  mean: 5/16 (0.3125)

  variable : p2
  #f: 13/16 (0.8125)
  #t: 3/16 (0.1875)
  mean: 3/16 (0.1875)


  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 : s
  Probabilities:
  0: 0.6901000000000000
  1: 0.1865000000000000
  2: 0.0775000000000000
  3: 0.0313000000000000
  4: 0.0146000000000000
  mean = 0.4938

  var : p
  Probabilities:
  false: 0.6901000000000000
  true: 0.3099000000000000
  mean = [false = 0.6901,true = 0.3099]

  var : 2
  Probabilities:
  false: 0.8135000000000000
  true: 0.1865000000000000
  mean = [false = 0.8135,true = 0.1865]

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean % ,
                           % show_percentiles,show_histogram,
                           % show_hpd_intervals,hpd_intervals=[0.84],
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  N = 6, % number of cars
  TaxiP = 1/2, % probability of a taxi

  Cars = bern_n(TaxiP,N),

  S = count_occurrences([1,1,1],sublists(Cars,3)), % [1 : T in sublists(Cars,3), T == [1,1,1]].sum,

  % Probability that there are at least 1 occurrence of 3 taxis in a row
  P = check(S > 0),

  % Probability that there are exactly 1 occurrence of 3 taxis in a row.
  P2 = check(S == 1),

  % add("cars",Cars),
  add("s",S),
  add("p",P),
  add("2",P2).