# TODO: Add comment
# 
# Author: E.Korsching Mar 16, 2012
###############################################################################


rnr <- function(
		x,					# pvalue matrix containing columns with pairwise comparison results (results of statistical tests)
		reg=NULL,			# comparison columns in which we want to see regulation (if given)  -  order in name matter e.g. a.d or d.a
		sorting=T,			# sort reg columns (regulation) before noreg columns (no regulation)
		palpha=0.05,		# alpha threshold
		return.logic=T		# return a vector denoting which rows (genes) remain after applying the threshold rules
){
	# Control of Type I error
	# return.logic: T: a data.frame  with row.names and a logical vec for candidates remaining after application of the rules
	#  F: remaining subset with their p values
	# example call: see bottom
	
	#ini
#	dsname <- deparse(substitute(x))
	nc <-dim(x)[2]
	nr <- dim(x)[1]
	if(nc<2 | nr<2){ cat("\n column and row number have to be greater 2"); break }
	cnames <- dimnames(x)[[2]]				#col names
	states <- rep(0, nc)						#reg noreg states, default setting: noreg
	
	# missing reg - noreg
	if(is.null(reg)){
		cat("\n Warning: All comparisons are set to NO regulation\n")
		#already default setting
	} else{
		m1 <- !is.na(match(cnames, reg))		#index to comparison cols which are defined as: reg
		states[m1] <- 1							# set to 1
		
		if(sorting){
			states.order <- order(states, decreasing=T)		#reg in the beginning
			x <- x[ ,states.order]
			states <- states[states.order]
			cnames <- cnames[states.order]
		}
	}
	
	#
	x1 <- ifelse(x <= palpha,T,F)		# change p-values into logical values (implicit to 0/1) according to alpha
	
	x2 <- apply(x1, 1, function(x1,states){ if(sum(x1==states)==nc){ return(T) }else{ return(F) } }, states=states)		#test where the rules apply
	x2.rules.true.rows <- sum(x2)
	if(x2.rules.true.rows==0){
		p.final <- NA				# make obvious that 0 is not a very small p value
	}else{
		p.final <- x2.rules.true.rows/nr		# rule based p value of the result 
	}
	
	x1.all.true <- sum(x1)				# significant results of the structure
	p.all.true <- x1.all.true/(nr*nc)		# rule free p value of the structure
	
	x3 <- vector(mode="numeric",length=nc)
	for(i in 1:nc){
		if(states[i]==1){
			x3[i] <- sum(x1[,i])
		}else{
			x3[i] <- sum(!x1[,i])
		}
	}
	x3p <- x3/nr				# rule free p value of each comparison
	
	x5 <- vector(mode="numeric",length=nc)
	x4 <- if(states[1]==1){ x1[,1] }else{ !x1[,1] }			# initialisation of procedure
	x5[1] <- sum(x4)											# p value step 1
	for(i in 2:nc){
		x4 <- x4 & ( if(states[i]==1){ x1[,i] }else{ !x1[,i] } )		# T&F -> F , add col by col and see what is remaining
		x5[i] <- sum(x4)
	}
	x5[x5==0] <- NA			# make obvious that 0 is not a very small p value
	x5p <- x5/nr				# rule based p value development by stepwise applying rules
	
	cat("\ncontrol of Type I error  -- comparisons (col):",nc,", factors (row):",nr,", p alpha:", palpha)
	cat("\nrule free p value of the complete matrix: ",round(p.all.true, 4))
	cat("\n n final: ",x2.rules.true.rows,", rule based p value of the result: ",round(p.final, 4))
	cat("\n stepwise development of the rule based p value: ", paste(round(x5p, 4), sep=" ") )
	cat("\nvalues per comparison:")
	cat("\n_________________________________\n")
	print( data.frame( test=cnames, reg=states, remaining=x3, probability=round(x3p, 4) ))
	cat("_________________________________\n")
	cat("\n")
	
	if(return.logic){ x2 }else{ data.frame(x[x2,]) }		#return full logical vector  or  all rows which passed the filtering
}


#test.rnr <- rnr(
#	x=data.frame(a.b=c(0.01394,0.17394,0.17394,0.17394,0.19204),
#			a.c=c(0.41246,0.04246,0.41246,0.41246,0.41246),
#			b.c=c(0.78472,0.78472,0.78472,0.78472,0.78472),
#			row.names=c("a","b","c","d","e")
#			),		#pvalues
#	reg=c("a.b","a.c"),
#	noreg=c("b.c"),
#	sorting=T,
#	palpha=0.05,
#	return.logic=T
#)