/* 

  My neighbour in Picat.

  This is a "real world" problem (happening a long time ago).

  My beautiful neighbour and I look out of our windows M and N times a day, respectively. 
  How many times do we look out at the same time?

  Say that there are 24*60 (=1440) timeslots a day, and we each look out 40 times a day
  (independently of each other). Then, on average, we would look out at the same time
  about once a day.


  If the probability of looking out are the same this is a matching problem. 
  For a simpler problem, let's say that we reduce this to looking out once in an hour.
  Picat> pdf_all($matching_dist(24),0.0001,0.9999999).print_list
  0 = 0.367879
  1 = 0.367879
  2 = 0.18394
  3 = 0.0613132
  4 = 0.0153283
  5 = 0.00306566
  6 = 0.000510944
  7 = 7.2992e-05
  8 = 9.12399e-06
  9 = 1.01378e-06
  10 = 1.01378e-07

  Picat> X=1/exp(1)
  X = 0.367879441171442


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

/*

  With some different lookout times:

  [m = 40,n = 40]
  var : same time
  Probabilities:
  1: 0.3566666666666667
  0: 0.3373333333333333
  2: 0.2073333333333333
  3: 0.0723333333333333
  4: 0.0206666666666667
  5: 0.0043333333333333
  6: 0.0013333333333333
  mean = 1.10067


  [m = 100,n = 40]
  var : same time
  Probabilities:
  2: 0.2483333333333333
  3: 0.2270000000000000
  1: 0.1716666666666667
  4: 0.1466666666666667
  5: 0.0820000000000000
  0: 0.0593333333333333
  6: 0.0376666666666667
  7: 0.0183333333333333
  8: 0.0056666666666667
  9: 0.0023333333333333
  10: 0.0010000000000000
  mean = 2.77667


  [m = 200,n = 40]
  var : same time
  Probabilities:
  5: 0.1813333333333333
  6: 0.1536666666666667
  4: 0.1436666666666667
  7: 0.1276666666666667
  3: 0.1096666666666667
  8: 0.0876666666666667
  9: 0.0576666666666667
  2: 0.0563333333333333
  10: 0.0310000000000000
  1: 0.0216666666666667
  11: 0.0133333333333333
  12: 0.0073333333333333
  0: 0.0040000000000000
  13: 0.0020000000000000
  14: 0.0016666666666667
  15: 0.0013333333333333
  mean = 5.59467

*/
go ?=>
  member([M,N],[[40,40],[100,40],[200,40]]),
  println([m=M,n=N]),
  reset_store,
  run_model(3_000,$model(M,N),[show_probs,mean % ,
                           % show_percentiles,show_histogram,
                           % show_hpd_intervals,hpd_intervals=[0.94],
                           % min_accepted_samples=1000,show_accepted_samples=true
                          ]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.
  
  
model(M,N) =>
  Mins = 24*60, % Number of minutes per day
  % M = 40, % Number of times I look out
  % N = 40, % Number of time my neighbour look out

  % The time slots that we look out
  Me        = [ bern(N/Mins) : _ in 1..Mins],
  Neighbour = [ bern(M/Mins) : _ in 1..Mins],  

  % Total number of looking out
  % SMe = Me.sum,
  % SNeighbour = Neighbour.sum,

  % Do we look out at the same time?
  SameTime = [cond( (Me[I] == 1, Neighbour[I] == 1), 1,0) : I in 1..Mins].sum,

  add("same time",SameTime).
  % add("s me",SMe),
  % add("s neighbour",SNeighbour).