/* 

  How many cards can you draw? in Picat.

  From Dr. Robert Kübler: "How Many Cards Can You Draw?"
  https://medium.com/science-spectrum/how-many-cards-can-you-draw-74fba0e679d3
  """
  Assume that you have a deck of n cards, numbered from 1 to n. You shuffle 
  it and turn it face down. Now, you keep drawing cards until you get 
  one that is smaller than the one you had before. What is the expected 
  number of cards you will draw from the pile?
  """

  For n = 10: The expected value is 2.7182818011463845,
  and - interestingly - that's E (as stated in the article).

  Maple:
  > identify(2.7182815);
                                                 exp(1)

 

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

/*
n = 10
var : a
Probabilities:
2: 0.4976000000000000
3: 0.3334800000000000
4: 0.1273200000000000
5: 0.0337600000000000
6: 0.0065400000000000
7: 0.0011400000000000
8: 0.0001200000000000
9: 0.0000400000000000
mean = 2.72226


n = 20
var : a
Probabilities:
2: 0.5006400000000000
3: 0.3320400000000000
4: 0.1249200000000000
5: 0.0337000000000000
6: 0.0071200000000000
7: 0.0014400000000000
9: 0.0000600000000000
8: 0.0000600000000000
10: 0.0000200000000000
mean = 2.7196


n = 100
var : a
Probabilities:
2: 0.4978200000000000
3: 0.3374400000000000
4: 0.1245600000000000
5: 0.0321800000000000
6: 0.0066200000000000
7: 0.0011600000000000
8: 0.0002000000000000
9: 0.0000200000000000
mean = 2.71672



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

f(C,Cs,N) = Res =>
  if Cs == [] then
    Res = N
  else
    T = Cs.first,
    if T < C then
      Res = N + 1
    else
      Res = f(Cs.first,Cs.rest,N+1)
    end
  end.


model(N) =>
  Cards = draw_without_replacement(N,1..N),
  % Cards = shuffle(1..N),
  A = f(Cards.first,Cards.rest,1),
  add("a",A).