/* 

  Random walk in Picat.

  https://brainstellar.com/puzzles/probability/215
  """
  You are initially located at the origin on the x-axis. You start a 
  random walk with an equal probability of moving left or right, one 
  step at a time. What is the probability that you will reach point a
  before reaching point -b.

  Assume a and b are natural numbers. a > 1; b > 1.

  Answer: (b / (a + b))
  ...
  [I]t takes on average a*b steps to end the game.
  """

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

/*

  [a = 6,b = 2,theoretical_num_steps = 12,theoretical_last_value_a = 0.25]
  var : num steps
  Probabilities (truncated):
  2: 0.2476000000000000
  4: 0.1211000000000000
  6: 0.0925000000000000
  8: 0.0783000000000000
  .........
  106: 0.0001000000000000
  98: 0.0001000000000000
  92: 0.0001000000000000
  90: 0.0001000000000000
  mean = 12.1888

  var : p
  Probabilities:
  false: 0.7454000000000000
  true: 0.2546000000000000
  mean = [false = 0.7454,true = 0.2546]

  var : last value
  Probabilities:
  -2: 0.7454000000000000
  6: 0.2546000000000000
  mean = 0.0368


  [a = 9,b = 4,theoretical_num_steps = 36,theoretical_last_value_a = 0.307692]
  var : num steps
  Probabilities (truncated):
  4: 0.0670000000000000
  6: 0.0619000000000000
  8: 0.0616000000000000
  10: 0.0477000000000000
  .........
  171: 0.0001000000000000
  169: 0.0001000000000000
  163: 0.0001000000000000
  140: 0.0001000000000000
  mean = 35.9348

  var : p
  Probabilities:
  false: 0.7026000000000000
  true: 0.2974000000000000
  mean = [false = 0.7026,true = 0.2974]

  var : last value
  Probabilities:
  -4: 0.7026000000000000
  9: 0.2974000000000000
  mean = -0.1338

[a = 7,b = 2,theoretical_num_steps = 14,theoretical_last_value_a = 0.222222]
var : num steps
Probabilities (truncated):
2: 0.2453000000000000
4: 0.1224000000000000
6: 0.0787000000000000
8: 0.0534000000000000
.........
94: 0.0001000000000000
89: 0.0001000000000000
81: 0.0001000000000000
74: 0.0001000000000000
mean = 14.156
Histogram (total 10000)
  2.000:  2453 ############################################################ (0.245 / 0.245)
  5.650:  2083 ###################################################          (0.208 / 0.454)
  9.300:  1271 ###############################                              (0.127 / 0.581)
 12.950:   745 ##################                                           (0.074 / 0.655)
 16.600:   789 ###################                                          (0.079 / 0.734)
 20.250:   595 ###############                                              (0.059 / 0.794)
 23.900:   367 #########                                                    (0.037 / 0.830)
 27.550:   370 #########                                                    (0.037 / 0.867)
 31.200:   299 #######                                                      (0.030 / 0.897)
 34.850:   162 ####                                                         (0.016 / 0.913)
 38.500:   165 ####                                                         (0.017 / 0.930)
 42.150:   100 ##                                                           (0.010 / 0.940)
 45.800:   151 ####                                                         (0.015 / 0.955)
 49.450:    99 ##                                                           (0.010 / 0.965)
 53.100:    66 ##                                                           (0.007 / 0.972)
 56.750:    67 ##                                                           (0.007 / 0.978)
 60.400:    55 #                                                            (0.005 / 0.984)
 64.050:    20                                                              (0.002 / 0.986)
 67.700:    37 #                                                            (0.004 / 0.989)
 71.350:    16                                                              (0.002 / 0.991)
 75.000:    11                                                              (0.001 / 0.992)
 78.650:    15                                                              (0.002 / 0.994)
 82.300:    11                                                              (0.001 / 0.995)
 85.950:     8                                                              (0.001 / 0.996)
 89.600:    10                                                              (0.001 / 0.997)
 93.250:     9                                                              (0.001 / 0.997)
 96.900:     4                                                              (0.000 / 0.998)
100.550:     4                                                              (0.000 / 0.998)
104.200:     5                                                              (0.001 / 0.999)
107.850:     4                                                              (0.000 / 0.999)
111.500:     2                                                              (0.000 / 0.999)
115.150:     1                                                              (0.000 / 0.999)
118.800:     1                                                              (0.000 / 0.999)
122.450:     0                                                              (0.000 / 0.999)
126.100:     1                                                              (0.000 / 1.000)
129.750:     0                                                              (0.000 / 1.000)
133.400:     1                                                              (0.000 / 1.000)
137.050:     2                                                              (0.000 / 1.000)
140.700:     0                                                              (0.000 / 1.000)
144.350:     1                                                              (0.000 / 1.000)

var : p
Probabilities:
false: 0.7750000000000000
true: 0.2250000000000000
mean = [false = 0.775,true = 0.225]
Histogram (total 10000)
false     :  7750 ############################################################ (0.775 / 0.775)
true      :  2250 #################                                            (0.225 / 1.000)

var : last value
Probabilities:
-2: 0.7750000000000000
7: 0.2250000000000000
mean = 0.025
Histogram (total 10000)
     -2:  7750 ############################################################ (0.775 / 0.775)
      7:  2250 #################                                            (0.225 / 1.000)

*/
go ?=>
  _ = random2(),
  A = random_integer1(10),
  B = random_integer1(10),
  println([a=A,b=B,theoretical_num_steps=(A*B),theoretical_last_value_a=(B/(A+B))]),
  reset_store,
  run_model(10_000,$model(A,B),[show_probs_trunc,mean,show_histogram]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  

f(Lst,A,B) = Res =>
  if membchk(A,Lst) ; membchk(-B,Lst) then
    Res = Lst
  else
    Res = f(Lst ++ [Lst.last + uniform_draw([-1,1])],A,B)
  end.

model(A,B) =>
  Lst = f([0],A,B),
  NumSteps = Lst.len - 1,
  LastValue = Lst.last,
  P = check(LastValue > 0),
  add("num steps",NumSteps),  
  add("p",P),
  add("last value",LastValue).