/* 

  Will the witches meet? in Picat.

  From Pascal Bercker: "Will The Witches Meet? A Reply To BL":
  https://medium.com/@pbercker/will-the-witches-meet-a-reply-to-bl-af79bf3fac96
  """
  Fellow Medium writer and puzzle enthusiast BL poses the following (simple) puzzle. 
  Two witches meet between midnight and 1:00 am, but each chooses their arrival time randomly.
  '''
  If both witches choose their arrival time at random between midnight and one, and
  each remains for precisely half an hour, what it the probability that their
  midnight paths cross and they meet inside the cofeehouse.
  '''
  """

  As Pascal mentions, the assumptions/priors when reading the problem can make
  difference of the answer of the probability the the witches will meet.

  Here we test:
  - uniform distribution (0..59)
  - random integers 0..59
  - triangular distribution(0,59)
  - normal distribution (normal(30,10) truncated/observed to be between >= 0 and < 60).
  - mean time the witches meet
  - more than 2 witches


  BL's Medium post: "How Often Do The Two Witches Meet At The Coffee Shop?"
  https://medium.com/math-games/how-often-do-the-two-witches-meet-at-the-coffee-shop-bed35bc8e33c


  Cf some other meeting problems:
  - ppl_meeting_collegues_at_office.pi
  - ppl_meeting_problem2.pi
  - ppl_meeting_problem.pi
  - ppl_meeting_under_the_clock.pi
  - ppl_meeting_under_the_clock2.pi


  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, ppl_common_utils.
import util.
% import ordset.

main => go.

/*

  var : w1 start
  Probabilities (truncated):
  59.999607196077527: 0.0001000000000000
  59.999123588203979: 0.0001000000000000
  59.993544258174275: 0.0001000000000000
  59.979135850434716: 0.0001000000000000
  .........
  0.024774279456946: 0.0001000000000000
  0.020866021523655: 0.0001000000000000
  0.009452272210946: 0.0001000000000000
  0.005394900220165: 0.0001000000000000
  mean = 29.9495
  HPD intervals:
  HPD interval (0.84): 8.40881749447846..58.55149414322874
  HPD interval (0.94): 2.08116981297786..58.10830196277625
  HPD interval (0.99): 0.02086602152366..59.31782168304446

  var : w1 end
  Probabilities (truncated):
  60: 0.4976000000000000
  59.996292530557277: 0.0001000000000000
  59.995601442640464: 0.0001000000000000
  59.992720219303258: 0.0001000000000000
  .........
  30.024774279456945: 0.0001000000000000
  30.020866021523656: 0.0001000000000000
  30.009452272210947: 0.0001000000000000
  30.005394900220164: 0.0001000000000000
  mean = 52.4638
  HPD intervals:
  HPD interval (0.84): 39.56160455456078..60.00000000000000
  HPD interval (0.94): 33.69628180921836..60.00000000000000
  HPD interval (0.99): 30.58513328460284..60.00000000000000

  var : w2 start
  Probabilities (truncated):
  59.999170498922084: 0.0001000000000000
  59.994706716385998: 0.0001000000000000
  59.980738414442513: 0.0001000000000000
  59.953915737547867: 0.0001000000000000
  .........
  0.015431195504699: 0.0001000000000000
  0.014497619129018: 0.0001000000000000
  0.011430764576155: 0.0001000000000000
  0.009953062986002: 0.0001000000000000
  mean = 29.9154
  HPD intervals:
  HPD interval (0.84): 6.18846380439981..56.56482894744948
  HPD interval (0.94): 3.16319419218376..59.51294266596108
  HPD interval (0.99): 0.00995306298600..59.36672248848096

  var : w2 end
  Probabilities (truncated):
  60: 0.4994000000000000
  59.998586154542203: 0.0001000000000000
  59.994228403081294: 0.0001000000000000
  59.991738558743023: 0.0001000000000000
  .........
  30.0154311955047: 0.0001000000000000
  30.014497619129017: 0.0001000000000000
  30.011430764576154: 0.0001000000000000
  30.009953062986003: 0.0001000000000000
  mean = 52.3952
  HPD intervals:
  HPD interval (0.84): 39.43286438911821..60.00000000000000
  HPD interval (0.94): 33.56600328514632..60.00000000000000
  HPD interval (0.99): 30.52181527042846..60.00000000000000

  var : p overlap
  Probabilities:
  true: 0.7429000000000000
  false: 0.2571000000000000
  mean = 0.7429
  HPD intervals:
  []

  var : diff
  Probabilities (truncated):
  59.636819651181263: 0.0001000000000000
  59.299453021632253: 0.0001000000000000
  59.153872271605707: 0.0001000000000000
  59.006508215799229: 0.0001000000000000
  .........
  0.009377086539466: 0.0001000000000000
  0.008901217956513: 0.0001000000000000
  0.008667698134051: 0.0001000000000000
  0.007027359682617: 0.0001000000000000
  mean = 20.1099
  HPD intervals:
  HPD interval (0.84): 0.00702735968262..36.22522829856035
  HPD interval (0.94): 0.00702735968262..45.52049100656085
  HPD interval (0.99): 0.00702735968262..53.77703650564748

  var : p diff
  Probabilities:
  true: 0.7429000000000000
  false: 0.2571000000000000
  mean = 0.7429


*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean,
                           show_hpd_intervals,hpd_intervals=[0.84,0.94,0.99]
                          ]),
  nl,
  nl.
go => true.
  
  
model() =>
  W1Start = uniform(0,60),
  W2Start = uniform(0,60),  
  
  W1End = min(W1Start+30,60),
  W2End = min(W2Start+30,60),

  Diff = abs(W1Start-W2Start),

  POverlap = check ( overlap(W1Start,W1End, W2Start,W2End)),
  PDiff = check(Diff <= 30),

  add("w1 start",W1Start),
  add("w1 end",W1End),
  add("w2 start",W2Start),
  add("w2 end",W2End),
  add("p overlap",POverlap),
  add("diff",Diff),
  add("p diff",PDiff).


/*
  Let's assume that the witches - for some strange reasons - only arrive 
  only at exact minutes.

  var : w1
  Probabilities (truncated):
  24: 0.0202000000000000
  19: 0.0201000000000000
  9: 0.0198000000000000
  35: 0.0196000000000000
  .........
  14: 0.0148000000000000
  10: 0.0147000000000000
  17: 0.0144000000000000
  2: 0.0144000000000000
  mean = 29.2758
  HPD intervals:
  HPD interval (0.84): 0.00000000000000..49.00000000000000
  HPD interval (0.94): 0.00000000000000..55.00000000000000
  HPD interval (0.99): 0.00000000000000..58.00000000000000

  var : w2
  Probabilities (truncated):
  7: 0.0201000000000000
  5: 0.0199000000000000
  12: 0.0198000000000000
  41: 0.0197000000000000
  .........
  44: 0.0152000000000000
  9: 0.0151000000000000
  28: 0.0145000000000000
  2: 0.0136000000000000
  mean = 28.9372
  HPD intervals:
  HPD interval (0.84): 0.00000000000000..49.00000000000000
  HPD interval (0.94): 0.00000000000000..55.00000000000000
  HPD interval (0.99): 0.00000000000000..58.00000000000000

  var : diff
  Probabilities (truncated):
  3: 0.0331000000000000
  5: 0.0328000000000000
  7: 0.0323000000000000
  1: 0.0314000000000000
  .........
  55: 0.0017000000000000
  56: 0.0016000000000000
  57: 0.0010000000000000
  58: 0.0007000000000000
  mean = 19.531
  HPD intervals:
  HPD interval (0.84): 0.00000000000000..35.00000000000000
  HPD interval (0.94): 0.00000000000000..44.00000000000000
  HPD interval (0.99): 0.00000000000000..53.00000000000000

  var : p
  Probabilities:
  true: 0.7694000000000000
  false: 0.2306000000000000
  mean = 0.7694



*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[show_probs_trunc,mean,
                           show_hpd_intervals,hpd_intervals=[0.84,0.94,0.99],
                           show_histogram, histogram_limit=60
                          ]),
  nl,
  nl.
go2 => true.
  
  
model2() =>
  W1 = random_integer(59),
  W2 = random_integer(59),  

  Diff = abs(W1-W2),
  P = check(Diff <= 30),

  add("w1",W1),
  add("w2",W2),  
  add("diff",Diff),
  add("p",P).



/*
  Pascal also tested triangular distribution (with some support of
  the validity of that from both Shakespeare and ChatGPT). 
  Let's do that as wells since Picat PPL happens to support this
  distribution.

  var : w1
  Probabilities (truncated):
  59.757739631729628: 0.0001000000000000
  59.176925132974219: 0.0001000000000000
  59.080491444111622: 0.0001000000000000
  58.968823137252521: 0.0001000000000000
  .........
  1.228343968594443: 0.0001000000000000
  1.190327093347633: 0.0001000000000000
  1.10397641537429: 0.0001000000000000
  1.044192951431626: 0.0001000000000000
  mean = 29.9973
  HPD intervals:
  HPD interval (0.84): 11.58852075637546..47.43170207753045
  HPD interval (0.94): 8.09340593561943..53.22836774252761
  HPD interval (0.99): 3.70493340484392..57.67531778190125

  var : w2
  Probabilities (truncated):
  59.661230139655146: 0.0001000000000000
  59.606216799817943: 0.0001000000000000
  59.211980347241713: 0.0001000000000000
  59.174017745138663: 0.0001000000000000
  .........
  0.759828680290664: 0.0001000000000000
  0.598296243679302: 0.0001000000000000
  0.520274945906542: 0.0001000000000000
  0.100533073408448: 0.0001000000000000
  mean = 29.9369
  HPD intervals:
  HPD interval (0.84): 11.30649459722487..46.94851037199403
  HPD interval (0.94): 7.59689215228351..52.69281676499671
  HPD interval (0.99): 2.51246526133445..56.49449876179284

  var : diff
  Probabilities (truncated):
  56.94884449258241: 0.0001000000000000
  54.35436164945169: 0.0001000000000000
  54.130177363280112: 0.0001000000000000
  53.097833135877515: 0.0001000000000000
  .........
  0.00473162700446: 0.0001000000000000
  0.003776256717018: 0.0001000000000000
  0.0036453857293: 0.0001000000000000
  0.002631733001778: 0.0001000000000000
  mean = 13.8718
  HPD intervals:
  HPD interval (0.84): 0.00263173300178..24.53538618233492
  HPD interval (0.94): 0.00263173300178..32.19498058122596
  HPD interval (0.99): 0.00263173300178..41.75954018016321

  var : p
  Probabilities:
  true: 0.9186000000000000
  false: 0.0814000000000000
  mean = 0.9186

*/
go3 ?=>
  reset_store,
  run_model(10_000,$model3,[show_probs_trunc,mean,
                           show_hpd_intervals,hpd_intervals=[0.84,0.94,0.99]
                          ]),
  nl,
  nl.
