/* 

  Breaking a stick in three pieces in Picat.

  https://math.stackexchange.com/questions/3872156/uniformly-at-random-break-a-unit-stick-in-two-places-what-is-the-probability-t
  """
  Uniformly at random, break a unit stick in two places. What is the probability that 
  the smallest piece is <= 1/5?

  I was asked this in an interview and wasn't sure how to solve it:

     Consider a stick of length 1
     * Select two points independently and uniformly at random on the stick. Break the stick at these 
       two points, resulting in 3 smaller pieces. What is the probability that the smallest of these 
       pieces is <= 1/5?
  """

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

/*

  var : b1
  Probabilities (truncated):
  0.996904653029938: 0.0002004409701343
  0.989278938150629: 0.0002004409701343
  0.969095710184935: 0.0002004409701343
  0.968248281613108: 0.0002004409701343
  .........
  0.000197638757619: 0.0002004409701343
  0.000175565481268: 0.0002004409701343
  0.000147870276192: 0.0002004409701343
  0.000142120290614: 0.0002004409701343
  mean = 0.32864

  var : b2
  Probabilities (truncated):
  0.997616666833692: 0.0002004409701343
  0.993859708306314: 0.0002004409701343
  0.979928455306184: 0.0002004409701343
  0.979143430469159: 0.0002004409701343
  .........
  0.00022991094749: 0.0002004409701343
  0.000141225289619: 0.0002004409701343
  0.000134701374981: 0.0002004409701343
  0.000036415644007: 0.0002004409701343
  mean = 0.33478

  var : piece1
  Probabilities (truncated):
  0.489010192215913: 0.0002004409701343
  0.488447134144812: 0.0002004409701343
  0.484844512066266: 0.0002004409701343
  0.481784746740844: 0.0002004409701343
  .........
  0.000142120290614: 0.0002004409701343
  0.000141225289619: 0.0002004409701343
  0.000134701374981: 0.0002004409701343
  0.000036415644007: 0.0002004409701343
  mean = 0.165083

  var : piece2
  Probabilities (truncated):
  0.997616666833692: 0.0002004409701343
  0.996904653029938: 0.0002004409701343
  0.993859708306314: 0.0002004409701343
  0.989278938150629: 0.0002004409701343
  .........
  0.017573696569341: 0.0002004409701343
  0.007526851728338: 0.0002004409701343
  0.001791840420008: 0.0002004409701343
  0.000522788148617: 0.0002004409701343
  mean = 0.498338

  var : piece3
  Probabilities (truncated):
  0.999342510476402: 0.0002004409701343
  0.996588667387417: 0.0002004409701343
  0.990109856701507: 0.0002004409701343
  0.974653874046474: 0.0002004409701343
  .........
  0.000222559552743: 0.0002004409701343
  0.000090186949861: 0.0002004409701343
  0.000076612457669: 0.0002004409701343
  0.000001869164408: 0.0002004409701343
  mean = 0.33658

  var : min piece
  Probabilities (truncated):
  0.328089852504474: 0.0002004409701343
  0.325504220242381: 0.0002004409701343
  0.323122162988001: 0.0002004409701343
  0.322687707525998: 0.0002004409701343
  .........
  0.000090186949861: 0.0002004409701343
  0.000076612457669: 0.0002004409701343
  0.000036415644007: 0.0002004409701343
  0.000001869164408: 0.0002004409701343
  mean = 0.110435

  var : p
  Probabilities:
  true: 0.8410503106835037
  false: 0.1589496893164963
  mean = [true = 0.84105,false = 0.15895]

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>
  [B1,B2] = uniform_n(0,1,2),

  Piece1 = min(B1,B2),
  Piece2 = max(B1,B2),
  Piece3 = 1 - Piece1 - Piece2,
  
  observe(Piece3 > 0),

  MinPiece = min([Piece1,Piece2,Piece3]),

  P = check(MinPiece < 1/5),
  
  if observed_ok then
    add_all([ ["b1",B1],
              ["b2",B2],
              ["piece1",Piece1],
              ["piece2",Piece2],
              ["piece3",Piece3],
              ["min piece",MinPiece],              
              ["p",P]])
  end.


/*
  Simpler model: using Dirichlet

  var : min piece
  Probabilities (truncated):
  0.329230016240484: 0.0001000000000000
  0.328566213962379: 0.0001000000000000
  0.326937098046785: 0.0001000000000000
  0.325474263974374: 0.0001000000000000
  .........
  0.000103532086083: 0.0001000000000000
  0.000087693154419: 0.0001000000000000
  0.000080913799495: 0.0001000000000000
  0.000037914044021: 0.0001000000000000
  mean = 0.111541

  var : p
  Probabilities:
  true: 0.8390000000000000
  false: 0.1610000000000000
  mean = [true = 0.839,false = 0.161]

*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go2 => true.
  
  
model2() =>
  Sticks = dirichlet_dist([1,1,1]),
  MinPiece = min(Sticks),
  P = check(MinPiece < 1/5),
  add("min piece",MinPiece),
  add("p",P).