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



proximity.bootstrap.sampling <- function(x, reiheRef, reiheTest, anzSamplings=100, method="pearson", output="list"){
	# draw bootstrap / permutation sample and calculate proximity in reference versus test cols
	# cross tabulation by a similarity measure of reference cols versus test cols
	#  every special approach is described next to the if clause
	# x: data.frame : will be processed col wise
	# reiheTest: order of the non-reference cols : rows of the resulting matrix - or selection
	# reiheRef: order of the reference cols : cols of the resulting matrix - or selection
	# anzSamplings: number of samplings (shuffling/bootstrap), 0: only unsampled result will be created
	# method: pearson, spearman or euclidean
	# out: stacked matrix of samples including the original unsampled object on place 1
	
	#ini
	if(!(method %in% c("euclidean","pearson","spearman"))){ cat("\n wrong similarity method "); return() }
	if(anzSamplings<0){ cat("\n anzSamplings has to be 0, 1, 2 and up (integer)"); return() }
	anzSamplings <- anzSamplings+1		#add the (unsampled) proximity table of the original data set - at least that will be produced
	
	nr <- length(reiheTest)
	nc <- length(reiheRef)
	nrx <- nrow(x)
	ncx <- ncol(x)
	if(nr<2 | nc<2 | nrx<2){ cat("\n check: reiheTest >=2 , reiheRef >=2 , data rows >=2"); break }
	
	var.list <- apply(x,2,var)
	var.list <- var.list[order(var.list)]	#calculate all ordered (0 ascending) variances
	cat("\n var ", var.list, "\n")		#output var values
	if(method=="pearson" & sum(var.list==0)>=1){ cat("\n warning -- var=0 -- special handling for 'cor' activated, check results \n") }		#give warning (especially for cor measure)
	
	#functions
	if(method=="euclidean"){
		ct <- function(x,y) apply(y,2,function(x,y) sqrt(sum((x-y)^2)) ,x)
	}
	
	if(method=="pearson"){
		ct <- function(x,y){
			cor.fn <- function(x,y){
				tmp <- cov(x,y)		# means: test if var(x) or var(y) =0 
				if(tmp==0){			# if cov 0 : -> set cor by definition to a value nearly 0 (randomly chosen  of +0.00000001 or -0.00000001 )
					tmp <- sample(x=c(0.00000001,-0.00000001), size=1)
				}else{
					tmp <- cor(x,y)		#calculate regular Pearson cor coefficient
				}
				return(tmp)
			}
			apply(y,2,cor.fn,x)
		}
	}

	if(method=="spearman"){
		ct <- function(x,y){
			cor.fn <- function(x,y){
#				tmp <- cov(x,y)		# means: test if var(x) or var(y) =0 			########## für spearman sinnvoll? #########
#				if(tmp==0){			# if cov 0 : -> set cor by definition to a value nearly 0 (randomly chosen  of +0.00000001 or -0.00000001 )
#					tmp <- sample(x=c(0.00000001,-0.00000001), size=1)
#				}else{
					tmp <- cor(x,y, method="spearman")		#calculate regular Spearman rank cor coefficient
#				}
				return(tmp)
			}
			apply(y,2,cor.fn,x)
		}
	}
	
	
	###work -- output always a matrix
	#def
	z <- matrix(0, nr*(anzSamplings), nc)
	name.x <- names(x)
	dimnames(z)[[2]] <- name.x[reiheRef]
	dimnames(z)[[1]] <- rep(name.x[reiheTest], times=(anzSamplings))
	
	y <- matrix(0, nrx, ncx)
	
	
	# normal proximity table has to be added first , all other tables are samplings
	j <- 1
	z[((j*nr)-(nr-1)):((j+1)*nr-nr),] <- apply(x[,reiheRef,drop=F], 2, ct, x[,reiheTest,drop=F])
	if(anzSamplings==1){ return(z) }
	
	
	#bootstrap sampling
	if(output=="boot.raw"){		# bootstrap sampling on the complete raw data
		
		for(j in 2:(anzSamplings)){
			y <- matrix(sample(as.vector(unlist(x)),replace=T), nrx, ncx)		# basic 0-hypothesis  (with replacement)
			z[((j*nr)-(nr-1)):((j+1)*nr-nr),] <- apply(y[,reiheRef,drop=F], 2, ct, y[,reiheTest,drop=F])
		}
		return(z)
	}
	
	if(output=="boot.prox"){		# bootstrap sampling on the complete proximity table
		
		for(j in 2:(anzSamplings)){
			z[((j*nr)-(nr-1)):((j+1)*nr-nr),] <- sample(z[((1*nr)-(nr-1)):((1+1)*nr-nr),],replace=T)		# sampling on z[j=1,]
		}
		return(z)
	}
	
	if(output=="boot.raw.col"){		# bootstrap sampling on the raw data (per factor col)
		
		for(j in 2:(anzSamplings)){
			for(i in 1:ncx){ y[,i] <- sample(x[,i],replace=T) }
			z[((j*nr)-(nr-1)):((j+1)*nr-nr),] <- apply(y[,reiheRef,drop=F], 2, ct, y[,reiheTest,drop=F])
		}
		return(z)
	}
	
	if(output=="boot.prox.col"){		# bootstrap sampling on the proximity table (per reference col)
		
		for(j in 2:(anzSamplings)){
			for(i in 1:nc){
				z[((j*nr)-(nr-1)):((j+1)*nr-nr),i] <- sample(z[((1*nr)-(nr-1)):((1+1)*nr-nr),i],replace=T)		# sampling on z[j=1,]
			}
		}
		return(z)
	}
	
	
	#permutation sampling
	if(output=="perm.raw"){		# permutation sampling on the complete raw data
		
		for(j in 2:(anzSamplings)){
			y <- matrix(sample(as.vector(unlist(x)),replace=F), nrx, ncx)		# basic 0-hypothesis  (without replacement)  [x is data.frame ergo a list obj]
			z[((j*nr)-(nr-1)):((j+1)*nr-nr),] <- apply(y[,reiheRef,drop=F], 2, ct, y[,reiheTest,drop=F])
		}
		return(z)
	}
	
	if(output=="perm.prox"){		# permutation sampling on the complete proximity table
		
		for(j in 2:(anzSamplings)){
			z[((j*nr)-(nr-1)):((j+1)*nr-nr),] <- sample(z[((1*nr)-(nr-1)):((1+1)*nr-nr),],replace=F)		# sampling on z[j=1,]
		}
		return(z)
	}
	
	if(output=="perm.raw.col"){		# permutation sampling on the raw data (per factor col)
		
		for(j in 2:(anzSamplings)){
			for(i in 1:ncx){ y[,i] <- sample(x[,i],replace=F) }
			z[((j*nr)-(nr-1)):((j+1)*nr-nr),] <- apply(y[,reiheRef,drop=F], 2, ct, y[,reiheTest,drop=F])
		}
		return(z)
	}
	
	if(output=="perm.prox.col"){		# permutation sampling on the proximity table (per reference col)
		
		for(j in 2:(anzSamplings)){
			for(i in 1:nc){
				z[((j*nr)-(nr-1)):((j+1)*nr-nr),i] <- sample(z[((1*nr)-(nr-1)):((1+1)*nr-nr),i],replace=F)		# sampling on z[j=1,]
			}
		}
		return(z)
	}
}


#aac2 <- proximity.bootstrap.sampling(He589.c,reiheRef=c(1,5,3,2,12,13),reiheTest=c(4,6:11,14:16),anzSamplings=4,method="cor",output="matrix")