go3 => true.

  
model3() =>
  W1 = triangular_dist(0,60),
  W2 = triangular_dist(0,60),  

  Diff = abs(W1-W2),
  P = check(Diff <= 30),

  add("w1",W1),
  add("w2",W2),  
  add("diff",Diff),
  add("p",P).


/*
  With (truncated) Normal distribution.

  var : w1
  Probabilities (truncated):
  59.701560725659832: 0.0001004318569850
  59.645275515641309: 0.0001004318569850
  59.537999636887491: 0.0001004318569850
  59.195972525020601: 0.0001004318569850
  .........
  0.675203949962942: 0.0001004318569850
  0.517263744771931: 0.0001004318569850
  0.344235450059692: 0.0001004318569850
  0.238997173358158: 0.0001004318569850
  mean = 29.9765
  HPD intervals:
  HPD interval (0.84): 15.97686775449353..43.57432537178706
  HPD interval (0.94): 11.18859921502655..47.84646253316886
  HPD interval (0.99): 4.74286327859867..53.75531566665359

  var : w2
  Probabilities (truncated):
  59.854642307750694: 0.0001004318569850
  59.372763395075218: 0.0001004318569850
  59.315067845037149: 0.0001004318569850
  59.200912563087066: 0.0001004318569850
  .........
  1.097116272679465: 0.0001004318569850
  1.016331411613336: 0.0001004318569850
  0.459625538346778: 0.0001004318569850
  0.008763729079146: 0.0001004318569850
  mean = 30.189
  HPD intervals:
  HPD interval (0.84): 16.09774879239472..43.53447780869855
  HPD interval (0.94): 11.52816604360732..47.99575844020627
  HPD interval (0.99): 5.02212701556905..54.18809802469507

  var : diff
  Probabilities (truncated):
  52.479209223561867: 0.0001004318569850
  47.507566313954854: 0.0001004318569850
  46.846224605425085: 0.0001004318569850
  46.724124132018034: 0.0001004318569850
  .........
  0.00456648637735: 0.0001004318569850
  0.003074797686708: 0.0001004318569850
  0.001633151744983: 0.0001004318569850
  0.001382564691074: 0.0001004318569850
  mean = 11.0742
  HPD intervals:
  HPD interval (0.84): 0.00163315174498..19.34264214287334
  HPD interval (0.94): 0.00138256469107..26.13101475611833
  HPD interval (0.99): 0.00138256469107..35.77451349174540

  var : p
  Probabilities:
  true: 0.9689665561916240
  false: 0.0310334438083760
  mean = 0.968967

*/
  
