/* 

  Fair coin from a biased coin in Picat.

  http://cplint.eu/e/von_neumann_trick.swinb
  """
  If you have a biased coin, how to obtain a fair coin, i.e., from a coin that lands heads with 
  probability p with p≠0.5, how to generate samples from (heads,tails) with P(heads)=P(tails)=0.5?

  John von Neuamnn gave the following solution

   1) Toss the coin twice.
   2) If you get HT return H, if you get TH return T
   3) Otherwise, return to 1 and start over

  In fact, if p is the probability of the biased coin landing heads, then the outcomes HT 
  and TH are equiprobable with probability p(1−p). However, with probability pp+(1−p)(1−p)
  these results are not obtained. In this case, we simply repeat the process. The probability 
  that we never get HT or TH is 0 so this procedure ends with certainty.

  See   
  - https://en.wikipedia.org/wiki/Fair_coin#Fair_results_from_a_biased_coin
  - von Neumann, John (1951). "Various techniques used in connection with random digits". 
    National Bureau of Standards Applied Math Series. 12: 36.
  """

  Cf 
  - ppl_unbiased_coin.pi
  - my Gamble model gamble_fair_coin_from_a_biased_coin.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.

/*

  var : biased coin 1
  Probabilities:
  head: 0.7946000000000000
  tail: 0.2054000000000000
  mean = [head = 0.7946,tail = 0.2054]

  var : coin 1
  Probabilities:
  head: 0.5004000000000000
  tail: 0.4996000000000000
  mean = [head = 0.5004,tail = 0.4996]

*/
go ?=>
  reset_store,
  run_model(10_000,$model,[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go => true.
  
biased_coin(C) = categorical([0.8,0.2],[head,tail]).
coin(C) = cond(BiasedCoinC != biased_coin(C+1),
               BiasedCoinC,
               coin(C+2)) =>
              BiasedCoinC = biased_coin(C).
model() =>
  Coin1 = coin(1),
  BiasedCoin1 = biased_coin(1),
  add("coin 1",Coin1),
  add("biased coin 1",BiasedCoin1).  


/*  
  Just for fun: Implementing von Neumann's method above:

  1) Toss the coin twice.
  2) If you get HT return H, if you get TH return T
  3) Otherwise, return to 1 and start over

  var : coin
  Probabilities:
  head: 0.5011000000000000
  tail: 0.4989000000000000
  mean = [head = 0.5011,tail = 0.4989]

*/
go2 ?=>
  reset_store,
  run_model(10_000,$model2,[show_probs_trunc,mean]),
  nl,
  % show_store_lengths,
  % fail,
  nl.
go2 => true.
  
model2() =>
  % 1. Toss the coin twice
  Cs = categorical_n([0.8,0.2],[head,tail],2),

  OK = false,
  Coin = _,
  while(OK == false)
    % if you get HT return H, if you get TH return T
    if Cs[1] != Cs[2] then
      Coin = Cs[1],
      OK := true
    else
      % otherwise, return to 1
      Cs := categorical_n([0.8,0.2],[head,tail],2),
    end
  end,

  add("coin",Coin).