# TODO: Add comment
# 
# Author: E.Korsching Mar 5, 2013
###############################################################################



estimation.sample.size <- function(x, refset=c(1,5,2,3,12,13), testset=c(4,6:11,14:16), method="pearson",
		k=30, ks=5, red.mode="ran", anzSamplings=1000, exp.name="froh.ck", plot1=F, fittype="both", bass=6, subfolder=""){
	# Scale number of experiments down and estimate how close the specific result is to the noise
	# x: data set
	# j: start size of all data
	# i: actual set size
	# k: size of reduction, k will become smaller in the vicinity of zero (see ks : shorter size)
	# l: loop/result counter
	# red.mode: reduction mode: "ran": uniform random numbers to pick lines, "chop": chop off from upper end of matrix
	# anzSamplings: 0 : only original, 1..many: original + bootstrap examples
	# exp.name: a suffix to the x data name, will be numbered (for each shortening)
	# plot1: plot all regression lines
	# fittype: "loess" or "supsmu" (with bass parameter) or "both"
	# bass: smoothing parameter : less smooth - 0..10 - very smooth
	# method: pearson, spearman or euclidean
	# subfolder: e.g. "/results/e1/" will be appended to work folder path
	
	#ini
	k1 <- k		#save k
	j <- nrow(x)	#sample size
	i <- j		#start condition
	rest <- (j %% k)		#modulo
	x.name <- deparse(substitute(x))
	cat("\n x.name",x.name)
	
	#run to estimate number of steps
	l <- 0		#counter start value
	while(i > 2){
		l <- l+1	#loop counter, number of results
		if(i<=rest){ k <- ks }
		i <- i-k	#shortening
	}
	
	#ini2
	erg1 <- data.frame(matrix(0, l, 9))
	names(erg1) <- c("n.samples","real.ssqg","real.sum.ssqg","boot.ssqg.min","boot.ssqg.max","boot.ssqg.median","boot.min.sum.ssqg","boot.max.sum.ssqg","boot.median.sum.ssqg")
	
	#production run
	l <- 0		#renew counter start value
	k <- k1		#renew k
	i <- j		#renew start condition
	while(i > 2){
		l <- l+1	#loop counter, number of results
		if(red.mode=="ran"){
			if(i!=j){
				pick <- runif(n=i, min=1, max=j)
			}else{
				pick <- seq(1,i,1)		#initial situation
			}
		}else{
			pick <- seq(1,i,1)
		}
		
		txt.0 <- paste(x.name,".",exp.name,l, sep="")
		assign(txt.0,
				enumeration.search.dep.boots.f90(x[pick,], anzSamplings=anzSamplings, refset=refset, testset=testset, method=method, sampleonly=F),
				envir=.GlobalEnv)
		txt.1 <- get(txt.0)		#make var
		
		#collect some results
		tmp0 <- anzSamplings+1
		tmp1 <- which.min(txt.1[[2]][2:tmp0,1]) +1		#+1 because we start with 2 and therefore get a relative position !
		tmp2 <- which.max(txt.1[[2]][2:tmp0,1]) +1
		erg1[l,1] <- i										#n.samples
		erg1[l,2] <- txt.1[[2]][1,1]						#real.ssqg
		erg1[l,3] <- txt.1[[2]][1,5]						#real.sum.ssqg
		erg1[l,4] <- txt.1[[2]][tmp1,1]						#boot.ssqg.min
		erg1[l,5] <- txt.1[[2]][tmp2,1]						#boot.ssqg.max
		#Fortran algo might produce some NaN during sampling : na.rm=T
		erg1[l,6] <- median(txt.1[[2]][2:tmp0,1],na.rm=T)	#boot.ssqg.median
		erg1[l,7] <- txt.1[[2]][tmp1,5]						#boot.min.sum.ssqg
		erg1[l,8] <- txt.1[[2]][tmp2,5]						#boot.max.sum.ssqg
		#Fortran algo might produce some NaN during sampling : na.rm=T
		erg1[l,9] <- median(txt.1[[2]][2:tmp0,5],na.rm=T)	#boot.median.sum.ssqg,  approximation here ! check if appropriate for normalisation
		
		if(i<=rest){ k <- ks }
		i <- i-k	#shortening
	}
	
	#plot -- regression line graphs
	if(plot1){
		pdf(file=paste(getwd(),subfolder,x.name,".",exp.name,".quality",".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""),
				width=11, height=7, onefile=T, title=paste(x.name,".",exp.name,"quality", sep=""), pointsize=12)
		for(ii in 1:l){
			txt.0 <- paste(x.name,".",exp.name,ii,sep="")
			txt.1 <- get( txt.0 )	#make var
			tmp1 <- length(testset)
			tmp0 <- nrow(txt.1[[2]])
			# plot original graph
			plot.sample.search.dependencies(x=txt.1[[3]][1:tmp1,,drop=F], o=txt.1[[1]][1,(tmp1+1):2],
					ylim=c(-0.7,0.7), y.at=NULL, ylab="correlation", text=T, description=txt.0, cex=0.8, mfrow=c(2,3) )
			# supplement with bootstrap values
			plot(x=1,y=1, type="n", xlim=c(0,10), ylim=c(0,10), axes=F, xlab="", ylab="")
			text(x=1,y=9, labels=paste("min.   : ", round( min(txt.1[[2]][2:tmp0,1]) ,digits=1), ",\n",
							"max.   : ", round( max(txt.1[[2]][2:tmp0,1]) ,digits=1), ",\n",
							"median : ", round( median(txt.1[[2]][2:tmp0,1]) ,digits=1), sep=""),
					adj=0, cex=1.2)
			title(main="bootstrap ssqg")
		}
		dev.off()
	}
	
	pdf(file=paste(getwd(),subfolder,x.name,".",exp.name,".main.quality",".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""),
			width=11, height=7, onefile=T, title=paste(x.name,".",exp.name,"main.quality", sep=""), pointsize=12)
	tmp0 <- c(0,j)
	tmp1 <- c(min(erg1[,c(2,4:6)]),max(erg1[,c(2,4:6)]))	#possible NA problems - see above
	plot(x=erg1[,1],y=erg1[,2], type="p", pch=19, ylim=tmp1, xlab="sample size", ylab="gobal ssq") #ssqg
	points(x=erg1[,1],y=erg1[,4], pch=2)	#min sampling
	points(x=erg1[,1],y=erg1[,5], pch=6)	#max sampling
	points(x=erg1[,1],y=erg1[,6], pch=5)	#median sampling
	
	if(fittype=="supsmu"){
		lines(supsmu(x=erg1[,1],y=erg1[,2], bass=bass), lty=2, col="blue")
		lines(supsmu(x=erg1[,1],y=erg1[,6], bass=bass), lty=2, col="red")
	}
	if(fittype=="loess"){
		lines(erg1[,1], predict(loess(erg1[,2]~erg1[,1])), lty=2, col="blue")		#lwd=2
		lines(erg1[,1], predict(loess(erg1[,6]~erg1[,1])), lty=2, col="red")
	}
	if(fittype=="both"){
		lines(supsmu(x=erg1[,1],y=erg1[,2], bass=bass), lty=2, col="blue")
		lines(supsmu(x=erg1[,1],y=erg1[,6], bass=bass), lty=2, col="red")
		lines(erg1[,1], predict(loess(erg1[,2]~erg1[,1])), lty=2, col="blue")
		lines(erg1[,1], predict(loess(erg1[,6]~erg1[,1])), lty=2, col="red")
	}
	
	legend(x=(tmp0[2]-((tmp0[2]-tmp0[1])*0.19)), y=(tmp1[2]-((tmp1[2]-tmp1[1])*0.85)),
			legend=c("ssqg","sampling, max ssqg","sampling, median ssqg","sampling, min ssqg"), col=c("black","black","black","black"), pch=c(19,6,5,2), cex=0.7)
	legend(x=(tmp0[2]-((tmp0[2]-tmp0[1])*0.38)), y=(tmp1[2]-((tmp1[2]-tmp1[1])*0.85)),
			legend=c("fit of ssqg","fit of sampling median"), col=c("blue","red"), lty=c(2,2), cex=0.7)
	if(red.mode=="ran"){ tmp2 <- "random pick" }else{ tmp2 <- "chop off from end" }
	title(main=paste("Sample size dependency / linear downsizing (",tmp2,") \n Basic global ssq versus sampling global ssq", sep=""))
	dev.off()
	
	return(erg1)
}


#date()
#gfro <- estimation.sample.size(x=fro.orig.1a, refset=c(1,2,16,17,18,19), testset=c(11,12,13,10,4,20,5,15),
#		k=20, ks=5, red.mode="chop", anzSamplings=1000, exp.name="ttestt", plot1=T, fittype="both", bass=4, subfolder="/results/")
#date()	# 4min