# TODO: Add comment
# 
# Author: E.Korsching  11-2016
###############################################################################



shuffling.p.jp1 <- function(x,y, parA=">=0", parB="<=0", parC=">2", n=100, sampling.method="permutation"){
	# calculate a shuffling p value on all rows of a data.frame
	#  shuffling is done column wise
	# x,y: two chips (see nik2ek.R), need unique row.names
	# n: number of shufflings
	# sampling.method: "permutation" (without replacement) , "bootstrap" (with replacement)
	
	nrx <- nrow(x)
	ncx <- ncol(x)
	nx <- row.names(x)
	nry <- nrow(y)
	ncy <- ncol(y)
	ny <- row.names(y)
	if(nrx!=nry){ stop("x and y have different row numbers") }
	if(sum(nx!=ny)!=0){ stop("x and y have different row names") }
	
	if(sampling.method=="shuffling"){ replace=F }else{ replace=T }
	
	# result df
	z1 <- data.frame(matrix(0,nrx,3), row.names=nx)
	names(z1) <- c("original.hits","sampling.hits","p.value")
	# ratio of overlap per shuffling with original hits
	rn <- vector("numeric",length=n)
	
	# calculate real values
	z <- calculateModel(c1=x, c2=y, parA=parA, parB=parB, parC=parC)	# return df with hits
	zn <- row.names(z)
	zn.l <- length(zn)
	z1[row.names(z), 1] <- 1
	# sampling
	for(i in 1:n){
		z3 <- x
		z4 <- y
		for(j in 1:ncx){
			z3[,j] <- sample(x=z3[,j], replace=replace)
			z4[,j] <- sample(x=z4[,j], replace=replace)
		}
		z <- calculateModel(c1=z3, c2=z4, parA=parA, parB=parB, parC=parC)	# return df with hits
		rn[i] <- sum( row.names(z) %in% zn ) / zn.l		# logical vector
		z1[row.names(z), 2] <- z1[row.names(z), 2] + 1
	}
	
	# calculate hit p
	z1[,3] <- z1[,2] / n
	# order
	z1 <- z1[ order(z1[,1],z1[,3], decreasing=T) , ]
	
	# calculate overlap features
	z5 <- vector("numeric",length=5)
	names(z5) <- c("min.overlap","max.overlap","mean.overlap","sd.overlap","number.hits.original")
	z5[1] <- min(rn)
	z5[2] <- max(rn)
	z5[3] <- mean(rn)
	z5[4] <- sd(rn)
	z5[5] <- zn.l
	
	return( list(dist.hits=z1, dist.overlap=rn, overlap.p=z5) )
}


#a1 <- date()
#a2 <- shuffling.p.jp1(x=chip03,y=chip04, parA=">=0", parB="<=0", parC=">2", n=50000, sampling.method="permutation")
#a1;date()	# 100 : 1s