# PowerAnalysis_script.R

# Created by Sajal Kumar
# Copyright (c) NMSU Song lab

# A script to compute the statistical power of sharma-song test in relation to sample size, given
# table size, effect size and false positive rate.

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
require(DiffXCoExpNet)

nthreads = 10

PowerAnalysis_script = function(){
  
  # 3 tables sizes
  r = c(3:5)
  # 100K simulations
  num = 100000
  
  # find common threshold for the three table sizes at a fixed sample size
  thres = rep(0, length(r))
  tpr_collect = vector("list", length = length(r))
  
  for(j in 1:length(r)){
    
    shso_dist =  SharmaSongEsDist(N = c(50,50), n_conditions = 2, nthreads = 10, 
                                  r = r[j], s = r[j], d = NULL, num = num)
    thres[j] = sort(shso_dist[,4])[round(0.6*length(shso_dist[,4]))]
  }
  
  thres = median(thres)
  
  # compute statistical power
  for(j in 1:length(r)){
    
    cat("Row :",r[j],"\n")
    
    # get valid samples
    n = c(seq(r[j], 19, 1), seq(20, 100, 5))
    
    tpr = rep(0, length(n))
    for(i in 1:length(n)){
      
      shso_dist =  SharmaSongEsDist(N = c(n[i], n[i]), n_conditions = 2, nthreads = 10, 
                                    r = r[j], s = r[j], d = NULL, num = num)
      tpr[i] = sum(shso_dist[,3] < 0.05 & shso_dist[,4] > thres)/sum(shso_dist[,4] > thres)
      
    }
    
    tpr_collect[[j]] = tpr
  }
  
  # plot results
  pdf("../Results/SharmaSongPowerAnalysis.pdf", width=10)
  
  col = c("red", "blue", "green")
  n = c(seq(r[1], 19, 1), seq(20, 100, 5))
  
  par(mar=c(5,5,3,3))
  plot(0, xlim = c(0, length(n)), ylim = c(0, 1), ylab = "Power", xlab = "Sample size", type = "n",
       cex.lab = 2, xaxt = "n", main = "", cex.axis=1.5)
  axis(side = 1, at = c(1:length(n)), labels = n, las=2, cex.axis=1.5)
  grid()
  ts = rep("", length(tpr_collect))
  for(j in 1:length(tpr_collect)){
    lines(tpr_collect[[j]]~c(j:length(n)), col=col[j], lwd = 3)
    ts[j] = paste0(r[j],"x",r[j])
  }
  legend("bottomright", legend = ts, fill = col, bty = "n",
         title = "Table size", cex=2)
  
  dev.off()
  
}