/* 

  Chinese restraurant process in Picat.

  https://en.wikipedia.org/wiki/Chinese_restaurant_process
  """
  In probability theory, the Chinese restaurant process is a discrete-time stochastic process, 
  analogous to seating customers at tables in a restaurant. Imagine a restaurant with an 
  infinite number of circular tables, each with infinite capacity. Customer 1 sits at the 
  first table. The next customer either sits at the same table as customer 1, or the next 
  table. This continues, with each customer choosing to either sit at an occupied table with a 
  probability proportional to the number of customers already there (i.e., they are more likely 
  to sit at a table with many customers than few), or an unoccupied table. At time n, the 
  n customers have been partitioned among m ≤ n tables (or blocks of the partition). The results 
  of this process are exchangeable, meaning the order in which the customers sit does not affect 
  the probability of the final distribution. This property greatly simplifies a number of problems 
  in population genetics, linguistic analysis, and image recognition.

  The restaurant analogy first appeared in a 1985 write-up by David Aldous,(1) where it was 
  attributed to Jim Pitman (who additionally credits Lester Dubins).(2)

  ...

  At time n+1 the element "n+1" is either
  1. added to one of the blocks of the partition Bn, where each block is chosen with probability
     |b|/(n+1) where |b| is the size of the block (i.e. the number of elements), or
  2. added to the partition Bn as new singleton block, with probability 1/(n+1)

  ...

  The probability assigned to any particular partition (ignoring the order in which customers 
  sit around any particular table) is

   Pr(Bn = B) =   prod (b in B) * (|b|-1)!
                  ------------------------   B in Pn
                        n! 

  where b is a block in the partition B and |b| is the size of b.
  [Pn is the partitions of n]

  """

  Cf 
  - ppl_chinese_restraurant_process2.pi 
    with a different model as well as showing 
    the probabillity distribution crp_dist.

  - Gamble model gamble_chinese_restaurant_process.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.

/*
  Here are simulation for some different values of max guests

  max_quests = 10
  var : a
  Probabilities (truncated):
  [8,1,1]: 0.0552000000000000
  [7,1,1,1]: 0.0511000000000000
  [7,2,1]: 0.0374000000000000
  [6,2,1,1]: 0.0331000000000000
  [9,1]: 0.0283000000000000
  [6,3,1]: 0.0265000000000000
  [5,3,1,1]: 0.0241000000000000
  [6,1,2,1]: 0.0208000000000000
  [6,1,1,1,1]: 0.0192000000000000
  [7,1,2]: 0.0188000000000000
  .........
  [1,1,2,2,1,3]: 0.0001000000000000
  [1,1,2,2,1,1,1,1]: 0.0001000000000000
  [1,1,2,1,5]: 0.0001000000000000
  [1,1,2,1,3,2]: 0.0001000000000000
  [1,1,2,1,2,2,1]: 0.0001000000000000
  [1,1,1,4,1,2]: 0.0001000000000000
  [1,1,1,3,3,1]: 0.0001000000000000
  [1,1,1,3,1,1,2]: 0.0001000000000000
  [1,1,1,1,3,2,1]: 0.0001000000000000
  [1,1,1,1,2,3,1]: 0.0001000000000000
  mean = []

  var : a[1]
  Probabilities:
  5: 0.1414000000000000
  4: 0.1395000000000000
  3: 0.1363000000000000
  6: 0.1336000000000000
  2: 0.1229000000000000
  7: 0.1172000000000000
  1: 0.1098000000000000
  8: 0.0686000000000000
  9: 0.0283000000000000
  10: 0.0024000000000000
  mean = 4.479

  var : max a
  Probabilities:
  5: 0.2371000000000000
  4: 0.2128000000000000
  6: 0.1917000000000000
  7: 0.1472000000000000
  3: 0.0925000000000000
  8: 0.0804000000000000
  9: 0.0308000000000000
  2: 0.0051000000000000
  10: 0.0024000000000000
  mean = 5.4494

  var : len a
  Probabilities:
  4: 0.3617000000000000
  3: 0.2989000000000000
  5: 0.1992000000000000
  2: 0.0746000000000000
  6: 0.0565000000000000
  7: 0.0064000000000000
  1: 0.0024000000000000
  8: 0.0003000000000000
  mean = 3.8773


  max_quests = 20
  var : a
  Probabilities (truncated):
  [16,1,1,1,1]: 0.0070000000000000
  [17,1,1,1]: 0.0059000000000000
  [15,2,1,1,1]: 0.0051000000000000
  [13,4,1,1,1]: 0.0051000000000000
  [16,2,1,1]: 0.0050000000000000
  [15,1,1,1,1,1]: 0.0044000000000000
  [15,3,1,1]: 0.0035000000000000
  [16,1,2,1]: 0.0034000000000000
  [12,5,1,1,1]: 0.0032000000000000
  [10,7,1,1,1]: 0.0032000000000000
  .........
  [1,1,5,4,3,3,3]: 0.0001000000000000
  [1,1,4,13,1]: 0.0001000000000000
  [1,1,4,1,2,8,1,2]: 0.0001000000000000
  [1,1,3,7,1,4,2,1]: 0.0001000000000000
  [1,1,2,14,2]: 0.0001000000000000
  [1,1,2,13,2,1]: 0.0001000000000000
  [1,1,2,10,3,3]: 0.0001000000000000
  [1,1,2,8,4,1,2,1]: 0.0001000000000000
  [1,1,1,5,3,4,3,1,1]: 0.0001000000000000
  [1,1,1,2,2,10,2,1]: 0.0001000000000000
  mean = []

  var : a[1]
  Probabilities:
  8: 0.0811000000000000
  10: 0.0759000000000000
  7: 0.0758000000000000
  6: 0.0724000000000000
  5: 0.0722000000000000
  11: 0.0702000000000000
  9: 0.0701000000000000
  4: 0.0700000000000000
  2: 0.0682000000000000
  3: 0.0664000000000000
  12: 0.0623000000000000
  13: 0.0569000000000000
  1: 0.0563000000000000
  14: 0.0400000000000000
  15: 0.0315000000000000
  16: 0.0193000000000000
  17: 0.0082000000000000
  18: 0.0029000000000000
  19: 0.0003000000000000
  mean = 7.8347

  var : max a
  Probabilities:
  8: 0.1377000000000000
  9: 0.1312000000000000
  7: 0.1162000000000000
  10: 0.1157000000000000
  11: 0.1038000000000000
  12: 0.0856000000000000
  6: 0.0777000000000000
  13: 0.0721000000000000
  14: 0.0498000000000000
  15: 0.0379000000000000
  5: 0.0311000000000000
  16: 0.0216000000000000
  17: 0.0094000000000000
  4: 0.0065000000000000
  18: 0.0031000000000000
  19: 0.0004000000000000
  3: 0.0002000000000000
  mean = 9.8419

  var : len a
  Probabilities:
  6: 0.2910000000000000
  5: 0.2781000000000000
  7: 0.1675000000000000
  4: 0.1418000000000000
  8: 0.0698000000000000
  3: 0.0289000000000000
  9: 0.0177000000000000
  10: 0.0032000000000000
  2: 0.0015000000000000
  11: 0.0005000000000000
  mean = 5.7211

  max_quests = 50
  var : a
  Probabilities (truncated):
  [40,3,2,1,1,1,1,1]: 0.0003000000000000
  [43,1,1,1,2,1,1]: 0.0002000000000000
  [40,2,2,1,3,1,1]: 0.0002000000000000
  [40,1,2,2,1,2,1,1]: 0.0002000000000000
  [40,1,2,1,1,3,1,1]: 0.0002000000000000
  [39,3,2,2,1,1,1,1]: 0.0002000000000000
  [39,2,2,3,1,1,1,1]: 0.0002000000000000
  [38,6,1,1,1,1,1,1]: 0.0002000000000000
  [38,2,4,2,1,1,1,1]: 0.0002000000000000
  [38,1,2,1,3,1,1,2,1]: 0.0002000000000000
  .........
  [1,2,12,20,1,8,2,2,1,1]: 0.0001000000000000
  [1,2,7,17,2,8,4,1,3,2,3]: 0.0001000000000000
  [1,2,3,21,16,2,1,1,1,1,1]: 0.0001000000000000
  [1,1,28,2,7,6,3,2]: 0.0001000000000000
  [1,1,24,10,8,2,1,1,1,1]: 0.0001000000000000
  [1,1,24,3,10,5,1,3,1,1]: 0.0001000000000000
  [1,1,17,8,12,3,1,1,2,3,1]: 0.0001000000000000
  [1,1,15,8,9,7,3,4,1,1]: 0.0001000000000000
  [1,1,2,26,3,8,1,3,1,1,1,2]: 0.0001000000000000
  [1,1,1,1,4,6,15,14,1,1,4,1]: 0.0001000000000000
  mean = []

  var : a[1]
  Probabilities (truncated):
  13: 0.0354000000000000
  16: 0.0350000000000000
  23: 0.0334000000000000
  15: 0.0334000000000000
  18: 0.0330000000000000
  11: 0.0329000000000000
  9: 0.0322000000000000
  12: 0.0318000000000000
  14: 0.0313000000000000
  21: 0.0310000000000000
  .........
  37: 0.0069000000000000
  38: 0.0064000000000000
  39: 0.0041000000000000
  40: 0.0030000000000000
  41: 0.0015000000000000
  43: 0.0009000000000000
  42: 0.0009000000000000
  44: 0.0003000000000000
  45: 0.0002000000000000
  46: 0.0001000000000000
  mean = 17.2837

  var : max a
  Probabilities (truncated):
  19: 0.0589000000000000
  18: 0.0584000000000000
  16: 0.0571000000000000
  21: 0.0555000000000000
  17: 0.0546000000000000
  23: 0.0541000000000000
  20: 0.0529000000000000
  22: 0.0524000000000000
  15: 0.0469000000000000
  24: 0.0463000000000000
  .........
  40: 0.0037000000000000
  9: 0.0019000000000000
  41: 0.0017000000000000
  43: 0.0010000000000000
  42: 0.0009000000000000
  8: 0.0007000000000000
  44: 0.0003000000000000
  45: 0.0002000000000000
  7: 0.0002000000000000
  46: 0.0001000000000000
  mean = 22.3739

  var : len a
  Probabilities:
  9: 0.2331000000000000
  10: 0.2101000000000000
  8: 0.1775000000000000
  11: 0.1371000000000000
  7: 0.0950000000000000
  12: 0.0684000000000000
  6: 0.0301000000000000
  13: 0.0293000000000000
  14: 0.0097000000000000
  5: 0.0066000000000000
  15: 0.0024000000000000
  4: 0.0005000000000000
  16: 0.0002000000000000
  mean = 9.3843

  max_quests = 100
  var : a
  Probabilities (truncated):
  [86,2,1,1,1,1,3,1,1,1,1,1]: 0.0001000000000000
  [85,3,3,1,1,1,2,2,1,1]: 0.0001000000000000
  [85,1,1,1,2,1,1,4,2,1,1]: 0.0001000000000000
  [84,3,6,1,1,1,1,1,1,1]: 0.0001000000000000
  [84,2,1,1,8,1,1,2]: 0.0001000000000000
  [84,1,4,3,1,1,2,1,1,1,1]: 0.0001000000000000
  [84,1,1,1,1,7,1,1,1,1,1]: 0.0001000000000000
  [83,4,1,2,2,1,1,2,2,1,1]: 0.0001000000000000
  [83,2,1,1,2,1,3,3,1,1,1,1]: 0.0001000000000000
  [83,1,2,3,4,2,1,1,1,1,1]: 0.0001000000000000
  .........
  [1,4,3,46,16,13,3,1,5,3,1,1,1,1,1]: 0.0001000000000000
  [1,3,38,12,2,5,9,4,8,2,3,4,2,1,1,2,2,1]: 0.0001000000000000
  [1,3,2,12,44,13,5,1,1,2,7,2,1,2,1,1,1,1]: 0.0001000000000000
  [1,2,61,5,2,11,1,12,1,1,1,1,1]: 0.0001000000000000
  [1,2,13,22,18,13,5,3,13,1,2,2,2,2,1]: 0.0001000000000000
  [1,2,9,53,19,1,2,2,2,1,3,1,2,1,1]: 0.0001000000000000
  [1,1,66,3,2,6,1,2,5,6,2,1,1,1,2]: 0.0001000000000000
  [1,1,58,1,8,10,2,3,4,2,2,1,1,1,1,2,1,1]: 0.0001000000000000
  [1,1,29,25,16,17,2,1,1,1,2,1,1,1,1]: 0.0001000000000000
  [1,1,22,25,37,4,6,1,1,1,1]: 0.0001000000000000
  mean = []

  var : a[1]
  Probabilities (truncated):
  28: 0.0193000000000000
  25: 0.0190000000000000
  27: 0.0184000000000000
  29: 0.0183000000000000
  38: 0.0181000000000000
  18: 0.0181000000000000
  11: 0.0180000000000000
  36: 0.0178000000000000
  1: 0.0178000000000000
  19: 0.0177000000000000
  .........
  76: 0.0020000000000000
  78: 0.0015000000000000
  79: 0.0007000000000000
  82: 0.0006000000000000
  84: 0.0004000000000000
  81: 0.0004000000000000
  80: 0.0004000000000000
  83: 0.0003000000000000
  85: 0.0002000000000000
  86: 0.0001000000000000
  mean = 32.5847

  var : max a
  Probabilities (truncated):
  34: 0.0337000000000000
  38: 0.0319000000000000
  40: 0.0310000000000000
  36: 0.0306000000000000
  42: 0.0300000000000000
  39: 0.0292000000000000
  41: 0.0291000000000000
  37: 0.0291000000000000
  35: 0.0285000000000000
  33: 0.0277000000000000
  .........
  80: 0.0005000000000000
  84: 0.0004000000000000
  83: 0.0004000000000000
  81: 0.0004000000000000
  16: 0.0004000000000000
  85: 0.0002000000000000
  14: 0.0002000000000000
  86: 0.0001000000000000
  15: 0.0001000000000000
  13: 0.0001000000000000
  mean = 42.7675

  var : len a
  Probabilities:
  13: 0.1881000000000000
  14: 0.1771000000000000
  12: 0.1567000000000000
  15: 0.1354000000000000
  11: 0.1023000000000000
  16: 0.0916000000000000
  10: 0.0501000000000000
  17: 0.0484000000000000
  18: 0.0204000000000000
  9: 0.0132000000000000
  19: 0.0088000000000000
  20: 0.0036000000000000
  8: 0.0028000000000000
  21: 0.0007000000000000
  7: 0.0004000000000000
  24: 0.0002000000000000
  23: 0.0001000000000000
  22: 0.0001000000000000
  mean = 13.5252

*/
go ?=>
  member(MaxGuests,[10,20,50,100]),
  MaxGuests = 100,
  println(max_quests=MaxGuests),
  reset_store,
  run_model(10_000,$model(MaxGuests),[show_probs_trunc,mean,truncate_size=10]),
  nl,
  fail,
  nl.
go => true.
  

% Generate the table of the next customer
% (or return a if all customers are placed).
f(A,MaxGuests) = Res =>
  Len = A.len,
  if A.sum >= MaxGuests then
    Res = A
  else
    if flip(1/(Len+1)) == true then
      % Pick a new empty table
      Res = f(A ++ [1],MaxGuests)      
    else
      % Pick a table proportion to the people already sitting there
      T = categorical(A,1..Len),
      NewA = [cond(I == T,A[I]+1,A[I]) : I in 1..Len],
      Res = f(NewA,MaxGuests)
    end
  end.
  
model(MaxGuests) =>
  A = f([1],MaxGuests),
  
  add("a",A),
  add("a[1]",A[1]),
  add("max a",A.max),
  add("len a",A.len).