/* 

  Rolling multiple dice and picking the highest value in Picat.

  https://www.youtube.com/watch?v=X_DdGRjtwAo&t=157s
  The unexpected logic behind rolling multiple dice and picking the highest

  Following Matt Parker's simulations.

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

/*
  @11:10 Parker shows that the formula for two dice is 
    1/(6*N)*(N+1)*(4*N-1)
  where N is the sides of the die.
*/
theoretical_2d(N) = (1 / (6*N)) * (N+1) * ((4*N) - 1).

/*

  * For two d(N) dice.

  [n = 20,m = 2,theoretical_2d = 13.824999999999999 (553 / 40)]
  var : v
  Probabilities:
  19: 0.0963000000000000
  20: 0.0954000000000000
  17: 0.0846000000000000
  18: 0.0820000000000000
  16: 0.0791000000000000
  15: 0.0752000000000000
  13: 0.0677000000000000
  14: 0.0663000000000000
  11: 0.0556000000000000
  12: 0.0542000000000000
  10: 0.0466000000000000
  9: 0.0429000000000000
  8: 0.0338000000000000
  7: 0.0324000000000000
  6: 0.0257000000000000
  5: 0.0226000000000000
  4: 0.0167000000000000
  3: 0.0130000000000000
  2: 0.0079000000000000
  1: 0.0020000000000000
  mean = 13.8559

  (Matt Parker's first Python simulation gives 13.829135)

  [n = 12,m = 2,theoretical_2d = 8.486111111111111 (611 / 72)]
  var : v
  Probabilities:
  12: 0.1570000000000000
  11: 0.1500000000000000
  10: 0.1325000000000000
  9: 0.1168000000000000
  8: 0.1050000000000000
  7: 0.0896000000000000
  6: 0.0799000000000000
  5: 0.0623000000000000
  4: 0.0434000000000000
  3: 0.0344000000000000
  2: 0.0232000000000000
  1: 0.0059000000000000
  mean = 8.4974


  [n = 6,m = 2,theoretical_2d = 4.472222222222221 (161 / 36)]
  var : v
  Probabilities:
  6: 0.3113000000000000
  5: 0.2447000000000000
  4: 0.1930000000000000
  3: 0.1382000000000000
  2: 0.0855000000000000
  1: 0.0273000000000000
  mean = 4.4762


  [n = 120,m = 2,theoretical_2d = 80.498611111111117 (57959 / 720)]
  var : v
  Probabilities (truncated):
  117: 0.0174000000000000
  118: 0.0173000000000000
  114: 0.0170000000000000
  119: 0.0165000000000000
  .........
  7: 0.0005000000000000
  4: 0.0004000000000000
  3: 0.0001000000000000
  2: 0.0001000000000000
  mean = 80.815


  * For three dice

  [n = 6,m = 3]
  var : v
  Probabilities:
  6: 0.4233000000000000
  5: 0.2866000000000000
  4: 0.1698000000000000
  3: 0.0830000000000000
  2: 0.0328000000000000
  1: 0.0045000000000000
  mean = 4.9711


  [n = 12,m = 3]
  var : v
  Probabilities:
  12: 0.2338000000000000
  11: 0.1878000000000000
  10: 0.1522000000000000
  9: 0.1228000000000000
  8: 0.0963000000000000
  7: 0.0773000000000000
  6: 0.0568000000000000
  5: 0.0361000000000000
  4: 0.0223000000000000
  3: 0.0106000000000000
  2: 0.0035000000000000
  1: 0.0005000000000000
  mean = 9.4599


  [n = 20,m = 3]
  var : v
  Probabilities:
  20: 0.1441000000000000
  19: 0.1295000000000000
  18: 0.1137000000000000
  17: 0.0990000000000000
  16: 0.0912000000000000
  15: 0.0821000000000000
  14: 0.0670000000000000
  13: 0.0593000000000000
  12: 0.0481000000000000
  11: 0.0422000000000000
  10: 0.0351000000000000
  9: 0.0276000000000000
  8: 0.0190000000000000
  7: 0.0170000000000000
  6: 0.0109000000000000
  5: 0.0074000000000000
  4: 0.0031000000000000
  3: 0.0024000000000000
  2: 0.0013000000000000
  mean = 15.5081


  [n = 120,m = 3]
  var : v
  Probabilities (truncated):
  118: 0.0279000000000000
  116: 0.0253000000000000
  114: 0.0234000000000000
  119: 0.0230000000000000
  .........
  10: 0.0001000000000000
  9: 0.0001000000000000
  8: 0.0001000000000000
  4: 0.0001000000000000
  mean = 90.1547


*/
go ?=>
  member([N,M],[[20,2],
                [12,2],
                [6,2],
                [120,2],
                [6,3],
                [12,3],
                [20,3],
                [120,3]]),

  if M == 2 then
    println([n=N,m=M,theoretical_2d=theoretical_2d(N).and_rat])
  else
    println([n=N,m=M])  
  end,
  reset_store,
  Options = cond(N <= 20,[show_probs,mean],[show_probs_trunc,mean]),
  run_model(10_000,$model(N,M),Options),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.
  
  
model(N,M) =>
  % Throw M dice
  Throws = dice_n(N,M),

  % What is the max value
  V = Throws.max,
  
  add("v",V).


/*
  Theoretical values for 2 dice:

  n: 1   theoretical (2 dice): 1.0000000000 (1 / 1)
  n: 2   theoretical (2 dice): 1.7500000000 (7 / 4)
  n: 3   theoretical (2 dice): 2.4444444444 (22 / 9)
  n: 4   theoretical (2 dice): 3.1250000000 (25 / 8)
  n: 5   theoretical (2 dice): 3.8000000000 (19 / 5)
  n: 6   theoretical (2 dice): 4.4722222222 (161 / 36)
  n: 7   theoretical (2 dice): 5.1428571429 (36 / 7)
  n: 8   theoretical (2 dice): 5.8125000000 (93 / 16)
  n: 9   theoretical (2 dice): 6.4814814815 (175 / 27)
  n: 10  theoretical (2 dice): 7.1500000000 (143 / 20)
  n: 11  theoretical (2 dice): 7.8181818182 (86 / 11)
  n: 12  theoretical (2 dice): 8.4861111111 (611 / 72)
  n: 13  theoretical (2 dice): 9.1538461538 (119 / 13)
  n: 14  theoretical (2 dice): 9.8214285714 (275 / 28)
  n: 15  theoretical (2 dice): 10.4888888889 (472 / 45)
  n: 16  theoretical (2 dice): 11.1562500000 (357 / 32)
  n: 17  theoretical (2 dice): 11.8235294118 (201 / 17)
  n: 18  theoretical (2 dice): 12.4907407407 (1349 / 108)
  n: 19  theoretical (2 dice): 13.1578947368 (250 / 19)
  n: 20  theoretical (2 dice): 13.8250000000 (553 / 40)
  n: 120 theoretical (2 dice): 80.4986111111 (57959 / 720)

*/
go2 => 
  foreach(N in 1..20++[120])
    T = theoretical_2d(N),
    printf("n: %-3w theoretical (2 dice): %0.10f (%w)\n",N,T,T.to_rat)
  end,
  nl.