/* 

  Rosalind problem Mendel's First Law in Picat.

  https://rosalind.info/problems/iprb/
  """
  Problem

  [Figure 2. The probability of any outcome (leaf) in a probability tree diagram is given by the product 
   of probabilities from the start of the tree to the outcome. For example, the probability that X is 
   blue and Y is blue is equal to (2/5)(1/4), or 1/10.
  ]
  Probability is the mathematical study of randomly occurring phenomena. We will model such a phenomenon 
  with a random variable, which is simply a variable that can take a number of different distinct 
  outcomes depending on the result of an underlying random process.

  For example, say that we have a bag containing 3 red balls and 2 blue balls. If we let X represent the 
  random variable corresponding to the color of a drawn ball, then the probability of each of the two 
  outcomes is given by Pr(X=red)=35 and Pr(X=blue)=25.

  Random variables can be combined to yield new random variables. Returning to the ball example, 
  let Y model the color of a second ball drawn from the bag (without replacing the first ball). 
  The probability of Y being red depends on whether the first ball was red or blue. To represent all 
  outcomes of X and Y, we therefore use a probability tree diagram. This branching diagram represents 
  all possible individual probabilities for X and Y, with outcomes at the endpoints ("leaves") of the 
  tree. The probability of any outcome is given by the product of probabilities along the path from 
  the beginning of the tree; see Figure 2 for an illustrative example.

  An event is simply a collection of outcomes. Because outcomes are distinct, the probability of an 
  event can be written as the sum of the probabilities of its constituent outcomes. For our colored 
  ball example, let A be the event "Y is blue." Pr(A) is equal to the sum of the probabilities of two 
  different outcomes: Pr(X=blue and Y=blue)+Pr(X=red and Y=blue), or 3/10+1/10=2/5 (see Figure 2 above).

  Given: Three positive integers k, m, and n, representing a population containing k+m+n organisms: 
    k individuals are homozygous dominant for a factor, 
    m are heterozygous, 
    and n are homozygous recessive.

  Return: The probability that two randomly selected mating organisms will produce an individual 
  possessing a dominant allele (and thus displaying the dominant phenotype). Assume that any two 
  organisms can mate.

  Sample Dataset
    2 2 2
  Sample Output
    0.78333

  Hint
    Consider simulating inheritance on a number of small test cases in order to check your solution.
  """

  (Maple: 
  > identify(    0.78333);
  47/60

  This model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

% import v3_utils.
% import util.
% import cp.


main => go.

% TODO: Make it a formula based on K,M,and N!
go ?=>
  K = 28, % individuals are homozygous dominant for a factor
  M = 18, % are heterozygous
  N = 20, % are homozygous recessive
  % 28 18 20


  Pop = K + M + N,
  % Simulation
  L = ["AA" : _ in 1..K] ++ ["Aa" : _ in 1..M] ++ ["aa" : _ in 1..N],
  println(l=L),
  Len = L.len,
  C = 0,
  Tot = 0,
  foreach(I in 1..Len, J in 1..Len, I != J)
    Tot := Tot + 4,
    foreach(L1 in L[I], L2 in L[J])
        if L1 == 'A' ; L2 == 'A' then
          C := C + 1
        end
    end   
  end,
  nl,
  println(tot=Tot),
  println(c=C),
  println(pct=(C/Tot)),
  nl,
  
  % We have Pop*Pop - Pop comparisons x 4.
  % Here (3*6 - 6) * 4 = (36 - 6) * 4 = 30 * 4 = 120
  %
  % Each allele will do (Pop - 1) * 4 = (6-1)*4 = 5*4 = 20 comparisons.
  % 
  % of these
  % - AA will dominate always
  % - Aa will dominate half the times
  % - aa will never dominate
  % So: how many combinations of aa are there? And then we take 1-p(aa)
  % For each pair ("pq" and "rs"), we always do 4 checks: pr,ps,qr,qs
  % and if any of these is an A then it dominates
  % This is :
  %   - AA-AA: 4 of 4
  %   - AA-Aa: 4 of 4
  %   - AA-aa: 4 of 4  
  %   - Aa-Aa: 3 of 4
  %   - Aa-aa: 2 of 4  
  %   - aa-aa: 0 of 4
  %

  % Prob
  P = 1.0 - (0.25 * M * (M-1) + M*N + N*(N-1)) / (Pop*(Pop-1)),
  println(p=P),
  nl.
go => true.