/* 

  Sqrt of max of 0..1 (is the same) in Picat.

  From Stand-up Maths
  "The unexpected probability result confusing everyone"
  https://www.youtube.com/watch?v=ga9Qk38FaHM

  The video shows that 
    (max (uniform 0 1) (uniform 0 1))
  is the same (i.e. behaves the same) as
    (sqrt (uniform 0 1))


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

/*

  var : diff
  Probabilities (truncated):
  0.92777558970511: 0.0001000000000000
  0.925259812430579: 0.0001000000000000
  0.916173837791425: 0.0001000000000000
  0.91605845457459: 0.0001000000000000
  .........
  -0.249998964174809: 0.0001000000000000
  -0.249999275683118: 0.0001000000000000
  -0.249999943109925: 0.0001000000000000
  -0.249999990554101: 0.0001000000000000
  mean = 0.00268449
  Percentiles:
  (0.001 -0.24999525271887879)
  (0.01 -0.24967659780529922)
  (0.025 -0.24749426122264953)
  (0.05 -0.2399787421285795)
  (0.1 -0.21789361397684093)
  (0.25 -0.1532682979951909)
  (0.5 -0.056375125808591942)
  (0.75 0.097472823319372903)
  (0.84 0.22148629709423551)
  (0.9 0.3336178811646473)
  (0.95 0.46735751941771292)
  (0.975 0.57558635062302943)
  (0.99 0.68633686285908357)
  (0.999 0.84847547466546236)
  (0.9999 0.92526006400830452)
  (0.99999 0.92752403713543274)

  var : v max
  Probabilities (truncated):
  0.999932869337468: 0.0001000000000000
  0.999878381844553: 0.0001000000000000
  0.999719490297008: 0.0001000000000000
  0.999719179235268: 0.0001000000000000
  .........
  0.016048504978441: 0.0001000000000000
  0.015230000957488: 0.0001000000000000
  0.013586132327833: 0.0001000000000000
  0.003950190732232: 0.0001000000000000
  mean = 0.670867
  Percentiles:
  (0.001 0.030909761711447388)
  (0.01 0.10933989699433554)
  (0.025 0.1599834571289753)
  (0.05 0.22837322830612455)
  (0.1 0.32164467257524126)
  (0.25 0.5099693031562349)
  (0.5 0.71012944109278231)
  (0.75 0.86944467219032573)
  (0.84 0.91707000449163378)
  (0.9 0.94916281027214733)
  (0.95 0.97482842392047331)
  (0.975 0.98751087650307956)
  (0.99 0.9954972558261348)
  (0.999 0.99942443348207721)
  (0.9999 0.99987838729330258)
  (0.99999 0.99992742113305133)

  var : v sqrt
  Probabilities (truncated):
  0.999966434105399: 0.0001000000000000
  0.999939189073292: 0.0001000000000000
  0.999859735311413: 0.0001000000000000
  0.999859579758712: 0.0001000000000000
  .........
  0.02256700103644: 0.0001000000000000
  0.016496889330736: 0.0001000000000000
  0.009220981047208: 0.0001000000000000
  0.006036620719949: 0.0001000000000000
  mean = 0.668183
  Percentiles:
  (0.001 0.037416409907645158)
  (0.01 0.10849453344382851)
  (0.025 0.16290387666264652)
  (0.05 0.23009153705968158)
  (0.1 0.32403901419504394)
  (0.25 0.50257153230609486)
  (0.5 0.7092381602711727)
  (0.75 0.86568145428667487)
  (0.84 0.9156571189432533)
  (0.9 0.94790053797527662)
  (0.95 0.97386603664375238)
  (0.975 0.98634654075050421)
  (0.99 0.99423415015283889)
  (0.999 0.99940581197464184)
  (0.9999 0.99993919179779545)
  (0.99999 0.99996370987463901)

  Note that the two 50th percentiles are 0.71012944109278231 and  
  0.7092381602711727 respectively, which are close to 
  sqrt(1/2) = 0.7071067811865476.


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,show_percentiles]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  X = uniform(0,1),
  Y = uniform(0,1),  

  VMax = max(X,Y),
  VSqrt = sqrt(X),

  Diff = VMax-VSqrt,
  add("v max",VMax),
  add("v sqrt",VSqrt),
  add("diff",Diff).


/*
  Using random_integer1(20) instead of uniform(0,1), 
  and for the sqrt version we have to use ceiling of sqrt(random1(N*N))

  n = 20
  var : diff
  Probabilities (truncated):
  0: 0.0689000000000000
  -1: 0.0622000000000000
  1: 0.0599000000000000
  -2: 0.0568000000000000
  .........
  -18: 0.0015000000000000
  18: 0.0006000000000000
  -19: 0.0004000000000000
  19: 0.0003000000000000
  mean = -0.0883
  Percentiles:
  (0.001 -18)
  (0.01 -15)
  (0.025 -13)
  (0.05 -11)
  (0.1 -9)
  (0.25 -5)
  (0.5 0)
  (0.75 4)
  (0.84 7)
  (0.9 9)
  (0.95 11)
  (0.975 13)
  (0.99 15)
  (0.999 17)
  (0.9999 19)
  (0.99999 19)

  var : v max
  Probabilities (truncated):
  20: 0.0960000000000000
  19: 0.0933000000000000
  18: 0.0886000000000000
  17: 0.0821000000000000
  .........
  4: 0.0184000000000000
  3: 0.0149000000000000
  2: 0.0077000000000000
  1: 0.0027000000000000
  mean = 13.7374
  Percentiles:
  (0.001 1)
  (0.01 2)
  (0.025 3)
  (0.05 5)
  (0.1 7)
  (0.25 10)
  (0.5 15)
  (0.75 18)
  (0.84 19)
  (0.9 19)
  (0.95 20)
  (0.975 20)
  (0.99 20)
  (0.999 20)
  (0.9999 20)
  (0.99999 20)

  var : v sqrt
  Probabilities (truncated):
  20: 0.1002000000000000
  19: 0.0917000000000000
  18: 0.0866000000000000
  17: 0.0805000000000000
  .........
  4: 0.0184000000000000
  3: 0.0115000000000000
  2: 0.0078000000000000
  1: 0.0040000000000000
  mean = 13.8257
  Percentiles:
  (0.001 1)
  (0.01 2)
  (0.025 4)
  (0.05 5)
  (0.1 7)
  (0.25 11)
  (0.5 15)
  (0.75 18)
  (0.84 19)
  (0.9 20)
  (0.95 20)
  (0.975 20)
  (0.99 20)
  (0.999 20)
  (0.9999 20)
  (0.99999 20)


*/

go2 ?=>
  N = 20,
  println(n=N),
  reset_store,
  run_model(10_000,$model2(N),[show_probs_trunc,mean,show_percentiles]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2(N) =>

  X = random_integer1(N),
  Y = random_integer1(N),
  Z = random_integer1(N*N),

  VMax = max(X,Y),
  VSqrt = ceiling(sqrt(Z)),

  Diff = VMax-VSqrt,
  add("v max",VMax),
  add("v sqrt",VSqrt),
  add("diff",Diff).