/* 

  Number of walks until no shoes in Picat.

  From 
  Gunnar Blom, Lars Holst, Dennis Sandell:
  "Problems and Snapshots from the World of Probability"
  Page 8f, Problem 1.6 Number of walks until no shoes

  """
  A has a house with one front door and one back door. He places
  two pairs of walking shoes at each door. For each walk, he selects
  one door at random, puts on a pair of shoes, returns after a walk
  to a randomly chosen door, and takes off the shoes at the door.
  We want to determine the expected number of finished walks until
  A discovers that no shoes are available at the door he has selected
  for his next walk.
  """

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

% The theoretical probability
%                       (- (+ (* 2 n) (* 4 n i)) (* 2 i i))
theoretical_prob(N) = (2*N) + (4*N*N) - (2*N*N).

/*

  [n = 2,theoretical = 12]
  var : a
  Probabilities (truncated):
  [right,right]: 0.0340000000000000
  [left,left]: 0.0336000000000000
  [left,left,left]: 0.0154000000000000
  [right,right,right]: 0.0142000000000000
  .........
  [left,left,left,left,left,left,left,left,right,right,right,left,right,left,right]: 0.0001000000000000
  [left,left,left,left,left,left,left,left,right,right,right,left,left,left,left,right,right,right,left,right,right,left,right,right,left,right,right,right,right,left,right,left,left,right,left,left,left,right,right]: 0.0001000000000000
  [left,left,left,left,left,left,left,left,right,left,right,left]: 0.0001000000000000
  [left,left,left,left,left,left,left,left]: 0.0001000000000000
  Mean (truncated)
  mean = [[right,right] = 0.034,[left,left] = 0.0336,[left,left,left] = 0.0154,[right,right,right] = 0.0142,[right,left,left] = 0.0098,[right,left,right] = 0.0077,[right,right,left] = 0.0075,[left,left,right] = 0.0073,[left,right,left] = 0.0071,[left,right,right] = 0.0067,[left,left,right,left] = 0.0066,[left,right,left,left] = 0.0064,[left,left,left,left] = 0.0062,[right,right,left,right] = 0.0061,[right,left,right,right] = 0.0057,[right,right,right,right] = 0.0054,[left,right,right,right] = 0.0051,[right,left,right,left] = 0.0047,[right,left,left,left] = 0.0045,[right,right,left,left] = 0.0044]

  var : len
  Probabilities (truncated):
  4: 0.0774000000000000
  3: 0.0757000000000000
  5: 0.0735000000000000
  2: 0.0676000000000000
  .........
  65: 0.0001000000000000
  64: 0.0001000000000000
  61: 0.0001000000000000
  57: 0.0001000000000000
  mean = 11.93


  [n = 3,theoretical = 24]
  var : a
  Probabilities (truncated):
  [left,left,left]: 0.0084000000000000
  [right,right,right]: 0.0079000000000000
  [left,left,left,left]: 0.0051000000000000
  [right,right,right,right]: 0.0048000000000000
  .........
  [left,left,left,left,left,left,left,left,right,left,right,right,left,left,left,right,right,left,right,left,left,left,right,right,right,left,left,right,right,right,left,right,right,right]: 0.0001000000000000
  [left,left,left,left,left,left,left,left,right,left,left,right,right,left,left,right,left,right,left,left,right,right,left,right,right,right,left,left,right,left,left,left,left,left,left,right,right,right,right,left,left,right,left,left,right,left,right,left]: 0.0001000000000000
  [left,left,left,left,left,left,left,left,left,right,left,left,right,right,right,left,right,left,left,left,left,right,right,left,right,right,left,left,left,left,left,right,right,right,right,right,left,left,left,left,right,left,left,left,right,left]: 0.0001000000000000
  [left,left,left,left,left,left,left,left,left,left,left,left,left]: 0.0001000000000000
  Mean (truncated)
  mean = [[left,left,left] = 0.0084,[right,right,right] = 0.0079,[left,left,left,left] = 0.0051,[right,right,right,right] = 0.0048,[right,right,right,right,right] = 0.0029,[right,right,right,left,right] = 0.0029,[left,left,left,left,left] = 0.0029,[left,right,left,left,left] = 0.0024,[left,left,left,right] = 0.0023,[left,left,left,right,left] = 0.0022,[right,left,left,left,left] = 0.0021,[right,left,right,right] = 0.002,[left,left,right,left,left] = 0.002,[left,left,left,left,right] = 0.002,[right,right,left,right] = 0.0019,[right,right,right,left] = 0.0018,[right,left,left,left] = 0.0018,[left,right,right,right] = 0.0018,[right,left,right,right,right,right] = 0.0017,[left,left,right,left] = 0.0017]

  var : len
  Probabilities (truncated):
  7: 0.0408000000000000
  6: 0.0404000000000000
  8: 0.0383000000000000
  9: 0.0376000000000000
  .........
  122: 0.0001000000000000
  120: 0.0001000000000000
  117: 0.0001000000000000
  113: 0.0001000000000000
  mean = 24.3467


  [n = 1,theoretical = 4]
  var : a
  Probabilities (truncated):
  [right]: 0.1261000000000000
  [left]: 0.1233000000000000
  [left,right]: 0.0639000000000000
  [right,left]: 0.0605000000000000
  .........
  [left,left,left,left,left,left,right,left]: 0.0001000000000000
  [left,left,left,left,left,left,right]: 0.0001000000000000
  [left,left,left,left,left,left,left]: 0.0001000000000000
  [left,left,left,left,left,left]: 0.0001000000000000
  Mean (truncated)
  mean = [[right] = 0.1261,[left] = 0.1233,[left,right] = 0.0639,[right,left] = 0.0605,[right,right] = 0.0316,[left,left] = 0.0289,[left,right,right] = 0.0239,[right,left,left] = 0.0232,[left,right,left] = 0.022,[right,left,right] = 0.0214,[left,left,right] = 0.0169,[right,right,left] = 0.0159,[right,left,left,right] = 0.0118,[right,left,right,left] = 0.0113,[left,right,right,right] = 0.0101,[right,left,left,left] = 0.009,[left,right,left,right] = 0.009,[left,right,right,left] = 0.0086,[right,left,right,right] = 0.0079,[left,right,left,left] = 0.0076]

  var : len
  Probabilities (truncated):
  1: 0.2494000000000000
  2: 0.1849000000000000
  3: 0.1362000000000000
  4: 0.1072000000000000
  .........
  28: 0.0002000000000000
  27: 0.0002000000000000
  35: 0.0001000000000000
  26: 0.0001000000000000
  mean = 4.0127


*/
go ?=>
  member(N,[2,3,1]),
  println([n=N,theoretical=theoretical_prob(N)]),
  reset_store,
  run_model(10_000,$model(N),[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,
  fail,
  nl.
go => true.


% Leave shoes at a random door
leave_shoes(A,Left,Right) = cond(flip() == true,
                                  take_shoes(A,Left+1,Right),
                                  take_shoes(A,Left,Right+1)).

% Take shoes from a random door
% Note: we only add to A when taking the shoes.
take_shoes(A,Left,Right) = Ret =>
  % Pick a door and check if there are any shoes
  Pick = cond(flip() == true,left,right),
  if (Pick == left, Left == 0) ;
     (Pick == right, Right == 0)
  then
    Ret = A
  else
    if Pick == left then
      Ret = leave_shoes(A ++ [Pick],Left-1,Right)
    else
      Ret = leave_shoes(A ++ [Pick],Left,Right-1)    
    end
  end.

model(N) =>  
  A = take_shoes([],N,N),
  add("a",A),
  add("len",A.len).