gMWT: Generalized Mann-Whitney Type Tests

# Changes:
# 30-06-2013, Added the keepPM option, DF

triple.gmw <- function(X,g,goi,type,nper,alternative,mc,PARAMETERS,output,alg, keepPM,order){

 res <- list()
 pm <- NULL
 if(keepPM) pm <- matrix(NA,ncol=ncol(X),nrow=nper)
 diffTests <- getComb(goi,"triple",order=order)

 METHOD <- c("********* Triple based Test *********")
 DNAME <- PARAMETERS[[1]]
 TEST  <- PARAMETERS[[2]]
 TYPE  <- PARAMETERS[[3]]
 ALTERNATIVE <- PARAMETERS[[4]]
 STATISTIC   <- PARAMETERS[[5]]
 PVAL        <- PARAMETERS[[6]]

 dimX      <- PARAMETERS[[7]]
 XisVector <- PARAMETERS[[8]]

## Case: X is vector
    if(XisVector){
##---------------------------------------------------------------------------------------------------------------------------------------
       if(alternative=="two.sided"){
	  if(type=="permutation"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: permutation, two sided, X is vector
	    for(testRun in 1:nrow(diffTests))
	    {
	      obsValue1 <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,1]],X[g==diffTests[testRun,2]],X[g==diffTests[testRun,3]]))
	      obsValue2 <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,3]],X[g==diffTests[testRun,2]],X[g==diffTests[testRun,1]]))
	      nullDist1 <- perm.triple(X[g==diffTests[testRun,1]],X[g==diffTests[testRun,2]],X[g==diffTests[testRun,3]],nper,algorithm=alg)
	      PVAL <- min(sum(nullDist1>=obsValue1)/nper,sum(nullDist1>=obsValue2)/nper)
	
	      names(PVAL) <- "p.value"
	      STATISTIC <- max(obsValue1,obsValue2)
	      names(STATISTIC) <- "obs.value"
	      ALTERNATIVE <- "two.sided"
	      resTemp<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
	      class(resTemp)<-"htest"
	      
              res[[testRun]] <- resTemp
	      names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," > 1/6 or P",diffTests[testRun,3],diffTests[testRun,2],diffTests[testRun,1]," > 1/6",sep="")
	    }
	    if(output=="min")
	    {
	      resMin <- matrix(NA,ncol=1,nrow=length(res))
	      colnames(resMin) <- "pValues"
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
		resMin[i,1] <- res[[i]]$p.value
	      }
	      res <- resMin
	    }
	  } else if(type=="asymptotic"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: asymptotic, two sided, X is vector
	    for(testRun in 1:nrow(diffTests)){
  	    nx <- length(X[g==diffTests[testRun,1]])
  	    ny <- length(X[g==diffTests[testRun,2]])
  	    nz <- length(X[g==diffTests[testRun,3]])
  	    obsValue1 <- (as.vector(estPI(X=c(X[g==diffTests[testRun,1]], X[g==diffTests[testRun,2]], X[g==diffTests[testRun,3]]),
  	                                  g=c(rep(1,nx), rep(2,ny), rep(3,nz)),
  	                                  type="triple")$probs)-1/6)/sqrt(PHatVar.asymp(nx,ny,nz))
  	    obsValue2 <- (as.vector(estPI(X=c(X[g==diffTests[testRun,3]], X[g==diffTests[testRun,2]], X[g==diffTests[testRun,1]]),
  	                                  g=c(rep(1,nz), rep(2,ny), rep(3,nx)),
  	                                  type="triple")$probs)-1/6)/sqrt(PHatVar.asymp(nz,ny,nx))

	    p1 <- 1-pnorm(obsValue1)
	    p2 <- 1-pnorm(obsValue2)
	    
	    PVAL <- 2*min(p1,p2)
	    	    
	    names(PVAL) <- "p.value"
	    STATISTIC <- max(obsValue1,obsValue2)
	    names(STATISTIC) <- "obs.value"
	    ALTERNATIVE <- "two.sided"
	    resTemp<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
	    class(resTemp)<-"htest"
	    
	    res[[testRun]] <- resTemp
	    names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," > 1/6 or P",diffTests[testRun,3],diffTests[testRun,2],diffTests[testRun,1]," > 1/6",sep="")
	    } 	   
	    if(output=="min"){
	      resMin <- matrix(NA,ncol=1,nrow=length(res))
	      colnames(resMin) <- "pValues"
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
	        resMin[i,1] <- res[[i]]$p.value
	      }
	      res <- resMin
	    }
	    
	  } else {
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: other options, two sided, X is vector
	    res <- c()
	    stop("We do not have this kind of type for the triple test!")
	  }