go4 ?=>
  reset_store,
  run_model(10_000,$model4,[show_probs_trunc,mean,
                           show_hpd_intervals,hpd_intervals=[0.84,0.94,0.99]
                          ]),
  nl,
  nl.
go4 => true.
  

model4() =>
  W1 = normal_dist(30,10),
  W2 = normal_dist(30,10),  

  observe( (W1 >= 0, W1 < 60, W2 >= 0, W2 < 60)),

  Diff = abs(W1-W2),
  P = check(Diff <= 30),

  if observed_ok then
    add("w1",W1),
    add("w2",W2),  
    add("diff",Diff),
    add("p",P)
  end.



/*
  Here is an another model which also checks the average meet time,
  assuming the witches disappear (unmeet) at one o'clock.
  Again the assumption is that they arrive (and depart) at exact minutes.

  var : s
  Probabilities:
  0: 0.2538000000000000
  3: 0.0342000000000000
  1: 0.0341000000000000
  2: 0.0329000000000000
  5: 0.0322000000000000
  7: 0.0309000000000000
  4: 0.0301500000000000
  8: 0.0301000000000000
  6: 0.0296000000000000
  10: 0.0280500000000000
  9: 0.0274500000000000
  11: 0.0271000000000000
  12: 0.0268500000000000
  15: 0.0256000000000000
  18: 0.0248000000000000
  20: 0.0247500000000000
  17: 0.0247500000000000
  13: 0.0245500000000000
  16: 0.0243000000000000
  14: 0.0234500000000000
  19: 0.0230000000000000
  21: 0.0222500000000000
  24: 0.0205000000000000
  22: 0.0205000000000000
  25: 0.0200500000000000
  23: 0.0200500000000000
  27: 0.0196000000000000
  28: 0.0191500000000000
  26: 0.0189500000000000
  29: 0.0186000000000000
  30: 0.0077500000000000
  mean = 10.1887

  var : p
  Probabilities:
  true: 0.7462000000000000
  false: 0.2538000000000000
  mean = 0.7462

*/
go5 ?=>
  reset_store,
  run_model(20_000,$model5,[show_probs_trunc,truncate_size=30,mean]),
  nl,
  nl.
