/* 

  Messing with envelopes in Picat.

  https://brainstellar.com/puzzles/probability/205
  """
  There are n letters and n envelopes. You put the letters randomly in the 
  envelopes so that each letter is in one envelope. (Effectively a random permutation of 
  n numbers chosen uniformly). Calculate the expected number of envelopes with the correct 
  letter inside them.

  Answer: 1
  """

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

/*
  n = 5
  var : p
  Probabilities:
  1: 0.3752000000000000
  0: 0.3638000000000000
  2: 0.1627000000000000
  3: 0.0897000000000000
  5: 0.0086000000000000
  mean = 1.0127

  n = 10
  var : p
  Probabilities:
  0: 0.3688000000000000
  1: 0.3620000000000000
  2: 0.1849000000000000
  3: 0.0643000000000000
  4: 0.0156000000000000
  5: 0.0036000000000000
  6: 0.0007000000000000
  8: 0.0001000000000000
  mean = 1.0101

  n = 100
  var : p
  Probabilities:
  0: 0.3717000000000000
  1: 0.3673000000000000
  2: 0.1806000000000000
  3: 0.0633000000000000
  4: 0.0138000000000000
  5: 0.0026000000000000
  6: 0.0006000000000000
  7: 0.0001000000000000
  mean = 0.9909

  n = 500
  var : p
  Probabilities:
  0: 0.3710000000000000
  1: 0.3608000000000000
  2: 0.1836000000000000
  3: 0.0660000000000000
  4: 0.0146000000000000
  5: 0.0035000000000000
  6: 0.0005000000000000
  mean = 1.0049

  n = 1000
  var : p
  Probabilities:
  1: 0.3687000000000000
  0: 0.3633000000000000
  2: 0.1859000000000000
  3: 0.0625000000000000
  4: 0.0159000000000000
  5: 0.0027000000000000
  6: 0.0009000000000000
  7: 0.0001000000000000
  mean = 1.0112


  Compare with Poisson(1):

  Picat> pdf_all($poisson_dist(1),0.00001,0.99999).printf_list
  0 0.367879441171443 
  1 0.367879441171442 
  2 0.183939720585721 
  3 0.06131324019524 
  4 0.01532831004881 
  5 0.003065662009762 
  6 0.000510943668294 
  7 0.000072991952613 
  8 0.000009123994077 

  And with 1/exp(1):
  Picat> X=1/exp(1)  
  X = 0.367879441171442


*/
go ?=>
  member(N,[5,10,100,500,1000]),
  println(n=N),
  reset_store,
  run_model(10_000,$model(N),[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,nl,
  fail,
  nl.
go => true.
  
  
model(N) =>
  
  % Letter = draw_without_replacement(N,1..N),
  Letter = shuffle(1..N),
  P = [cond(Letter[I] == I,1,0) : I in 1..N].sum, 
  add("p",P).