##---------------------------------------------------------------------------------------------------------------------------------------
       } else if(alternative=="greater"){
	    if(type=="permutation"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: permutation, greater, X is vector
	    for(testRun in 1:nrow(diffTests))
	    {
	      obsValue <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,1]],X[g==diffTests[testRun,2]],X[g==diffTests[testRun,3]]))
	      nullDist <- perm.triple(X[g==diffTests[testRun,1]],X[g==diffTests[testRun,2]],X[g==diffTests[testRun,3]],nper,algorithm=alg)
	      PVAL <- sum(nullDist>=obsValue)/nper
	
	      names(PVAL) <- "p.value"
	      STATISTIC <- obsValue
	      names(STATISTIC) <- "obs.value"
	      ALTERNATIVE <- "greater"
	      resTemp<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
	      class(resTemp)<-"htest"
	      
              res[[testRun]] <- resTemp
	      names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," > 1/6",sep="")
	    }
	    if(output=="min")
	    {
	      resMin <- matrix(NA,ncol=1,nrow=length(res))
	      colnames(resMin) <- "pValues"
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
		resMin[i,1] <- res[[i]]$p.value
	      }
	      res <- resMin
	    }

	  } else if(type=="asymptotic"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: asymptotic, greater, X is vector
	    for(testRun in 1:nrow(diffTests)){
	      nx <- length(X[g==diffTests[testRun,1]])
	      ny <- length(X[g==diffTests[testRun,2]])
	      nz <- length(X[g==diffTests[testRun,3]])
	      obsValue1 <- (as.vector(estPI(X=c(X[g==diffTests[testRun,1]], X[g==diffTests[testRun,2]], X[g==diffTests[testRun,3]]),
	                                    g=c(rep(1,nx), rep(2,ny), rep(3,nz)),
	                                    type="triple")$probs)-1/6)/sqrt(PHatVar.asymp(nx,ny,nz))
	      
	      p1 <- 1-pnorm(obsValue1)
	      
	      PVAL <- p1
	      
	      names(PVAL) <- "p.value"
	      STATISTIC <- obsValue1
	      names(STATISTIC) <- "obs.value"
	      ALTERNATIVE <- "greater"
	      resTemp<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
	      class(resTemp)<-"htest"
	      
	      res[[testRun]] <- resTemp
	      names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," > 1/6",sep="")
	    } 	   
	    if(output=="min"){
	      resMin <- matrix(NA,ncol=1,nrow=length(res))
	      colnames(resMin) <- "pValues"
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
	        resMin[i,1] <- res[[i]]$p.value
	      }
	      res <- resMin
	    }
          } else {
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: other options, greater, X is vector
	    res <- c()
	    stop("We do not have this kind of type for the triple test!")
	  }
       } else if(alternative=="smaller"){
##---------------------------------------------------------------------------------------------------------------------------------------
	   if(type=="permutation"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: permutation, smaller, X is vector
	    for(testRun in 1:nrow(diffTests))
	    {
	      obsValue <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,1]],X[g==diffTests[testRun,2]],X[g==diffTests[testRun,3]]))
	      nullDist <- perm.triple(X[g==diffTests[testRun,1]],X[g==diffTests[testRun,2]],X[g==diffTests[testRun,3]],nper,algorithm=alg)
	      PVAL <- sum(nullDist<obsValue)/nper
	
	      names(PVAL) <- "p.value"
	      STATISTIC <- obsValue
	      names(STATISTIC) <- "obs.value"
	      ALTERNATIVE <- "smaller"
	      resTemp<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
	      class(resTemp)<-"htest"
	      
              res[[testRun]] <- resTemp
	      names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," < 1/6",sep="")
	    }
	    if(output=="min")
	    {
	      resMin <- matrix(NA,ncol=1,nrow=length(res))
	      colnames(resMin) <- "pValues"
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
		resMin[i,1] <- res[[i]]$p.value
	      }
	      res <- resMin
	    }

	  } else if(type=="asymptotic"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: asymptotic, smaller, X is vector
	    for(testRun in 1:nrow(diffTests)){
	      nx <- length(X[g==diffTests[testRun,1]])
	      ny <- length(X[g==diffTests[testRun,2]])
	      nz <- length(X[g==diffTests[testRun,3]])
	      obsValue1 <- (as.vector(estPI(X=c(X[g==diffTests[testRun,1]], X[g==diffTests[testRun,2]], X[g==diffTests[testRun,3]]),
	                                    g=c(rep(1,nx), rep(2,ny), rep(3,nz)),
	                                    type="triple")$probs)-1/6)/sqrt(PHatVar.asymp(nx,ny,nz))
	      
	      p1 <- pnorm(obsValue1)
	      
	      PVAL <- p1
	      
	      names(PVAL) <- "p.value"
	      STATISTIC <- obsValue1
	      names(STATISTIC) <- "obs.value"
	      ALTERNATIVE <- "greater"
	      resTemp<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
	      class(resTemp)<-"htest"
	      
	      res[[testRun]] <- resTemp
	      names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," < 1/6",sep="")
	    } 	   
	    if(output=="min"){
	      resMin <- matrix(NA,ncol=1,nrow=length(res))
	      colnames(resMin) <- "pValues"
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
	        resMin[i,1] <- res[[i]]$p.value
	      }
	      res <- resMin
	    }
          } else {
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: other options, one sided, X is vector
	    res <- c()
	    stop("We do not have this kind of type for the triple test!")
	  }
       } else {
	    res <- c()
	    stop("There is no other option than small, greater or two-sided...")
       }
