#!/usr/local/bin/Rscript
args <- commandArgs(TRUE)
if(length(args) < 1){
    stop("Please please pass in the name of the dataset(zhang, moran, or scherzer) as an argument")
}
if(args[1] != "zhang" & args[1] != "moran" & args[1] != "scherzer"){
    stop("This script only take argument zhang, moran, or scherzer")
}
dataname <- args[1]
cached = FALSE
## check if cached data exists and has the right number of gene sets
if(file.exists(paste0("./pval_comparison_", dataname, ".rda"))){
    A <- readRDS(paste0("./pval_comparison_", dataname, ".rda"))
    if(dim(A)[1] != 6180){
        warning(paste("The cached file",paste0("./pval_comparison_", dataname, ".rda") ,"doesn't contain all gene sets, recomputing."))
    }else{
        print(paste("loading from", paste0("./pval_comparison_", dataname, ".rda")))
        cached = TRUE
    }
}

if(!cached){
    library(pipeGS)
    load(paste0("./", dataname, "_all.RData"))
    pval.linear <- rep(NA, length(gs))
    #'@param X binary indicator
    #'@param Y the standardized genome measurement matrix
    #'@param eset1 The original geneset measurement data, with gene names extractable by calling featureNames(eset1)
    #'@param gs the gene set definitions with entrezids
    #'@param i The index of the gene set, can range from 1 to length(gs)
    #'@param entrezids.map defines the gene name and entrizid map.
    get.pval.linear.geneset.phat13 <- function(X, Y, eset1, gs, i, entrezids.map){
        Y.geneset <- Y[, which(colnames(Y) %in% entrezids.map[which(entrezids.map[,2] %in% gs[[i]]),1])]
        time.saddle1 <- system.time(phat.saddle1 <- get.phat.linear(X, Y.geneset, method = "saddlepoint")$saddle1)[[3]]
        time.phat1 <- system.time(phat1 <- get.phat.linear(X, Y.geneset, method = "phat1")$phat1)[[3]]
        time.phat2 <- system.time(phat2 <- get.phat.linear(X, Y.geneset, method = "phat2")$phat2)[[3]]
        time.phat3 <- system.time(phat3 <- get.phat.linear(X, Y.geneset, method = "phat3")$phat3)[[3]]
        return(c(phat.saddle1 = phat.saddle1, phat1 = phat1, phat2 = phat2, phat3 = phat3, time.saddle1 = time.saddle1, time.phat1 = time.phat1, time.phat2 = time.phat2, time.phat3 = time.phat3))
    }
    A <- c()
    t0 <- proc.time()[[3]]
    for(i in 1:length(gs)){
        print(paste(dataname, "geneset", i))
        A <- rbind(A, get.pval.linear.geneset.phat13(X, Y, eset1, gs, i = i, entrezids.map))
        if(i%%10 == 1){
            print(paste0("Done ",floor(i/length(gs)*10000)/100, "%. ", "Estimated remaining time ", ceiling((proc.time()[[3]]-t0)/i*(length(gs) - i)/60), " min"))
        }
    }
    A <- data.frame(A)
    A["phat.mc"] <- phat.mc
    saveRDS(data.frame(A), paste0("./pval_comparison_", dataname, ".rda"))
}

## filter phat.mc to be larger than 0 and less than 0.05
print(sort(A$phat.mc)[1:10])
A.small <- A[A$phat.mc < 0.05,]
if(dataname == "moran"){
    A.tiny <- A[A$phat.mc < 10^-4,]
}else if(dataname == "scherzer"){
    A.tiny <- A[A$phat.mc < 10^-3,]
}else if(dataname == "zhang"){
    A.tiny <- A.small
    print("zhang dataset pvalue matrix")
    print(A.small[,c("phat.saddle1", "phat1", "phat2", "phat3", "phat.mc")])
}
corr.mat <- data.frame(matrix(0,4,4))
estimates <- c("phat1", "phat2", "phat3", "phat.saddle1")
for(i in 1:4){
    corr.mat[1,i] <- cor(A.small[estimates[i]], A.small["phat.mc"], use = "complete", method = "pearson")
    corr.mat[2,i] <- cor(A.small[estimates[i]], A.small["phat.mc"], use = "complete", method = "kendall")
    corr.mat[3,i] <- cor(A.tiny[estimates[i]], A.tiny["phat.mc"], use = "complete", method = "pearson")
    corr.mat[4,i] <- cor(A.tiny[estimates[i]], A.tiny["phat.mc"], use = "complete", method = "kendall")
}
colnames(corr.mat) <- c("phat1", "phat2", "phat3", "phat.saddle")
corr.mat["dataset"] <- c(paste(dataname, "small"), paste(dataname, "small"), paste(dataname, "tiny"), paste(dataname, "tiny"))
corr.mat["correlation"] <- c("Pearson", "Kendall", "Pearson", "Kendall")
corr.mat["# genesets"] <- c(dim(A.small)[1], dim(A.small)[1], dim(A.tiny)[1], dim(A.tiny)[1])
corr.mat <- corr.mat[,c(5,6,7,1,2,3,4)]
options(digits = 4)
print("Correlation with phat.mc")
print(corr.mat)
print("average running time")
sapply(A, mean, 1)[5:8]