############################FIRST PART################################
###################Function for simulation############################
######################################################################
SIMULATION <- function(numc, numobs,pdiff,lfc,pzerodiff,n.sim,
                       mu.prob=NULL, pzero.prob=NULL){
  for (k in 1:n.sim){
    
    #############Get simulation data################
    print(k)  
    set.seed(SEEDS[k])
    simdata<- Getsample(numc=numc,numobs=numobs,pdiff=pdiff,lfc=lfc,pzerodiff=pzerodiff,simpool)  
    
    ####################Get simulation parameters###############
    paradata<- simdata$Parameters
    geneid <- c(rep(1,(numc*pdiff)),rep(0,round(numc*(1-pdiff))))
    geneid1 <- c(rep(1,round(numc*pdiff/2)),rep(-1,round(numc*pdiff/2)),rep(0,round(numc*(1-pdiff))))
    
    ##########Data part###############
    ldata<- simdata$simdata
    idint<- rowSums(ldata> 0) > 2
    ldata <- ldata[idint,] 
    paradata <- paradata[idint,] 
    geneid<- geneid[idint]
    geneid1<-geneid1[idint]
    numint <- c(1:length(geneid))
    indall <- c()
    for(i in 1:nrow(ldata)){
      ind <- (!all(ldata[i,][1:numobs]!=0)         #exclude all zero or all non zero data in each group to adapt the other method. 
              &!all(ldata[i,][1:numobs]==0)
              &!all(ldata[i,][(1+numobs):(2*numobs)]!=0)
              &!all(ldata[i,][(1+numobs):(2*numobs)]==0))
      indall<- c(indall, ind)
    }
    paradata <- cbind(paradata,geneid,geneid1,numint)
    paradata <- paradata[indall,]
    data2<-ldata[indall,]
    datafile <- paste("simdata_N_",numobs,"_",k,".Rdata",sep="")
    save(data2,file=datafile)
    fnamep <- paste("paradata_N_",numobs,"_",k,".Rdata",sep="")
    save(paradata,file=fnamep)
    
    #####################Our method####################
    #sim.result <- MF(indata=data2,covar=c(rep(1,numobs),rep(0,numobs)), nc=numobs, nd=numobs,min.non0n=5)
    sim.result <- DASEV.DA(indata=data2, cov.matrix=cbind(rep(1, 2*numobs),c(rep(0,numobs),rep(1,numobs))), 
                           test_cov=2, min.non0n=3, requiredn=10, requiredn2=30, 
                           DL_method= "Fixed Difference", DL_value=0.1, maxit_MLE=100, 
                           maxit=10000, test_model=c("Both", "Mean", "Pzero"), massat=-Inf)
    fname <- paste("simresult_N_",numobs,"_",k,".Rdata",sep="")
    save(sim.result,file=fname)
    
    #####################Other methods####################
    t.pvalue <- c()
    t.var <- c()
    t.mucontrol <- c()
    t.mucase <- c()
    t.pcontrol <- c()
    t.pcase <- c()
    
    t.pvalue.m <- c()
    t.var.m <- c()
    t.mucontrol.m <- c()
    t.mucase.m <- c()
    t.pcontrol.m <- c()
    t.pcase.m <- c()
    
    t.pvalue.p <- c()
    t.var.p <- c()
    t.mucontrol.p <- c()
    t.mucase.p <- c()
    t.pcontrol.p <- c()
    t.pcase.p <- c()
    
    t2.pvalue <- c()
    t2.mucontrol <- c()
    t2.mucase <- c()
    
    twoT.pvalueX1binary <- c()
    twoT.pvalueX1nonpmv <- c()
    twoT.pvalueX2 <- c()
    
    for (i in 1:nrow(data2)){
      data <- as.matrix(data2[i,])
      covar <- as.matrix(c(rep(0,numobs), rep(1,numobs)))
      try <- Mixture(data, covar)
      try2 <- AFT(data, covar)
      trytwoW <- TwoPart(data, covar, test="wilcoxon", point.mass=0)
      trytwoT <- TwoPart(data, covar, test="t.test", point.mass=0)
      
      tpvalue <- try$P.value
      tvar <- try$Estimates[5]^2
      tmucontrol <- try$Estimates[1]
      tmucase <- try$Estimates[3]
      tpcontrol <- 1-try$Estimates[2]
      tpcase <- 1-try$Estimates[4]
      
      tpvalue.m <-  try$P.value.m
      tvar.m <-  try$Estimates.m[5]^2
      tmucontrol.m <-  try$Estimates.m[1]
      tmucase.m <-  try$Estimates.m[3]
      tpcontrol.m <-  1-try$Estimates.m[2]
      tpcase.m <-  1-try$Estimates.m[4]
      
      tpvalue.p <-  try$P.value.p
      tvar.p <-  try$Estimates.p[5]^2
      tmucontrol.p <-  try$Estimates.p[1]
      tmucase.p <-  try$Estimates.p[3]
      tpcontrol.p <-  1-try$Estimates.p[2]
      tpcase.p <-  1-try$Estimates.p[4]
      
      t2pvalue <- pchisq(try2$chi,1,lower=F)
      t2mucontrol <- try2$coefficients[1]
      t2mucase <- try2$coefficients[1]+try2$coefficients[2]
      
      t.pvalue <- c(t.pvalue,tpvalue)
      t.var <- c(t.var,tvar)
      t.mucontrol <- c(t.mucontrol,tmucontrol)
      t.mucase <- c(t.mucase,tmucase)
      t.pcontrol <- c(t.pcontrol,tpcontrol)
      t.pcase <- c(t.pcase,tpcase)
      
      t.pvalue.m <- c(t.pvalue.m,tpvalue.m)
      t.var.m <- c(t.var.m,tvar.m)
      t.mucontrol.m <- c(t.mucontrol.m,tmucontrol.m)
      t.mucase.m <- c(t.mucase.m,tmucase.m)
      t.pcontrol.m <- c(t.pcontrol.m,tpcontrol.m)
      t.pcase.m <- c(t.pcase.m,tpcase.m)
      
      t.pvalue.p <- c(t.pvalue.p,tpvalue.p)
      t.var.p <- c(t.var.p,tvar.p)
      t.mucontrol.p <- c(t.mucontrol.p,tmucontrol.p)
      t.mucase.p <- c(t.mucase.p,tmucase.p)
      t.pcontrol.p <- c(t.pcontrol.p,tpcontrol.p)
      t.pcase.p <- c(t.pcase.p,tpcase.p)
      
      t2.pvalue <- c(t2.pvalue,t2pvalue)
      t2.mucontrol <- c(t2.mucontrol,t2mucontrol)
      t2.mucase <- c(t2.mucase,t2mucase)
      
      twoT.pvalueX1binary <- c(twoT.pvalueX1binary,trytwoT$X1pvalue[1])
      twoT.pvalueX1nonpmv <- c(twoT.pvalueX1nonpmv,trytwoT$X1pvalue[2])
      twoT.pvalueX2 <- c(twoT.pvalueX2,trytwoT$X2pvalue)
    }
    
    other.result <- cbind(t.pvalue,t.var,t.mucontrol,t.mucase,t.pcontrol,t.pcase,
                          t.pvalue.m,t.var.m,t.mucontrol.m,t.mucase.m,t.pcontrol.m,t.pcase.m,
                          t.pvalue.p,t.var.p,t.mucontrol.p,t.mucase.p,t.pcontrol.p,t.pcase.p,
                          t2.pvalue,t2.mucontrol,t2.mucase,
                          twoT.pvalueX1binary,twoT.pvalueX1nonpmv,twoT.pvalueX2)
    rownames(other.result) <- c(1:nrow(data2))
    
    fname1 <- paste("otherresult_N_",numobs,"_",k,".Rdata",sep="")
    save(other.result,file=fname1)
  }
  
}

