## genetic model tests
SNPs.test <- function(x.freq)
{
    fCase  <- x.freq[1:3]
    fTotal <- x.freq[7:9]

    res.rec <- prop.trend.test(fCase, fTotal, score=c(0,0,1)) # recessive model
    res.add <- prop.trend.test(fCase, fTotal, score=c(0,1,2)) # additive  model
    res.dom <- prop.trend.test(fCase, fTotal, score=c(0,1,1)) # dominant  model

    T.rec <- res.rec$statistic
    T.add <- res.add$statistic
    T.dom <- res.dom$statistic

    max(c(T.rec, T.add, T.dom))
}


## Calculate statistics
cal.stats <- function(x.freq)
{
    res.temp <- SNPs.test(x.freq)
    norm.obs <- sqrt(res.temp)

    p  <- x.freq[7:9]/sum(x.freq[7:9])
    p0 <- ifelse(p[1]<1e-10, 1e-10, p[1])
    p0 <- ifelse(p[1]>1-1e-10, 1-1e-10, p[1])
    p1 <- ifelse(p[2]<1e-10, 1e-10, p[2])
    p1 <- ifelse(p[2]>1-1e-10, 1-1e-10, p[2])
    p2 <- ifelse(p[3]<1e-10, 1e-10, p[3])
    p2 <- ifelse(p[3]>1-1e-10, 1-1e-10, p[3]) 
    cor.rec.add <- p0*(p1+2*p2)/sqrt(p0*(1-p0))/sqrt(p0*(p1+2*p2)+p2*(p1+2*p0))
    cor.rec.dom <- p0*p2/sqrt(p0*(1-p0))/sqrt(p2*(1-p2))    
    cor.add.dom <- p2*(2*p0+p1)/sqrt(p2*(1-p2))/sqrt(p0*(p1+2*p2)+p2*(p1+2*p0))    
    cor.vector  <- c(1, cor.rec.add, cor.rec.dom, cor.rec.add, 1, cor.add.dom, cor.rec.dom, cor.add.dom, 1)
    cor.matrix  <- matrix(data=cor.vector, nrow=3, ncol=3, byrow=TRUE)    
    
    L <- rep(norm.obs,3)

    p.exact <- 1 - pmvnorm(lower=-L, upper=L, mean=c(0,0,0), corr=cor.matrix, abseps = 1e-8)

    c(norm.obs, p.exact)
}


## generate n*6 table
generate.table <- function(n.case, n.ctrl, maf, lambda)
{
    loc.f <- function(x)
    {
        x/sum(x)
    }

    loc.g <- function(y, m)
    {
        as.vector( rmultinom(1, size=m, prob=y) )
    }

    p <- maf
    n <- length(p)

    p.ctrl <- cbind((1-p)^2, 2*p*(1-p), p^2)
    temp   <- p.ctrl*lambda
    p.case <- t( apply(temp, 1, loc.f) )
    
    f.case <- t( apply(p.case, 1, loc.g, n.case) )
    f.ctrl <- t( apply(p.ctrl, 1, loc.g, n.ctrl) )

    cbind(f.case, f.ctrl, f.case+f.ctrl+1e-10)
}


## 1 REC, 1 ADD and 1 DOM
MAX.PV.3 <- function(n.snp=1e5, n.case=1000, n.ctrl=1000, RR=1.2, outFile="output_1.txt", seed=1)
{
    CAD_CTRL <- read.table(file="//data//liq//MAX_Rank//prog//CAD_CTRL_MAF.txt")
    CAD_CASE <- read.table(file="//data//liq//MAX_Rank//prog//CAD_CASE_MAF.txt")

    m <- 3
    n <- n.snp-m
    maf <- c(CAD_CASE[1:m,7], CAD_CTRL[1:n,7])

    # randomly generate table
    set.seed(seed)
    lambda <- matrix(c(1,1,RR,1,RR,2*RR-1,1,RR,RR,rep(1,n*3)), nrow=n.snp, byrow=TRUE)
    F <- generate.table(n.case, n.ctrl, maf, lambda)

    qet  <- t(apply(F, 1, cal.stats))
    temp <- cbind(sort(qet[,1],decreasing=TRUE), sort(qet[,2]))
    res  <- matrix(data=NA, nrow=m, ncol=2)
    for (i in 1:m)
    {
	for (j in 1:2) res[i,j] <- which(temp[,j]==qet[i,j])[1]
    }
    
    write(t(res), file=outFile, ncolumns=ncol(res), sep="\t")
}


## 1 REC, 2 ADD, 2 MUL and 1 DOM
MAX.PV.6 <- function(n.snp=1e5, n.case=1000, n.ctrl=1000, RR=1.2, outFile="output_1.txt", seed=1)
{
    CAD_CTRL <- read.table(file="//data//liq//MAX_Rank//prog//CAD_CTRL_MAF.txt")
    CAD_CASE <- read.table(file="//data//liq//MAX_Rank//prog//CAD_CASE_MAF.txt")

    m <- 6
    n <- n.snp-m
    maf <- c(CAD_CASE[1:m,7], CAD_CTRL[1:n,7])

    # randomly generate table
    set.seed(seed)
    lambda <- matrix(c(1,1,RR,1,RR,2*RR-1,1,RR,2*RR-1,1,RR,RR^2,1,RR,RR^2,1,RR,RR,rep(1,n*3)), nrow=n.snp, byrow=TRUE)
    F <- generate.table(n.case, n.ctrl, maf, lambda)

    qet  <- t(apply(F, 1, cal.stats))
    temp <- cbind(sort(qet[,1],decreasing=TRUE), sort(qet[,2]))
    res  <- matrix(data=NA, nrow=m, ncol=2)
    for (i in 1:m)
    {
	for (j in 1:2) res[i,j] <- which(temp[,j]==qet[i,j])[1]
    }
    
    write(t(res), file=outFile, ncolumns=ncol(res), sep="\t")
}


## 2 REC, 2 ADD, 2 MUL and 2 DOM
MAX.PV.8 <- function(n.snp=1e5, n.case=1000, n.ctrl=1000, RR=1.2, outFile="output_1.txt", seed=1)
{
    CAD_CTRL <- read.table(file="//data//liq//MAX_Rank//prog//CAD_CTRL_MAF.txt")
    CAD_CASE <- read.table(file="//data//liq//MAX_Rank//prog//CAD_CASE_MAF.txt")

    m <- 8
    n <- n.snp-m
    maf <- c(CAD_CASE[1:m,7], CAD_CTRL[1:n,7])

    # randomly generate table
    set.seed(seed)
    lambda <- matrix(c(1,1,RR,1,1,RR,1,RR,2*RR-1,1,RR,2*RR-1,1,RR,RR^2,1,RR,RR^2,1,RR,RR,1,RR,RR,rep(1,n*3)), nrow=n.snp, byrow=TRUE)
    F <- generate.table(n.case, n.ctrl, maf, lambda)

    qet  <- t(apply(F, 1, cal.stats))
    temp <- cbind(sort(qet[,1],decreasing=TRUE), sort(qet[,2]))
    res  <- matrix(data=NA, nrow=m, ncol=2)
    for (i in 1:m)
    {
	for (j in 1:2) res[i,j] <- which(temp[,j]==qet[i,j])[1]
    }
    
    write(t(res), file=outFile, ncolumns=ncol(res), sep="\t")
}