# TODO: Add comment
# 
# Author: E.Korsching 18.03.2012
###############################################################################


###### start sampling process, bootstrap or permutation test ######
## e.g. 10000 permutations
## to avoid a loop along all permutations, compute first a matrix with 10000 permutated rows like 
##					3	7	2	6	1	...
##					1	8	3	10	2	...
## compute a second matrix , repeat one gene row 10000 times and apply the permutation matrix to it
## use contrast matrices to compute all possible comparisons 
##	(6 groups, => k(k-1)/=15) group means, standard errors ... for all permutations at once
##	1	-1	0	0	0	0	0...
##	1	0	-1	0	0	0	0...
##	...
##	0	1	-1	0	0	0	0...
## group.mean.matrix %*% contrast.matrix (Lmat)	=> 10000*15 differences of group means of type xbar-ybar
## repeat this m genes times
## so we have a loop over genes with a lot of dot products in
## to avoid dynamic memory problems, the process should be splitted by arrays with more than 10000 genes 
## calling the above routine 5...10 times and merge the results together - block wise processing of genes
## this is performed by PTEST.batch which is called by the function dialog.PTEST (split option =5)


if(permutations > 1){
	if(!plot.it){ par(mfrow=c(2,1)) }else{ par(mfrow=c(2,2)) }
	## start t-test   permutations
#		if(is.vector(x)){ m <- 1 }else{ m <- dim(x)[1] }				# number of genes
#		m <- nr
#		if(is.vector(x)){ n <- 1:length(x) }else{ n <- dim(x)[2] }		# number of variable group labels       n-nc new defined ?????? x as vec ?????
	if(!diagnostic){		# trace the elapsed time
		start.time <- date()
		cat("\nstart.time: ", start.time)
		plot(x=1, y=1, type="n", xlim=c(0,200), ylim=c(0,round(nr/1000)),
				xlab=paste("Multiple t-permutation test on",nr,"genes (",permutations," random samples,seed=",seed,")"), ylab="Gene" )
#			sub <- paste("\nData:",xname)
		title(sub=paste("\nPermutation start:", start.time))	
	}
	cat("\nStarting permutations test with ",permutations,"random samples drawn from ", nc," expression values for each gene\n\n")
	pboot <- matrix(0, nrow=nr, ncol=dim(t0)[2])					# initialize matrix tvalues.1 >= statistic comparisons
	categories.1 <- categories										# sample the lables
	cat("\n categories ",categories.1)
	categories.1 <- as.factor(categories.1)						# group labels
	group.label <- unique(categories.1)										# safe the category names as col labels
	
	ni <- as.vector(table(match(categories.1,group.label)))
	nmat <- matrix(rep(ni,permutations),nrow=permutations,byrow=T)	# matrix with ni
	
	Lmat <- contrast.matrix(categories.1, comparison)
	
	set.seed(seed)		# initialize with a 'seed' to establish a reproducable environment for the permutation
	
	for(k in 1:nr){		# nr genes with t0 vectors
		if(!diagnostic){ points(k%%1000, floor(k/1000), col=4+k%%10, pch=15) }		# elapsed time point
		
		xp <- matrix(rep(t(x[k,]), permutations), nrow=permutations, byrow=T)		#template sampling matrix
		xperm <- t(apply(xp, 1, sample))		#sampling
#			print(summary(xperm))
		within.vars <- NULL	
		within.means <- NULL
		tvalues.1 <- NULL
		
		for(i in 1:length(group.label)){
			m1 <- !is.na(match(categories.1,group.label[i]))					# pattern match between groups and experiment names
			y <- xperm[,m1]												# take only xcols which correspond to a unique category
			within.means <- cbind(within.means,rowMeans(y))			# means within groups
			within.vars <- cbind(within.vars,rowVars(y))				# variances		
		}
		
		estimate.1 <- as.matrix(within.means) %*% Lmat
		stderr.1  <- NULL
		
		if(var.equal){
			df <- as.matrix(nmat) %*% abs(Lmat) -2						# n1+n2-2
			stderr.1 <- sqrt( ( as.matrix(within.vars*(nmat-1)) %*% abs(Lmat) )/df * as.matrix(1/nmat) %*% abs(Lmat) )
		}else{
			var.n <- within.vars/nmat										# vars/ni
			stderr.1 <- sqrt(as.matrix(var.n) %*% abs(Lmat))		# Welch method: sqrt( var.x/nx + var.y/ny )
			# df <- (varx/nx+vary/ny)^2/((varx/nx)^2/(nx-1)+(vary/ny)^2/(ny-1))
			df <- (as.matrix(var.n) %*% abs(Lmat))^2/as.matrix((var.n)^2/(nmat-1)) %*% abs(Lmat)
		}
		
		tvalues.1 <- estimate.1/stderr.1
		pval.1 <- 2*(1-pt(abs(tvalues.1), df))
		pvalues.1 <- matrix( pval.1, nrow=dim(estimate.1)[1], ncol=dim(estimate.1)[2], byrow=F, dimnames=list(row.names(x), comparison) )
		
		# update
		if(diagnostic){ pp0 <- matrix(0, nrow=permutations, ncol=dim(t0)[2], dimnames=list(NULL, dimnames(t0)[[2]])) }
		ci95 <- NULL
		diff.1 <- NULL
		p.null <- NULL
		lower.sig <- NULL
		upper.sig <- NULL
		
		n1 <- dim(t0)[2]
		diff.1 <- colMeans(estimate.1)
		
		for(i in 1:dim(t0)[2]){					# i=1...p comparisons(e.g.15) , k=1...m genes , [,i]= all 1...1000 permutations
			p.null[i] <- sum(ifelse(abs(estimate.1[,i]) >= abs(diff.1[i]),1,0))/permutations
			if(maxt){
				tt.boot <- ifelse(abs(tvalues.1[,i]) >= abs(t0[k,i]),1,0)		# max T: compare t-values, se corrected diffs
			}else{	
				tt.boot <- ifelse(pvalues.1[,i] <= p0[k,i],1,0)					# min P: compare p-values min(p <= p0,1)
			}	
			pboot[k,i] <- sum(tt.boot)											# count the false positives (Westfall & Young)
			lower.sig[i] <- as.vector(quantile(estimate.1[,i],0.025))
			upper.sig[i] <- as.vector(quantile(estimate.1[,i],0.975))
			
			if(diagnostic){ pp0[,i] <- cumsum(tt.boot)/(1:length(tt.boot)) }
		}
		
		pvalue.sig <- pboot[k,]	/permutations			# actual permutation p-values for gene k
		if(diagnostic){					# plot diagnostics and statistics from bootstrap process
			sub <- paste("\nFig.",k,". P-values from ",permutations," random samples",sep="")
			sub <- paste(sub," for ",p," specified comparisons.",sep="")
			if(!(missing(gin))){ sub <- paste(sub,"\nGene info:",genome.info[k,]) }
			plot.rows(t(pp0),type="l",gin="",legend=T,cex=0.55,xlab="Samples",ylab="p-values",box=F)
			title(sub=sub,cex=0.8,adj=0)
			boxplot(data.frame(pp0), boxcol=c(12,14), medcol=8, medlwd=1, outline=F, confint=T, style.bxp="old", ylab="p-values                 ",srt=90,cex=cex,adj=1)
			plot.pvalues(pp0)			# number of rejected hyptheses vs. Type I error sum(p <= palpha), palpha 0...1
			
			###### plot t0 -test data ######
			figure <- k
			main <- rn[k]
			aov.pvalue <- round(anova$fp.value[k], 5)
#				mse.sig <- round(anova$within.var[k],5) # se=stderr[k,]
			stat <- data.frame(t(rbind(estimate=diff.1, lower=lower.sig, upper=upper.sig, pvalue=p.null, fc=fold.change[k,], se=stderr[k])))
			row.names(stat) <- dn
			stat$SIG <- factor(cut(stat$pvalue, breaks = c(0, 0.0001, 0.001, 0.01, 0.05, 100), labels = c("****", "***", "**", "*", ""), include.lowest = T))
			
			sub <- paste("Fig.",figure,". Permutation-test (null distribution).",sep="")
			sub <- paste(sub,"\n95 % non-simultaneous confidence intervals",sep="")
			sub <- paste(sub,"\nand p-values for ",p," specified comparisons ",sep="")
			
			if(permutations > 1){ sub <- paste(sub,"\np-values are based on",permutations,"random permutations") }
			if(proc != "noadjust"){ sub <- paste(sub,"\nstep down correction method:",proc) }
			sub.0 <- sub
			stat.0 <- stat
			if(!(missing(gin))){
				lg <- nchar(genome.info[i,])
				sub <- paste(sub,"\ninfo:",substring(genome.info[k,],1,lg))
			}
			par(cex=0.7, adj=1)		#cex changed ek
			title(main=main)
			par(cex=cex)
			
			sub <- paste("Fig.",figure,".Densities of t0 and t-permutation distributions with alpha-quantiles",sep="")
			sub <- paste(sub,"\nt0 is derived from ",p," comparisons, alpha=",palpha)
			sub <- paste(sub,"\nand the t-permutation distribution is based on",permutations," random samples")
			
			tv <- unlist(tvalues.1)
			xlim <- range(tv)
			qt <- quantile(tv,c(palpha/2,1-palpha/2))
			dy1 <- density(tv)
			xlim1 <- range(dy1$x)
#				xlim <- c(-20,20)
			ylim1 <- range(dy1$y)
			ylim <- c(0, max(0.3, ylim1) )
			
			plot.density(x=tv[abs(tv) < qt[2]], normfac=1-palpha/2, col=c(1,16), ylim=1*ylim, xlim=xlim, sub="", xlab="t-values and palpha-quantiles", cex=cex, new=T)
			plot.density(x=statistic[main,], col=c(3,10), new=F)			# density=10,
			plot.density(x=tv[abs(tv) < qt[2]], normfac=1-palpha/2, col=c(1,16), angle=-45, density=10, new=F)
			
			par(adj=0)
			title(sub=sub)
			plot.density(x=tv[tv <= qt[1]], normfac=palpha/2, col=c(5,11), new=F)	# angle=45,density=10,
			plot.density(x=tv[tv >= qt[2]], normfac=palpha/2, col=c(6,13), new=F)	# angle=-45,density=10,
			par(col=5)
			rug(tv[tv <= qt[1]])
			par(col=6)
			rug(tv[tv >= qt[2]])
			par(col=16)
			rug(tv[abs(tv) < qt[2]])
			par(col=3)
			rug(statistic[main,])
			par(col=1)
			abline(v=quantile(tv, c(0.0,0.005,0.025,0.5,0.975,0.995,1.0)), lty=2)
			text(x=quantile(tv, c(0.0,0.005,0.025,0.5,0.975,0.995,1.0)), y=rep(1.2*ylim[2],7), labels=c(0,0.005,0.025,0.5,0.975,0.995,1), srt=90)
			legend(x=7, y=0.9*ylim[2], bty="n",
					legend=c( paste("t0 (",p," comparisons)"),"qt(t,a/2) < t < qt(t,1-a/2)", "t <= qt(t,a/2)", "t >= qt(t,1-a/2)" ),
					density=10,
					angle=c(45,-45,45,-45),
					fill=c(3,1,5,6),
					col=c(3,1,5,6), 
					cex=cex
			)
			####### END of plot t0-test data #######
		}
	}	
	
	
	if(proc!="noadjust"){
		cat("\nAdjusting pvalues using step down procedures, method:",proc,"\n")
		pvalues  <- apply(pboot, 1, adjust.p, proc=proc)		# step down adjustments on rows (not on cols)
	}
	
}

if(!diagnostic){
	end.time <- date()		# plot the elapsed time on the trace plot
	title(sub=paste("\n\n\nend:",end.time))	
}

####### END of sampling process #######