/* 

  Runs in Picat.

  Some experiments with runs.

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

% Loop until the runs of xs are "settled", normally a single number
% Returns the last xs and the number of iterations
runs_loop(A) = runs_loop1(A,0,false).
runs_loop(A,Trace) = runs_loop1(A,0,Trace).
runs_loop1(A,I,Trace) = Res =>
  if Trace == true then
    println("runs_loop"=A=A.len)
  end,
  NewA = get_runs_lens(A),
  if A == NewA then
    Res = [A,I]
  else
    Res = runs_loop1(NewA,I+1,Trace)
  end.


/*
  Here is an example of runs-loop on a binary list of 200 elements. 
  This evaluates the given list with get-runs until the list is settled (does not change).

  [1,1,0,1,1,0,1,0,1,0,1,1,1,0,0,1,0,0,1,0,0,1,1,1,1,0,0,1,1,0,0,0,1,0,0,0,0,1,1,1,1,0,0,1,1,1,0,1,1,1,1,0,0,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,1,0,0,1,1,1,0,0,0,0,1,1,1,0,1,1,1,1,0,1,0,1,1,0,0,0,1,1,0,0,0,0,1,1,1,0,0,1,0,0,1,1,1,0,1,1,1,0,0,1,0,0,1,1,0,1,1,1,1,0,0,1,0,1,0,1,1,0,0,1,0,1,0,0,1,1,1,0,0,1,1,0,1,0,1,0,0,1,1,1,1,1,0,1,0,1,0,1,1,0,0,0,1,0,0,0,0,1,0,0,0,1,0,1,0,1,1,0,0,1,1,1,0,1,0,1,0,0,0,0,0,1]
  get_runs = [[1,1],[0],[1,1],[0],[1],[0],[1],[0],[1,1,1],[0,0],[1],[0,0],[1],[0,0],[1,1,1,1],[0,0],[1,1],[0,0,0],[1],[0,0,0,0],[1,1,1,1],[0,0],[1,1,1],[0],[1,1,1,1],[0,0,0],[1,1],[0],[1,1],[0],[1,1],[0],[1,1],[0],[1,1,1],[0,0],[1,1,1],[0,0,0,0],[1,1,1],[0],[1,1,1,1],[0],[1],[0],[1,1],[0,0,0],[1,1],[0,0,0,0],[1,1,1],[0,0],[1],[0,0],[1,1,1],[0],[1,1,1],[0,0],[1],[0,0],[1,1],[0],[1,1,1,1],[0,0],[1],[0],[1],[0],[1,1],[0,0],[1],[0],[1],[0,0],[1,1,1],[0,0],[1,1],[0],[1],[0],[1],[0,0],[1,1,1,1,1],[0],[1],[0],[1],[0],[1,1],[0,0,0],[1],[0,0,0,0],[1],[0,0,0],[1],[0],[1],[0],[1,1],[0,0],[1,1,1],[0],[1],[0],[1],[0,0,0,0,0],[1]]
  [2,1,2,1,1,1,1,1,3,2,1,2,1,2,4,2,2,3,1,4,4,2,3,1,4,3,2,1,2,1,2,1,2,1,3,2,3,4,3,1,4,1,1,1,2,3,2,4,3,2,1,2,3,1,3,2,1,2,2,1,4,2,1,1,1,1,2,2,1,1,1,2,3,2,2,1,1,1,1,2,5,1,1,1,1,1,2,3,1,4,1,3,1,1,1,1,2,2,3,1,1,1,1,5,1]
  runs_loop = [1,1,0,1,1,0,1,0,1,0,1,1,1,0,0,1,0,0,1,0,0,1,1,1,1,0,0,1,1,0,0,0,1,0,0,0,0,1,1,1,1,0,0,1,1,1,0,1,1,1,1,0,0,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,1,0,0,1,1,1,0,0,0,0,1,1,1,0,1,1,1,1,0,1,0,1,1,0,0,0,1,1,0,0,0,0,1,1,1,0,0,1,0,0,1,1,1,0,1,1,1,0,0,1,0,0,1,1,0,1,1,1,1,0,0,1,0,1,0,1,1,0,0,1,0,1,0,0,1,1,1,0,0,1,1,0,1,0,1,0,0,1,1,1,1,1,0,1,0,1,0,1,1,0,0,0,1,0,0,0,0,1,0,0,0,1,0,1,0,1,1,0,0,1,1,1,0,1,0,1,0,0,0,0,0,1] = 200
  runs_loop = [2,1,2,1,1,1,1,1,3,2,1,2,1,2,4,2,2,3,1,4,4,2,3,1,4,3,2,1,2,1,2,1,2,1,3,2,3,4,3,1,4,1,1,1,2,3,2,4,3,2,1,2,3,1,3,2,1,2,2,1,4,2,1,1,1,1,2,2,1,1,1,2,3,2,2,1,1,1,1,2,5,1,1,1,1,1,2,3,1,4,1,3,1,1,1,1,2,2,3,1,1,1,1,5,1] = 105
  runs_loop = [1,1,1,5,1,1,1,1,1,1,1,2,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,4,2,3,1,1,2,4,1,1,5,1,1,1,1,1,1,4,2,1,4,1,1] = 75
  runs_loop = [3,1,7,1,2,1,20,1,13,1,3,1,1,1,2,1,1,2,1,6,1,1,1,1,2] = 25
  runs_loop = [1,1,1,1,1,1,1,1,1,1,1,3,1,2,1,1,1,4,1] = 19
  runs_loop = [11,1,1,1,3,1,1] = 7
  runs_loop = [1,3,1,2] = 4
  runs_loop = [1,1,1,1] = 4
  runs_loop = [4] = 1
  runs_loop = [1] = 1
  [[1],9]

  I.e. it took 9 iterations to settle down. And it seems to always end with a single 1.
*/
go ?=>
  _ = random2(),
  Throws = random_integer_n(2,200), % 200 0/1s
  % Throws = random_integer_n(2,200), % 200 0/1s
  % Throws = [1,0,1,0,1,0,1,0],
  println(Throws),
  println(get_runs=get_runs(Throws)),
  println(get_runs_lens(Throws)),
  println(runs_loop(Throws,true)),
  nl.