go5 => true.
  

model5() =>
  W1 = random_integer1(60),
  W2 = random_integer1(60),
  % Generate 1s from the arrival + 30 minutes
  W1L = [ cond( (I >= W1, I< W1+30),1,0) : I in 1..60],
  W2L = [ cond( (I>= W2, I< W2+30),1,0) : I in 1..60],

  % When did they meet?
  Same = [cond( (W1L[I] == 1,W2L[I] == 1),1,0) : I in 1..60],
  S = Same.sum,
  P = check(S > 0),
  
  add("s",S),
  add("p",P).


/*
  And let's check for more witches.

  num_witches = 2
  var : s
  mean = 10.1796

  var : p
  mean = 0.7452


  num_witches = 3
  var : s
  mean = 4.8228

  var : p
  mean = 0.4883


  num_witches = 4
  var : s
  mean = 2.2733

  var : p
  mean = 0.2978


  num_witches = 5
  var : s
  mean = 1.1037

  var : p
  mean = 0.176


  num_witches = 6
  var : s
  mean = 0.509

  var : p
  mean = 0.0967


  num_witches = 7
  var : s
  mean = 0.2813

  var : p
  mean = 0.0571


  num_witches = 8
  var : s
  mean = 0.1301

  var : p
  mean = 0.0313


  num_witches = 9
  var : s
  mean = 0.062

  var : p
  mean = 0.0173


  num_witches = 10
  var : s
  mean = 0.0311

  var : p
  mean = 0.0093

*/

go6 ?=>
  member(NumWitches,2..10),
  println(num_witches=NumWitches),
  reset_store,
  run_model(10_000,$model6(NumWitches),[mean]),
  nl,
  fail,
  nl.
go6 => true.
  

model6(NumWitches) =>
  Ws = random_integer1_n(60,NumWitches),
  WL = [ [ cond( (T >= Ws[W], T< Ws[W]+30),1,0) : T in 1..60] :  W in 1..NumWitches],
  
  Same = [ cond( [cond(WL[W,T] == 1,1,0) : W in 1..NumWitches].sum == NumWitches,1,0) : T in 1..60],
  S = Same.sum,
  P = check(S > 0),
  
  add("s",S),
  add("p",P).