######################################################################





###########################SECOND PART################################
#######################Preform Simualtions############################
######################################################################
load("simpool.Rdata")
RNGversion("3.5.3") ######Set R version to replicate DASEV results.
set.seed(689515614) 
SEEDS <- sample(11111111:99999999, 100, replace=FALSE)

###################Simulation scenario one############################
setwd("C:/DASEV/S1")
SIMULATION(numc=5000,numobs=10,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=NULL,n.sim = 100,mu.prob=NULL, pzero.prob=NULL)
SIMULATION(numc=5000,numobs=20,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=NULL,n.sim = 100,mu.prob=NULL, pzero.prob=NULL)
SIMULATION(numc=5000,numobs=100,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=NULL,n.sim = 100,mu.prob=NULL, pzero.prob=NULL)
SIMULATION(numc=5000,numobs=200,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=NULL,n.sim = 100,mu.prob=NULL, pzero.prob=NULL)
######################################################################

###################Simulation scenario two############################
setwd("C:/DASEV/S2")
SIMULATION(numc=5000,numobs=10,pdiff=0.1,lfc = NULL,
           pzerodiff=c(2,0.5),n.sim = 100,mu.prob=NULL,pzero.prob=NULL)
SIMULATION(numc=5000,numobs=20,pdiff=0.1,lfc = NULL,
           pzerodiff=c(2,0.5),n.sim = 100,mu.prob=NULL,pzero.prob=NULL)