go => true.


/*
  Number of iterations for a repeated run on a binary list of different lengths (N)

  [n = 10,m = 2]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]

  var : num iterations
  Probabilities:
  6: 0.3356000000000000
  4: 0.2562000000000000
  5: 0.2440000000000000
  7: 0.1521000000000000
  3: 0.0058000000000000
  8: 0.0045000000000000
  2: 0.0018000000000000
  mean = 5.3801


  [n = 200,m = 2]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]

  var : num iterations
  Probabilities:
  10: 0.3572000000000000
  11: 0.2758000000000000
  9: 0.1576000000000000
  8: 0.1426000000000000
  12: 0.0626000000000000
  6: 0.0016000000000000
  13: 0.0013000000000000
  7: 0.0013000000000000
  mean = 9.9518


  [n = 1000,m = 2]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]

  var : num iterations
  Probabilities:
  12: 0.3543000000000000
  13: 0.3136000000000000
  14: 0.1167000000000000
  11: 0.1103000000000000
  10: 0.1010000000000000
  15: 0.0036000000000000
  9: 0.0003000000000000
  8: 0.0002000000000000
  mean = 12.2438



*/
go2 ?=>
  % member(M,[2,3]),
  M = 2,
  member(N,[10,200,1000]),
  println([n=N,m=M]),
  reset_store,
  run_model(10_000,$model2(N,M),[show_probs_trunc,mean % ,
                           % show_percentiles,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go2 => true.


model2(N,M) =>
  Throws = random_integer_n(M,N),

  [LastRun,NumIterations] = runs_loop(Throws,false),
  
  add("last run",LastRun),
  add("num iterations",NumIterations).


/*
  Same length (N=200) but with different number of possible integers (M=2..10).
  The mean of num_iterations are - perhaps unsurprisingly - decreasing when
  increasing the number of different values.

  [n = 200,m = 2]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : num iterations
  Probabilities:
  10: 0.3513000000000000
  11: 0.2786000000000000
  9: 0.1587000000000000
  8: 0.1414000000000000
  12: 0.0667000000000000
  6: 0.0016000000000000
  13: 0.0012000000000000
  7: 0.0005000000000000
  mean = 9.9662
  HPD intervals:
  HPD interval (0.94): 8.00000000000000..12.00000000000000


  [n = 200,m = 3]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : num iterations
  Probabilities:
  10: 0.3042000000000000
  11: 0.2915000000000000
  9: 0.2701000000000000
  12: 0.0950000000000000
  8: 0.0232000000000000
  7: 0.0134000000000000
  13: 0.0026000000000000
  mean = 10.1326
  HPD intervals:
  HPD interval (0.94): 9.00000000000000..12.00000000000000


  [n = 200,m = 4]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : num iterations
  Probabilities:
  9: 0.3753000000000000
  10: 0.2825000000000000
  11: 0.1990000000000000
  7: 0.0901000000000000
  12: 0.0404000000000000
  8: 0.0119000000000000
  13: 0.0008000000000000
  mean = 9.6128
  HPD intervals:
  HPD interval (0.94): 7.00000000000000..11.00000000000000


  [n = 200,m = 5]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : num iterations
  Probabilities:
  9: 0.3787000000000000
  10: 0.2687000000000000
  7: 0.2015000000000000
  11: 0.1180000000000000
  8: 0.0176000000000000
  12: 0.0144000000000000
  5: 0.0010000000000000
  13: 0.0001000000000000
  mean = 9.1237
  HPD intervals:
  HPD interval (0.94): 7.00000000000000..11.00000000000000


  [n = 200,m = 6]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : num iterations
  Probabilities:
  7: 0.3376000000000000
  9: 0.3244000000000000
  10: 0.2264000000000000
  11: 0.0691000000000000
  8: 0.0332000000000000
  5: 0.0056000000000000
  12: 0.0037000000000000
  mean = 8.6449
  HPD intervals:
  HPD interval (0.94): 7.00000000000000..11.00000000000000


  [n = 200,m = 7]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : num iterations
  Probabilities:
  7: 0.4426000000000000
  9: 0.2759000000000000
  10: 0.1749000000000000
  8: 0.0471000000000000
  11: 0.0408000000000000
  5: 0.0174000000000000
  12: 0.0013000000000000
  mean = 8.2585
  HPD intervals:
  HPD interval (0.94): 7.00000000000000..10.00000000000000


  [n = 200,m = 8]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : num iterations
  Probabilities:
  7: 0.4974000000000000
  9: 0.2321000000000000
  10: 0.1386000000000000
  8: 0.0583000000000000
  5: 0.0503000000000000
  11: 0.0226000000000000
  12: 0.0007000000000000
  mean = 7.9316
  HPD intervals:
  HPD interval (0.94): 7.00000000000000..11.00000000000000


  [n = 200,m = 9]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : num iterations
  Probabilities:
  7: 0.5489000000000001
  9: 0.1849000000000000
  10: 0.1046000000000000
  5: 0.0810000000000000
  8: 0.0692000000000000
  11: 0.0112000000000000
  12: 0.0002000000000000
  mean = 7.6366
  HPD intervals:
  HPD interval (0.94): 5.00000000000000..10.00000000000000


  [n = 200,m = 10]
  var : last run
  Probabilities:
  [1]: 1.0000000000000000
  mean = [[1] = 1.0]
  HPD intervals:
  show_hpd_intervals: data is not numeric

  var : num iterations
  Probabilities:
  7: 0.5530000000000000
  9: 0.1587000000000000
  5: 0.1216000000000000
  8: 0.0814000000000000
  10: 0.0784000000000000
  11: 0.0068000000000000
  12: 0.0001000000000000
  mean = 7.4185
  HPD intervals:
  HPD interval (0.94): 5.00000000000000..10.00000000000000


*/
go3 ?=>
  member(M,2..10),
  N = 200,
  println([n=N,m=M]),
  reset_store,
  run_model(10_000,$model2(N,M),[show_probs_trunc,mean,
                           % show_percentiles,
                           show_hpd_intervals,hpd_intervals=[0.94]
                           % show_histogram,
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go3 => true.