##############################################################
# fbsthw:
#   FBST for Hardy-Weinberg Equilibrium Test 
#
# Description: 
#   Performs the test of Hardy-Weinberg Equilibrium (HWE) in a locus with
#   several alleles A_1, A_2,..., A_k, k>=2. 
#    
# Usage:
#   fbsthw (x, prec=0.005, nmcmin=3000*k, nmcmax=30000*k)  
#  
# Arguments:
#   x: vector of genotype counts in a lower triangular form:
#          x = c(x_11,
#                x_21, x_22,
#                x_31, x_32, x_33
#                ...
#                x_k1, x_k2, .... x_kk)
#   prec: parameter for 99% confidence interval for the evidence measure:
#                       CI99% = evid+-prec                 (default=5e-3)
#   nmcmin: minimum number of iterations in MC integration (default=3000*k)
#   nmcmax: maximum number of iterations in MC integration (default=30000*k)
#   
# Details: 
#   Given k alleles A_1, A_2, ... A_k in a locus, and their genotype counts 
#   in lower triangular form x = (x_{ij}), 1<=j<=i<=k, fbstwh computes the
#   evidence for the Hardy-Weinberg Equilibrium) HWE. 
# 
#   The parameter of interest, denoted by theta_{ij}, is the population relative 
#   frequencies of genotypes A_{i}A_{j} (1<=j<=i<=k). Assuming a multinomial
#   model for x and a Dirichlet(1,1,...,1) priori for theta, the posterior p.d.f
#   is f(theta | x) = Dirichlet( x_{11}+1, x_{12}+1, ..., x_{1k}+1 )
#
#   Parameter space:
#   THETA = {theta = (theta_{ij}), 1<=j<=i<=k | 0<=theta_{ij}<=1, sum(theta)=1} 
#   
#   Parameter space under HWE: 
#   H = {theta in THETA | there exists a vector p=(p_1, p_2, ..., p_k), 
#            0<=p<=1, sum(p)=1, such that:  
#            theta_{ij} = (p_i)^2 if i=j; theta_{ij} = 2*p_i*p_j if i<>j}
#
#   The evidence measure against the hypothesis is computed in two steps:
#      1) compute f_0 = max_H f(theta|x)
#      2) compute the integral   
#            evb = int_T f(theta|x) dtheta / int_THETA f(theta|x) dtheta  
#         where T = {theta in THETA | f(theta | x) > f_0}
#   The evidence measure suporting HWE is:
#                   evid = 1 - evb
# 
#   Further details are found in Lauretto et al (2009) (Reference 1). 
#
# Value:
#   A list containing the following components:
#     evid: the evidence measure supporting HWE hypothesis (see Reference 1)
#     p.value: the significance (p-value) for the test. 
#
# Dependencies: 
#   The program fbsthw.exe must be copied in the same directory as this module
#   fbsthw.r.
#
# References:
#    1. M.S.Lauretto, F.Nakano, S.R.Faria Jr, C.A.B.Pereira, J.M.Stern. 
#       A straightforward multiallelic significance test for the Hardy-Weinberg
#       equilibrium law. Genetics and Molecular Biology 32(3): 619-625, 2009.
#
#    2. C.A.B.Pereira, J.M.Stern. Evidence and credibility: full Bayesian significance 
#       test for precise hypotheses. Entropy Journal 1: 69-80, 1999.
#
# Examples:
#
#   # 1) Louis and Dempster (1987)
#   counts = c(
#     0,			
#     3,  1,		
#     5,  18, 1,	
#     3,  7,  5,  2)
#
#   EvFBST = fbsthw(counts)
#   print(EvFBST$p.value)
#
#   The result should be approx. 0.0034
#
#
#   # 2) Guo and Thompson (1992)
#   counts = c(
#      3,
#      4, 2,
#      2, 2, 2,
#      3, 3, 2, 1,
#      0, 1, 0, 0, 0,
#      0, 0, 0, 0, 0, 1,
#      0, 0, 1, 0, 0, 0, 0,
#      0, 0, 0, 2, 1, 0, 0, 0)
#
#   EvFBST = fbsthw(counts)
#   print(EvFBST$p.value)
#
#   The result should be approx. 0.5570
#
#
#   # 3) Example 3 Wordsworth et al (1992) 
#   counts = c(
#      5,			
#      40, 12,		
#      6,  32,  2,	
#      30, 55, 15, 33)
#
#   EvFBST = fbsthw(counts)
#   print(EvFBST$p.value)
#
#   The result will be approx. 5.483e-05 
#
##############################################################

cumEvidAg = function(evidag, dimens, dimensh) {
   p = pchisq(qchisq(evidag, dimens),dimens-dimensh)
   return (p)
}


fbsthw <- function(x, prec=5e-3, nmcmin=0, nmcmax=0){
   #  Funcao para chamada do FBST para teste de equilibrio de hardy-weinberg 
   #  x eh um vetor representando as frequencias genotipicas (na forma de 
   #  uma matriz triangular inferior ou seja, x=c(n_11,n_21,n_31,...,n_k1,n_22,n_32,...)

   dimens = length(x)
 
   if (dimens>3) {
     qtalleles = suppressWarnings(as.integer(polyroot(c(-2*dimens, 1, 1))[1]))
   } else { qtalleles = 2 }
   
   if (nmcmin==0) nmcmin <- 5000*dimens
   if (nmcmax==0) nmcmax <- 30000*dimens;  
 
 
#  %%%%%%%%%%%%%%%%%%%%%
#  %
#  %	Parametros para as rotinas C
#  %   qtalelos
#  %   n11 n21 n31 ... nk1 n22 n32 ... nkk
#  %   prec nmcmin nmcmax
#  %
#  %  Saida: Evidencia a favor da hipotese 
#  %

   cat(qtalleles, "\n", file="drvinp.txt")
   cat(x, "\n", file="drvinp.txt", append=TRUE)
   cat(prec, format(c(nmcmin, nmcmax), scientific=FALSE), "\n", file="drvinp.txt", append=TRUE)

   cat(-1, "\n", file="drvout.txt")

   #shell("fbsthw drvinp.txt drvout.txt")       
   try(system("./fbsthw drvinp.txt drvout.txt"))
      
   evid = scan("drvout.txt",what=double(0), nmax=1) 
   p.value = 1 - cumEvidAg(1-evid, dimens, qtalleles)
   
   return(list(evid=evid, p.value=p.value))

}