## Case: X is a matrix
    } else{
##----------------------------------------------------------------------------------------------------------------------------------------
#Preparational things for the case that X is a matrix
    # First, restrict the cores to maximum of possible tests
    if(mc>detectCores()){
	mc <- detectCores()
	warning("You do not have so many cores on this machine! I automatically reduced it to your maximum number ",mc)
    }
    mc <- min(dimX[2],mc)

    if(alternative=="two.sided"){
       	  if(type=="permutation"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: permutation, two sided, X is matrix
	    innerLoop <- function(i,testRun){
             nullDist1 <- perm.triple(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i],nper,algorithm=alg)
             obsValue1 <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i]))
             obsValue2 <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,3],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,1],i]))
             pValue <- min(sum(nullDist1>=obsValue1)/nper,sum(nullDist1>=obsValue2)/nper)
             return(list(pValue=pValue,obsValue=max(obsValue1,obsValue2)))
            }

	    innerLoopPM <- function(i,testRun){
             nullDist1 <- perm.triple(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i],nper,algorithm=alg)
             obsValue1 <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i]))
             obsValue2 <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,3],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,1],i]))
             pValue <- min(sum(nullDist1>=obsValue1)/nper,sum(nullDist1>=obsValue2)/nper)
             return(list(pValue=pValue,obsValue=max(obsValue1,obsValue2), nullDist=nullDist1))
            }

	    if(keepPM){
	        nullDistRES <- list()
		STATISTIC <- list()
		for(i in 1:nrow(diffTests)){
		  nullDistRES[[i]] <- matrix(0, ncol=dimX[2],nrow=nper)
		  STATISTIC[[i]] <- c(rep(-1,dimX[2]))
		}
	    }	 

	    for(testRun in 1:nrow(diffTests))
	    { 
	      resTemp <- list()

	      if(keepPM==TRUE){
   	        resInner <-  unlist(mclapply(c(1:dimX[2]),innerLoopPM,testRun=testRun,mc.cores=mc))
		#nullDistRES <- matrix(0, ncol=dimX[2],nrow=nper)
              } else {
   	        resInner <- unlist(mclapply(c(1:dimX[2]),innerLoop,testRun,mc.cores=mc))
              }

	      for(i in 1:dimX[2])
	      {
		if(keepPM==TRUE){
                  PVAL <- resInner[nper*(i-1) + 2*(i) - 1]
                  STATISTIC[[testRun]][i] <- resInner[nper*(i-1) + 2*i]
                  nullDistRES[[testRun]][,i] <- resInner[(nper*(i-1) + 2*i + 1):(nper*i + 2*i)]
                } else {
		  PVAL <- resInner[2*i-1]
                  STATISTIC <- resInner[2*i]
		}
		obsValue <- STATISTIC
		names(PVAL) <- "p.value"
		ALTERNATIVE <- "two.sided"
		names(STATISTIC) <- "obs.value"
		resTemp[[i]]<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
		class(resTemp[[i]])<-"htest"	    
	      }
	     res[[testRun]] <- resTemp
	     names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," > 1/6 or P",diffTests[testRun,3],diffTests[testRun,2],diffTests[testRun,1]," > 1/6",sep="")
	    }
	    if(output=="min")
	    {
	      resMin <- matrix(NA,ncol=dimX[2],nrow=length(res))
	      colnames(resMin) <- colnames(X)
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
		for(j in 1:dimX[2])
		{
		  resMin[i,j] <- res[[i]][[j]]$p.value
		}
	      }
	      res <- resMin
	    }
	  } else if(type=="asymptotic"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: asymptotic, two sided, X is matrix
	   innerLoop <- function(i,testRun){
	     
	      nx <- length(X[g==diffTests[testRun,1],i])
	      ny <- length(X[g==diffTests[testRun,2],i])
	      nz <- length(X[g==diffTests[testRun,3],i])
	      
	      obsValue1 <- (as.vector(estPI(X=c(X[g==diffTests[testRun,1],i], X[g==diffTests[testRun,2],i], X[g==diffTests[testRun,3],i]),
	                                    g=c(rep(1,nx), rep(2,ny), rep(3,nz)),
	                                    type="triple")$probs)-1/6)/sqrt(PHatVar.asymp(nx,ny,nz))
	      obsValue2 <- (as.vector(estPI(X=c(X[g==diffTests[testRun,3],i], X[g==diffTests[testRun,2],i], X[g==diffTests[testRun,1],i]),
	                                    g=c(rep(1,nz), rep(2,ny), rep(3,nx)),
	                                    type="triple")$probs)-1/6)/sqrt(PHatVar.asymp(nz,ny,nx))

	      p1 <- 1-pnorm(obsValue1)
	      p2 <- 1-pnorm(obsValue2)
	      
	      pValue <- min(p1,p2)
	      
	      return(list(pValue=pValue,obsValue=max(obsValue1,obsValue2)))
	    }
	    
	    for(testRun in 1:nrow(diffTests))
	    { 
	      resTemp <- list()

	      resInner <- unlist(mclapply(c(1:dimX[2]),innerLoop,testRun,mc.cores=mc))
	      
	      for(i in 1:dimX[2])
	      {
	        PVAL <- resInner[2*i-1]
	        STATISTIC <- resInner[2*i]
	        obsValue <- STATISTIC
	        names(PVAL) <- "p.value"
	        ALTERNATIVE <- "two.sided"
	        names(STATISTIC) <- "obs.value"
	        resTemp[[i]]<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
	        class(resTemp[[i]])<-"htest"	    
	      }
	      res[[testRun]] <- resTemp
	      names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," > 1/6 or P",diffTests[testRun,3],diffTests[testRun,2],diffTests[testRun,1]," > 1/6",sep="")
	    }
	    if(output=="min")
	    {
	      resMin <- matrix(NA,ncol=dimX[2],nrow=length(res))
	      colnames(resMin) <- colnames(X)
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
	        for(j in 1:dimX[2])
	        {
	          resMin[i,j] <- res[[i]][[j]]$p.value
	        }
	      }
	      res <- resMin
	    }
          } else {
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: other options, two sided, X is vector
	    res <- c()
	    stop("We do not have this kind of type for the triple test!,O,T,M")
	  }
    } else if(alternative=="greater"){
	  if(type=="permutation"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: permutation, greater, X is matrix
	      # Define the function, that is performed for column i (important for parallelization)
	   innerLoop <- function(i,testRun){
             nullDist <- perm.triple(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i],nper,algorithm=alg)
             obsValue <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i]))
             pValue <- sum(nullDist>=obsValue)/nper
	     return(list(pValue=pValue,obsValue=obsValue))
            }

	   innerLoopPM <- function(i,testRun){
             nullDist <- perm.triple(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i],nper,algorithm=alg)
             obsValue <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i]))
             pValue <- sum(nullDist>=obsValue)/nper
	     return(list(pValue=pValue,obsValue=obsValue, nullDist=nullDist))
            }

	    if(keepPM){
	        nullDistRES <- list()
		STATISTIC <- list()
		for(i in 1:nrow(diffTests)){
		  nullDistRES[[i]] <- matrix(0, ncol=dimX[2],nrow=nper)
		  STATISTIC[[i]] <- c(rep(-1,dimX[2]))
		}
	    }

	    for(testRun in 1:nrow(diffTests))
	    { 
	      resTemp <- list()

	      if(keepPM==TRUE){
   	        resInner <-  unlist(mclapply(c(1:dimX[2]),innerLoopPM,testRun=testRun,mc.cores=mc))
		#nullDistRES <- matrix(0, ncol=dimX[2],nrow=nper)
              } else {
   	        resInner <-  unlist(mclapply(c(1:dimX[2]),innerLoop,testRun=testRun,mc.cores=mc))
              }

	      for(i in 1:dimX[2])
	      {

		if(keepPM==TRUE){
                  PVAL <- resInner[nper*(i-1) + 2*(i) - 1]
                  STATISTIC[[testRun]][i] <- resInner[nper*(i-1) + 2*i]
                  nullDistRES[[testRun]][,i] <- resInner[(nper*(i-1) + 2*i + 1):(nper*i + 2*i)]
                } else {
		  PVAL <- resInner[2*i-1]
                  STATISTIC <- resInner[2*i]
		}
		obsValue <- STATISTIC
		names(PVAL) <- "p.value"
		ALTERNATIVE <- "greater"
		names(STATISTIC) <- "obs.value"
		resTemp[[i]]<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
		class(resTemp[[i]])<-"htest"	    
	      }
	     res[[testRun]] <- resTemp
	     names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," > 1/6",sep="")
	    }
	    if(output=="min")
	    {
	      resMin <- matrix(NA,ncol=dimX[2],nrow=length(res))
	      colnames(resMin) <- colnames(X)
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
		for(j in 1:dimX[2])
		{
		  resMin[i,j] <- res[[i]][[j]]$p.value
		}
	      }
	      res <- resMin
	    }
	  } else if(type=="asymptotic"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: asymptotic, greater, X is matrix
	    innerLoop <- function(i,testRun){
	      
	      nx <- length(X[g==diffTests[testRun,1],i])
	      ny <- length(X[g==diffTests[testRun,2],i])
	      nz <- length(X[g==diffTests[testRun,3],i])
	      
	      obsValue1 <- (as.vector(estPI(X=c(X[g==diffTests[testRun,1],i], X[g==diffTests[testRun,2],i], X[g==diffTests[testRun,3],i]),
	                                    g=c(rep(1,nx), rep(2,ny), rep(3,nz)),
	                                    type="triple")$probs)-1/6)/sqrt(PHatVar.asymp(nx,ny,nz))
	      
	      pValue <- 1-pnorm(obsValue1)
	     
	      return(list(pValue=pValue,obsValue=obsValue1))
	    }
	    
	    for(testRun in 1:nrow(diffTests))
	    { 
	      resTemp <- list()
	      
	      resInner <- unlist(mclapply(c(1:dimX[2]),innerLoop,testRun,mc.cores=mc))
	      
	      for(i in 1:dimX[2])
	      {
	        PVAL <- resInner[2*i-1]
	        STATISTIC <- resInner[2*i]
	        obsValue <- STATISTIC
	        names(PVAL) <- "p.value"
	        ALTERNATIVE <- "greater"
	        names(STATISTIC) <- "obs.value"
	        resTemp[[i]]<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
	        class(resTemp[[i]])<-"htest"	    
	      }
	      res[[testRun]] <- resTemp
	      names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," > 1/6",sep="")
	    }
	    if(output=="min")
	    {
	      resMin <- matrix(NA,ncol=dimX[2],nrow=length(res))
	      colnames(resMin) <- colnames(X)
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
	        for(j in 1:dimX[2])
	        {
	          resMin[i,j] <- res[[i]][[j]]$p.value
	        }
	      }
	      res <- resMin
	    }
	  } else if(type=="asymptotic2"){
	    #----------------------------------------------------------------------------------------------------------------------------------------
	    # Case: asymptotic, greater, X is matrix
	    innerLoop <- function(i,testRun){
	      
	      
	      obsValue1 <- (as.vector(estPI(X=c(X[g==diffTests[testRun,1],i], X[g==diffTests[testRun,2],i], X[g==diffTests[testRun,3],i]),
	                                    g=c(rep(1,nx), rep(2,ny), rep(3,nz)),
	                                    type="triple")$probs)-1/6)/sqrt(PHatVar.asymp(nx,ny,nz))
	      
	      pValue <- 1-pnorm(obsValue1)
	      
	      return(list(pValue=pValue,obsValue=obsValue1))
	    }
	    
	    for(testRun in 1:nrow(diffTests))
	    { 
	      resTemp <- list()
	      nx <- sum(g==diffTests[testRun,1])
	      ny <- sum(g==diffTests[testRun,2])
	      nz <- sum(g==diffTests[testRun,3])
	      
	      obsValue1 <- (as.vector(estPI(X=X, g=g, type="triple")$probs)-1/6)/sqrt(PHatVar.asymp(nx,ny,nz))
	      
	      pValue <- 1-pnorm(obsValue1)
	      
	      resInner <- c(rbind(pValue, obsValue1))	      

	      for(i in 1:dimX[2])
	      {
	        PVAL <- resInner[2*i-1]
	        STATISTIC <- resInner[2*i]
	        obsValue <- STATISTIC
	        names(PVAL) <- "p.value"
	        ALTERNATIVE <- "greater"
	        names(STATISTIC) <- "obs.value"
	        resTemp[[i]]<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
	        class(resTemp[[i]])<-"htest"	    
	      }
	      res[[testRun]] <- resTemp
	      names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," > 1/6",sep="")
	    }
	    if(output=="min")
	    {
	      resMin <- matrix(NA,ncol=dimX[2],nrow=length(res))
	      colnames(resMin) <- colnames(X)
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
	        for(j in 1:dimX[2])
	        {
	          resMin[i,j] <- res[[i]][[j]]$p.value
	        }
	      }
	      res <- resMin
	    }
	  } else {
	    res <- c()
	    stop("We do not have this kind of type for the UIT!,O,G,M")
	  }
    } else if(alternative=="smaller"){
	  if(type=="permutation"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: permutation, smaller, X is matrix
	    innerLoop <- function(i,testRun){
             nullDist <- perm.triple(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i],nper,algorithm=alg)
             obsValue <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i]))
             pValue <- sum(nullDist<obsValue)/nper
	     return(list(pValue=pValue,obsValue=obsValue))
            }

	    innerLoopPM <- function(i,testRun){
             nullDist <- perm.triple(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i],nper,algorithm=alg)
             obsValue <- as.numeric(getP.Cnaive(X[g==diffTests[testRun,1],i],X[g==diffTests[testRun,2],i],X[g==diffTests[testRun,3],i]))
             pValue <- sum(nullDist<obsValue)/nper
	     return(list(pValue=pValue,obsValue=obsValue, nullDist=nullDist))
            }

	    if(keepPM){
	        nullDistRES <- list()
		STATISTIC <- list()
		for(i in 1:nrow(diffTests)){
		  nullDistRES[[i]] <- matrix(0, ncol=dimX[2],nrow=nper)
		  STATISTIC[[i]] <- c(rep(-1,dimX[2]))
		}
	    }

	    for(testRun in 1:nrow(diffTests))
	    { 
	      resTemp <- list()

	      if(keepPM==TRUE){
   	        resInner <-  unlist(mclapply(c(1:dimX[2]),innerLoopPM,testRun=testRun,mc.cores=mc))
		#nullDistRES <- matrix(0, ncol=dimX[2],nrow=nper)
              } else {
   	        resInner <-  unlist(mclapply(c(1:dimX[2]),innerLoop,testRun=testRun,mc.cores=mc))
              }

	      for(i in 1:dimX[2])
	      {

		if(keepPM==TRUE){
                  PVAL <- resInner[nper*(i-1) + 2*(i) - 1]
                  STATISTIC[[testRun]][i] <- resInner[nper*(i-1) + 2*i]
                  nullDistRES[[testRun]][,i] <- resInner[(nper*(i-1) + 2*i + 1):(nper*i + 2*i)]
                } else {
		  PVAL <- resInner[2*i-1]
                  STATISTIC <- resInner[2*i]
		}
		obsValue <- STATISTIC
		names(PVAL) <- "p.value"
		ALTERNATIVE <- "smaller"
		names(STATISTIC) <- "obs.value"
		resTemp[[i]]<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
		class(resTemp[[i]])<-"htest"	    
	      }
	     res[[testRun]] <- resTemp
	     names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," < 1/6",sep="")
	    }
	    if(output=="min")
	    {
	      resMin <- matrix(NA,ncol=dimX[2],nrow=length(res))
	      colnames(resMin) <- colnames(X)
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
		for(j in 1:dimX[2])
		{
		  resMin[i,j] <- res[[i]][[j]]$p.value
		}
	      }
	      res <- resMin
	    }
	  } else if(type=="asymptotic"){
#----------------------------------------------------------------------------------------------------------------------------------------
# Case: asymptotic, smaller, X is matrix
	    innerLoop <- function(i,testRun){
	      
	      nx <- length(X[g==diffTests[testRun,1],i])
	      ny <- length(X[g==diffTests[testRun,2],i])
	      nz <- length(X[g==diffTests[testRun,3],i])
	      
	      obsValue1 <- (as.vector(estPI(X=c(X[g==diffTests[testRun,1],i], X[g==diffTests[testRun,2],i], X[g==diffTests[testRun,3],i]),
	                                    g=c(rep(1,nx), rep(2,ny), rep(3,nz)),
	                                    type="triple")$probs)-1/6)/sqrt(PHatVar.asymp(nx,ny,nz))
	      
	      pValue <- pnorm(obsValue1)
	      
	      return(list(pValue=pValue,obsValue=obsValue1))
	    }
	    
	    for(testRun in 1:nrow(diffTests))
	    { 
	      resTemp <- list()
	      
	      resInner <- unlist(mclapply(c(1:dimX[2]),innerLoop,testRun,mc.cores=mc))
	      
	      for(i in 1:dimX[2])
	      {
	        PVAL <- resInner[2*i-1]
	        STATISTIC <- resInner[2*i]
	        obsValue <- STATISTIC
	        names(PVAL) <- "p.value"
	        ALTERNATIVE <- "smaller"
	        names(STATISTIC) <- "obs.value"
	        resTemp[[i]]<-c(list(method=METHOD,data.name=DNAME,alternative=ALTERNATIVE,statistic=STATISTIC,test=TEST,p.value=PVAL,type=TYPE))
	        class(resTemp[[i]])<-"htest"	    
	      }
	      res[[testRun]] <- resTemp
	      names(res)[testRun] <- paste("H1: P",diffTests[testRun,1],diffTests[testRun,2],diffTests[testRun,3]," < 1/6",sep="")
	    }
	    if(output=="min")
	    {
	      resMin <- matrix(NA,ncol=dimX[2],nrow=length(res))
	      colnames(resMin) <- colnames(X)
	      rownames(resMin) <- names(res)
	      for(i in 1:length(res))
	      {
	        for(j in 1:dimX[2])
	        {
	          resMin[i,j] <- res[[i]][[j]]$p.value
	        }
	      }
	      res <- resMin
	    }
          } else {
	    res <- c()
	    stop("We do not have this kind of type for the triple test!,O,S,M")
	  }
    } else {
	    res <- c()
	    stop("There are no other alternatives possible, sorry! All other....")
     }
  }
  if(type=="permutation"){
    ifelse(keepPM,res <- list(p.values=res, nullDist=nullDistRES, obsValue=obsValue), res <- list(p.values=res))
  } else {
    res <- list(p.values=res)
  }
  res
}