SIMULATION(numc=5000,numobs=100,pdiff=0.1,lfc = NULL,
           pzerodiff=c(2,0.5),n.sim = 100,mu.prob=NULL,pzero.prob=NULL)
SIMULATION(numc=5000,numobs=200,pdiff=0.1,lfc = NULL,
           pzerodiff=c(2,0.5),n.sim = 100,mu.prob=NULL,pzero.prob=NULL)
######################################################################

###################Simulation scenario three##########################
setwd("C:/DASEV/S3")
SIMULATION(numc=5000,numobs=10,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=NULL, pzero.prob=NULL)
SIMULATION(numc=5000,numobs=20,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=NULL, pzero.prob=NULL)
SIMULATION(numc=5000,numobs=100,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=NULL, pzero.prob=NULL)
SIMULATION(numc=5000,numobs=200,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=NULL, pzero.prob=NULL)
######################################################################

###################Simulation scenario four##########################
setwd("C:/DASEV/S4")
SIMULATION(numc=5000,numobs=10,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=c(0.9,0.1), 
           pzero.prob=c(0.9,0.1))
SIMULATION(numc=5000,numobs=20,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=c(0.9,0.1), 
           pzero.prob=c(0.9,0.1))
SIMULATION(numc=5000,numobs=100,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=c(0.9,0.1), 
           pzero.prob=c(0.9,0.1))
SIMULATION(numc=5000,numobs=200,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=c(0.9,0.1), 
           pzero.prob=c(0.9,0.1))
######################################################################

###################Simulation scenario five##########################
setwd("C:/DASEV/S5")
SIMULATION(numc=5000,numobs=10,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=c(0.9,0.1), 
           pzero.prob=c(0.1,0.9))
SIMULATION(numc=5000,numobs=20,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=c(0.9,0.1), 
           pzero.prob=c(0.1,0.9))
SIMULATION(numc=5000,numobs=100,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=c(0.9,0.1), 
           pzero.prob=c(0.1,0.9))
SIMULATION(numc=5000,numobs=200,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=c(2,0.5),n.sim=100,mu.prob=c(0.9,0.1), 
           pzero.prob=c(0.1,0.9))
######################################################################

###################Simulation scenario six##########################
setwd("C:/DASEV/S6")
simpool[,1] <- -Inf ###Set DL to be -Inf for no TPMVs
SIMULATION(numc=5000,numobs=10,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=NULL,n.sim=100,mu.prob=NULL,pzero.prob=NULL)
SIMULATION(numc=5000,numobs=20,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=NULL,n.sim=100,mu.prob=NULL,pzero.prob=NULL)
SIMULATION(numc=5000,numobs=100,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=NULL,n.sim=100,mu.prob=NULL,pzero.prob=NULL)
SIMULATION(numc=5000,numobs=200,pdiff=0.1,lfc = c(log(2),-log(2)),
           pzerodiff=NULL,n.sim=100,mu.prob=NULL,pzero.prob=NULL)
######################################################################