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



#options(width=160)

####  packages to load
## install.packages(c("curl","zoo","random"))
## install.packages(c("gridExtra","rmarkdown"))

base::search()
base::searchpaths()
rlang::search_envs()

## download source package			stats not separately available
# https://svn.r-project.org/R/trunk/src/nmath/snorm.c		in rnorm -> norm_rand() called 

#download.packages(pkgs="name", destdir="/home/korschi/R")

# Or get the url manually
#myPackage <- "ggplot2"
#aaa <- available.packages()
#myPackageUrl <- paste(aaa[ rownames(aaa) == myPackage, "Repository"], "/", myPackage, "_", aaa[ rownames(aaa) == myPackage, "Version"], ".tar.gz")
## [1] "https://cran.rstudio.com/src/contrib/ggplot2_2.2.1.tar.gz"
#
## then download
#download.file(url=myPackageUrl, destfile=paste0("/path/to/my/libs", "/", basename(myPackageUrl)))






####  general purpose functions

##  fetch output in pdf
term2pdf <- function(f.name="out", path="results02pr/"){
	# capture output on terminal and save to pdf
	# many fn calls
	# needs: sink(file="tmptxt.txt") before first fn call
	require(rmarkdown)
	fname <- paste(f.name,".Rmd",sep="")
	sink()
	tmp.text <- readLines("tmptxt.txt")
	cat("---\n","\\header-includes: |\n","    \\usepackage{ascii}\n","    \\usepackage[T1]{fontenc}\n","    \\fontsize{10pt}\\selectfont\n","...\n\n",
			sep="", file=fname )
	cat("\\large ", sep="  \n", file=fname, append=T)
	cat(tmp.text, sep="  \n", file=fname, append=T)
	render(fname, output_format="pdf_document", output_file=paste(path,f.name,".pdf",sep=""))
	# cleanup
	file.remove(fname,"tmptxt.txt")
	return()	# join associated pdf files
}


termout2pdf <- function(complete.call, f.name="out", path="results02pr/"){
	# capture output on terminal and save to pdf
	# one fn call
	# complete.call: existing function with all specific parameters
	require(rmarkdown)
	fname <- paste(f.name,".Rmd",sep="")
	sink(file="tmptxt.txt")
	a <- eval(parse(text=complete.call))
	sink()
	tmp.text <- readLines("tmptxt.txt")
	cat("---\n","\\header-includes: |\n","    \\usepackage{ascii}\n","    \\usepackage[T1]{fontenc}\n","    \\fontsize{10pt}\\selectfont\n","...\n\n",
			sep="", file=fname )
	cat("\\large ", sep="  \n", file=fname, append=T)
	cat(tmp.text, sep="  \n", file=fname, append=T)
	render(fname, output_format="pdf_document", output_file=paste(path,f.name,".pdf",sep=""))
	# cleanup
	file.remove(fname,"tmptxt.txt")
	return(a)	# join associated pdf files
}


##  find peaks

# counts
find.peak <- function(x, lowlim=0, v=F){
	# find peaks in a sequence of numbers -- with some additional work for grouping flanks around peaks
	f1 <- function(x, lowlim){# threshold data
		x[x<=lowlim] <- 0
		return(x)
	}
	f2 <- function(dx, xlen){# flank(s) of peak
		dx <- dx[2:xlen]
		# group vector
		sx <- vector("integer",(xlen-1))
		k <- 0	# labels
		a <- 0; b <- 0; c <- 0		# activity
		for(i in 1:(xlen-1)){
			if(dx[i]<0){
				if(a){
					sx[i] <- k
				}else{
					k <- k+1
					a<-1; b<-0; c<-0
					sx[i] <- k
				}
			}else if(dx[i]==0){
				if(b){
					sx[i] <- k
				}else{
					k <- k+1
					a<-0; b<-1; c<-0
					sx[i] <- k
				}
			}else if(dx[i]>0){
				if(c){
					sx[i] <- k
				}else{
					k <- k+1
					a<-0; b<-0; c<-1
					sx[i] <- k
				}
			}
		}
		# join flanks if possible
		msx <- max(sx)
		j <- T; i <- 1
		while(j){
			k <- 0
			if(sum(dx[sx==i])>0 & sum(dx[sx==(i+1)])<0){
				k <- i+1
				sx[sx==i] <- k
			}else if(sum(dx[sx==i])>0 & sum(dx[sx==(i+1)])==0 & sum(dx[sx==(i+2)])<0){
				k <- i+2
				sx[sx==i] <- k
				sx[sx==(i+1)] <- k
			}
			if(k==0){ i <- i+1 }else{ i <- k }
			if(i>=msx){ j <- F }
		}
		sx <- c(sx,sx[length(sx)])	# duplicate last group to correspond to x input
		return(sx)
	}
	if(lowlim>0){ x <- f1(x, lowlim) }
	# duplicate first and last number
	xlen <- length(x)
	x <- c(x[1],x,x[xlen])
	dx <- diff(x)
	sx <- f2(dx, xlen)
	sg <- sign(dx)
	dsg <- diff(sg)
	erg <- which(dsg < 0)
	if(v){
		cat("\nx ",x)
		cat("\ndx",dx)
		cat("\nsg",sg)
		cat("\ndsg",dsg)
		cat("\ngroups",max(sx),"\n")
	}
	return( list(peak.pos=erg, groups=sx, gnumber=length(unique(sx))) )
}
#find.peak(c(9,5,1,2,3,4,3,2,3), v=T)
#find.peak(c(9,5,1,2,3,4,3,2,3), lowlim=3, v=T)
#find.peak(c(0,0,1,2,3,4,3,2,2), v=T)
#find.peak(c(0,2,4,5,5,5,3,2,2), v=T)


# by stas_g
find.peaks <- function (x, m=3){
	# peaks - by name stas_g
	shape <- diff(sign(diff(x, na.pad = FALSE)))
	pks <- sapply(which(shape < 0), FUN = function(i){
				z <- i - m + 1
				z <- ifelse(z > 0, z, 1)
				w <- i + m + 1
				w <- ifelse(w < length(x), w, length(x))
				if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
			})
	pks <- unlist(pks)
	pks		# index
}
# find.peaks(x=c(1,2,3,4,3,2,3,4,4,3,2,1), m=3)
#plot(sin(seq(1,10,.1)))		# index plot
#a <- find.peaks(x=sin(seq(1,10,.1)), m=3)
#points(a, sin(seq(1,10,.1))[a], col="red", pch=16)


# by izmirlig
find.peaks2 <- function (x, thresh=0.05, span=0.25, lspan=0.05, noisey=TRUE){
	# peaks - by name izmirlig
	n <- length(x)
	y <- x
	mu.y.loc <- y
	if(noisey)
	{
		mu.y.loc <- (x[1:(n-2)] + x[2:(n-1)] + x[3:n])/3
		mu.y.loc <- c(mu.y.loc[1], mu.y.loc, mu.y.loc[n-2])
	}
	y.loess <- loess(x~I(1:n), span=span)
	y <- y.loess[[2]]
	sig.y <- var(y.loess$resid, na.rm=TRUE)^0.5
	DX.1 <- sign(diff(mu.y.loc, na.pad = FALSE))
	pks <- which(diff(DX.1, na.pad = FALSE) < 0 & DX.1[-(n-1)] > 0) + 1
	out <- pks
	if(noisey)
	{
		n.w <- floor(lspan*n/2)
		out <- NULL
		for(pk in pks)
		{
			inner <- (pk-n.w):(pk+n.w)
			outer <- c((pk-2*n.w):(pk-n.w),(pk+2*n.w):(pk+n.w))
			mu.y.outer <- mean(y[outer])
			if(!is.na(mu.y.outer)) 
				if (mean(y[inner])-mu.y.outer > thresh*sig.y) out <- c(out, pk)
		}
	}
	out		# index
}
# find.peaks2(x=c(1,2,3,4,3,2,3,4,4,3,2,1), thresh=0.05, span=0.25, lspan=0.05, noisey=F)
#plot(sin(seq(1,10,.1)))		# index plot
#a <- find.peaks2(x=sin(seq(1,10,.1)), thresh=0.05, span=0.25, lspan=0.05, noisey=F)
#points(a, sin(seq(1,10,.1))[a], col="red", pch=16)


# by huber
find.peaks3 <- function(x, y, w=1, span) {
	# peaks - by name huber
	# tune the original curve and find peaks
	# return x_peak positions, index_peak positions, all smoothed y values [n]
	#   y.max [n-2*w], delta [n-2*w]
	require(zoo)
	n <- length(x)
	y.smooth <- loess(y ~ x, span=span)$fitted
	y.max <- rollapply(zoo(y.smooth), 2*w+1, max, align="center")
	delta <- y.max - y.smooth[-c(1:w, n+1-1:w)]
	i.max <- which(delta <= 0) + w
	list(x=x[i.max], i=i.max, y.hat=y.smooth)
}
#a <- find.peaks3(x=x1, y=y1, w=3, span=0.05)


test.find.peaks3 <- function(x, y, w, span) {
	# test
	peaks <- find.peaks3(x=x, y=y, w=w, span=span)
	
	plot(x, y, cex=0.75, col="gray", main=paste("w = ",w,", span = ",span,sep=""))	# original curve
	lines(x, peaks$y.hat, lwd=2)	# smoothed curve
	y.min <- min(y)		# lower base line
	sapply(peaks$i, function(i){ lines(c(x[i],x[i]), c(y.min, peaks$y.hat[i]), col="Red", lty=2) } )	# vertical peak lines
	points(x[peaks$i], peaks$y.hat[peaks$i], col="Red", pch=19, cex=1.25)		# cicles on top of peak lines
}

# test data
#x1 <- 1:1000 / 100 - 5		# -5..5 step 0.01 = 1000 x values
#y1 <- exp(abs(x1)/20) * sin(2*x1 + (x1/5)^2) + cos(10*x1)/5 + rnorm(length(x1), sd=0.05)	# 1000 y values
#
#par(mfrow=c(3,2))
#test.find.peaks3(x1,y1,  2, 0.05)
#test.find.peaks3(x1,y1, 30, 0.05)
#test.find.peaks3(x1,y1,  2, 0.2 )
#test.find.peaks3(x1,y1,  2, 0.4 )
#test.find.peaks3(x1,y1,  2, 0.6 )

#par(mfrow=c(1,1))
#a <- find.peaks3(x=seq(1,10,.1), y=sin(seq(1,10,.1)), w=1, span=0.05)
#plot(seq(1,10,.1), sin(seq(1,10,.1)))
#lines(seq(1,10,.1), a$y.hat, lwd=2, col="green")
#points(a$x, sin(seq(1,10,.1))[a$x], col="red", pch=16)






####  diagnostics on R random number types
# and finally create data sets


# old usage -- but bi-modal peak in z score ... wrong
#  seems to be an wrong combination of random(rnd) functions or internal dependencies ...
#  --> wrong thinking - too many ideas at once
#...
#for(i in 1:samplN){
#	mean.te <- runif(a, b1, b2)		# runif(n, min, max)
#	sd.te <- runif(a, c1, c2)
#	tmp <- rnorm(a, mean.te, sd.te)		# rnorm(n, mean, sd)
#	erg[i,] <- tmp
#	erg1[i,] <- (tmp - mean(tmp))/fnO1(x=tmp, p=p[1])	# z
#	...
#}
#...


## ... solve the basic problem


## r-unif-natural -> r-norm-natural by Box-Muller
bmt.test <- function(n, ymax=1){
	# n : 2*number of inputs and results
	source("../0functions/0general/hist.plot.color.R")
	#	samples <- matrix(ncol=2,nrow=n)
	u1 <- runif(n)
	u2 <- runif(n)
	R <- sqrt(-2*log(1-u1))
	theta <- 2*pi*u2
	X <- R*cos(theta)
	Y <- R*sin(theta)
	
	norm01 <- rnorm(n)
	df1 <- cbind(u1=u1, u2=u2, x=X, y=Y, norm=norm01)
	par(mfrow=c(3,2))
	hist.plot.color(df1[,1], bin.num=20, lty=3, lwd=1, xlab="u1 runif", ylab="norm=1", x.at=c(0,1), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(0,1), x.range=c(0,1), y.max=0.1, col.f="chartreuse1", col.b="chartreuse4", cex=1, h.title="", add=F)
	hist.plot.color(df1[,2], bin.num=20, lty=3, lwd=1, xlab="u2 runif", ylab="norm=1", x.at=c(0,1), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(0,1), x.range=c(0,1), y.max=0.1, col.f="chartreuse1", col.b="chartreuse4", cex=1, h.title="", add=F)
	
	hist.plot.color(df1[,5], bin.num=20, lty=3, lwd=1, xlab="rnorm,Y,X", ylab="norm=1", x.at=c(-3,-2,-1,0,1,2,3), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(-3,3), x.range=c(-4,4), y.max=ymax, col.f="cadetblue2", col.b="blue", cex=1, h.title="", add=F)
	hist.plot.color(df1[,4], bin.num=20, lty=3, lwd=1, xlab="", ylab="", x.at=c(-3,-2,-1,0,1,2,3), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(-3,3), x.range=c(-4,4), y.max=ymax, col.f="gold1", col.b="red", cex=1, h.title="", add=T)
	hist.plot.color(df1[,3], bin.num=20, lty=3, lwd=1, xlab="", ylab="", x.at=c(-3,-2,-1,0,1,2,3), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(-3,3), x.range=c(-4,4), y.max=ymax, col.f="navajowhite2", col.b="darkorange2", cex=1, h.title="", add=T)
	
	hist.plot.color(df1[,5], bin.num=20, lty=3, lwd=1, xlab="rnorm", ylab="norm=1", x.at=c(-3,-2,-1,0,1,2,3), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(-3,3), x.range=c(-4,4), y.max=ymax, col.f="cadetblue2", col.b="blue", cex=1, h.title="", add=F)
	hist.plot.color(df1[,4], bin.num=20, lty=3, lwd=1, xlab="Y", ylab="norm=1", x.at=c(-3,-2,-1,0,1,2,3), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(-3,3), x.range=c(-4,4), y.max=ymax, col.f="gold1", col.b="red", cex=1, h.title="", add=F)
	hist.plot.color(df1[,3], bin.num=20, lty=3, lwd=1, xlab="X", ylab="norm=1", x.at=c(-3,-2,-1,0,1,2,3), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(-3,3), x.range=c(-4,4), y.max=ymax, col.f="navajowhite2", col.b="darkorange2", cex=1, h.title="", add=F)
	par(mfrow=c(1,1))
}
bmt.test(n=100, ymax=0.25)
bmt.test(n=1000, ymax=0.25)
bmt.test(n=10000, ymax=0.25)

bmt.ext <- function(n, mu=0, sd=1){
	# Box-Muller transform
	# n : vector: even number of unif data items expected
	xlen <- length(n)
	if((xlen %% 2)!=0){ stop("n should be even") }
	ulen <- xlen/2
	u1 <- n[1:ulen]
	u2 <- n[(ulen+1):xlen]
	R <- sqrt(-2*log(1-u1))
	theta <- 2*pi*u2
	X <- sd* R*cos(theta) +mu
	Y <- sd* R*sin(theta) +mu
	erg <- c(X,Y)
	return(erg)
}

#a <- bmt.ext(n=runif(8), mu=c(1,2,3,4,5,6,7,8,9,10), sd=c(1,1,1,1,1,1,1,1,1,1))			# 0.61113897 0.16230650 0.32680038 0.59328912 ...
#summary(a)
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# -2.69942 -0.72883  0.12784  0.03156  0.85700  2.44242 

range(aa_runif_n60000_g1_gs20_a)	# 1.000028 299.999938
a <- bmt.ext(n=adaptScale(aa_runif_n60000_g1_gs20_a[1:100,1], minS=1, maxS=300, minT=0, maxT=1))		# 4.878527 282.890638  57.766401  21.619363 ...
summary(a)
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# -2.50453 -0.57696  0.03618  0.09217  0.71972  2.43966 


bmt.d.plot <- function(n, ymax=1, title=""){
	# diagnostic plots
	xlen <- length(n)
	if((xlen %% 2)!=0){ stop("n should be even") }
	ulen <- xlen/2
	u1 <- n[1:ulen]
	u2 <- n[(ulen+1):xlen]
	R <- sqrt(-2*log(1-u1))
	theta <- 2*pi*u2
	X <- R*cos(theta)
	Y <- R*sin(theta)
	
	df1 <- cbind(u1=u1, u2=u2, x=X, y=Y)
	par(mfrow=c(3,2))
	hist.plot.color(df1[,1], bin.num=20, lty=3, lwd=1, xlab="u1", ylab="norm area=1", x.at=c(0,1), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(0,1), x.range=c(0,1), y.max=0.1, col.f="chartreuse1", col.b="chartreuse4", cex=1, h.title="", add=F)
	hist.plot.color(df1[,2], bin.num=20, lty=3, lwd=1, xlab="u2", ylab="norm area=1", x.at=c(0,1), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(0,1), x.range=c(0,1), y.max=0.1, col.f="chartreuse1", col.b="chartreuse4", cex=1, h.title="", add=F)
	
	hist.plot.color(df1[,c(3,4)], bin.num=20, lty=3, lwd=1, xlab="Y,X", ylab="norm area=1", x.at=c(-3,-2,-1,0,1,2,3), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(-3,3), x.range=c(-4,4), y.max=ymax, col.f="gray80", col.b="red", cex=1, h.title="", add=F)
#	hist.plot.color(df1[,3], bin.num=20, lty=3, lwd=1, xlab="", ylab="", x.at=c(-3,-2,-1,0,1,2,3), y.at=NULL, bar.width=0.7, offset=0, digits=2,
#			norm=T, xlim=c(-3,3), x.range=c(-4,4), y.max=ymax, col.f="navajowhite2", col.b="darkorange2", cex=1, h.title="", add=T)
	
	hist.plot.color(df1[,4], bin.num=20, lty=3, lwd=1, xlab="Y", ylab="norm area=1", x.at=c(-3,-2,-1,0,1,2,3), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(-3,3), x.range=c(-4,4), y.max=ymax, col.f="gold1", col.b="red", cex=1, h.title="", add=F)
	hist.plot.color(df1[,3], bin.num=20, lty=3, lwd=1, xlab="X", ylab="norm area=1", x.at=c(-3,-2,-1,0,1,2,3), y.at=NULL, bar.width=0.7, offset=0, digits=2,
			norm=T, xlim=c(-3,3), x.range=c(-4,4), y.max=ymax, col.f="navajowhite2", col.b="darkorange2", cex=1, h.title="", add=F)
	
	plot(1,1, type="n", xlim=c(0,1),ylim=c(0,1),xlab="",ylab="",axes=F)
	text(x=0.1,y=0.9, labels=title, adj=c(0,0.5), col="black")
	par(mfrow=c(1,1))
	
	return()
}

bmt.d.plot(n=runif(100), ymax=1, title="R runif")

bmt.d.plot(n=adaptScale(aa_runif_n60000_g1_gs20_a[1:100,1], minS=1, maxS=300, minT=0, maxT=1), ymax=1, title="other rnd")




## distribution sub-functions
# create data sets

# 1)
# random distribution of pseudo gene expression values
# random data set : runif based
# one group with nrow: gene/row number, ncol: column number
rnd.a1 <- function(nr, nc, min=10, max=300){
	# runif based
	rnd <- matrix(runif(nr*nc, min=min, max=max), nrow=nr, ncol=nc)
	return(rnd)
}
#rnd.a1(nr=3, nc=3, min=10, max=300)

rnd.a2 <- function(nat.unif, nr, nc, min=10, max=300){
	# nat.unif based
	nat.u.len <- length(nat.unif)
	needed <- nr*nc
	if(needed>nat.u.len){ stop("nat.unif too small - needed: ",needed," input: ",nat.u.len) }
	rnd <- matrix(adaptScale.E(nat.unif, minT=min, maxT=max), nrow=nr, ncol=nc)
	return(rnd)
}
#rnd.a2(nat.unif=a[1:9], nr=3, nc=3, min=10, max=300)


# 2)
rnd.b <- function(n, g, gs, min1=10, max1=300, min2=0.2, max2=2.5){
	# random data set : rnorm(mean,sd) based,  mean,sd: runif  based
	# random distribution of gene expression values, one to several groups
	# n: gene number, g: group number, gs: group size (columns)
	# min1,max1: mean range,  min2,max2: sd range
	rnd <- matrix(0,n,g*gs)
	for(i in 1:g){	# group
		mean.runif <- runif(n, min1, max1)
		#cat("mean ",mean.runif,"\n")
		sd.runif <- runif(n, min2, max2)
		#cat("sd ",sd.runif,"\n")
		for(j in 1:n){	# genes, row wise for all group members
			rnd[j, ((i-1)*gs+1):(i*gs)] <- rnorm(gs, mean=mean.runif[j], sd=sd.runif[j])
		}
	}
	# check&correct
	a <- sum(rnd<0)
	a1 <- n*g*gs
	cat("% <0 :",a/a1," numbers: ",a," of ",a1,"\n")
	rnd <- abs(rnd)
	return(rnd)
}
#rnd.b(n=3, g=3, gs=3, min1=10, max1=300, min2=0.2, max2=2.5)

#g <- sample(g)		# additional shuffling (?)
# sparse data?


rnd.c <- function(nat.unif, n, g, gs, min1=10, max1=300, min2=0.2, max2=2.5){
	# random data set : norm box-muller unif based,  mean,sd: unif  based
	# random distribution of gene expression values, one to several groups
	# n: gene number, g: group number, gs: group size (columns: needs to be even->bmt)
	# nat.unif: vector length: mean:sd:event  n*g + n*g + n*g*gs
	# min1, max1: mean range,  min2, max2: sd range
	# bmt needs an even group size 2,4,..
	nat.u.len <- length(nat.unif)
	needed <- n*g + n*g + n*g*gs
	if(needed>nat.u.len){ stop("nat.unif too small - needed: ",needed," input: ",nat.u.len) }
	nu.mean <- adaptScale.E(nat.unif[1:(n*g)], minT=min1, maxT=max1)
	nu.sd <- adaptScale.E(nat.unif[(n*g+1):(n*g*2)], minT=min2, maxT=max2)
	nu.event <- adaptScale.E(nat.unif[(n*g*2+1):(n*g*2+n*g*gs)], minT=0, maxT=1)
	# take one runif sample and split
	rnd <- matrix(0,n,g*gs)
	for(i in 1:g){	# per group
		# per gene/row
		mean.runif <- nu.mean[((i-1)*n+1):(i*n)]
		#cat("mean ",mean.runif,"\n")
		sd.runif <- nu.sd[((i-1)*n+1):(i*n)]
		#cat("sd ",sd.runif,"\n")
		for(j in 1:n){	# genes, row wise for all group members
			tmp <- nu.event[(((i-1)*n+j)*gs-(gs-1)):(((i-1)*n+j)*gs)]
			#cat("tmp ",tmp," mean ",mean.runif[j]," sd ",sd=sd.runif[j],"\n")
			rnd[j, ((i-1)*gs+1):(i*gs)] <- bmt.ext(tmp, mu=mean.runif[j], sd=sd.runif[j])
		}
	}
	# check&correct
	a <- sum(rnd<0)
	a1 <- n*g*gs
	cat("% <0 :",a/a1," numbers: ",a," of ",a1,"\n")
	rnd <- abs(rnd)
	return(rnd)
}

#a <- readQRNG(file="/home/korschi/on1/rndberlin/0-100MB_a.bin", n=10, byte=4)
#range(a)	# -2147483451  2147483559 length 26,214,400
#hist(a); plot(density(a))

#a1 <- rnd.c(a, n=3, g=3, gs=4, min1=10, max1=300, min2=0.2, max2=2.5)

#sd(a1[1,1:4]); sd(a1[1,5:8]); sd(a1[1,9:12])
#abs( matrix(c(1,2,-3,4,-5,6,7,8,9),3,3) )

#par(mfrow=c(2,2));
#hist.plot.color(rnorm(1000,0,0.3), bin.num=20, lty=3, lwd=1, xlab="", ylab="", x.at=c(0,1), y.at=NULL, bar.width=0.7, offset=0, digits=2,
#		norm=T, xlim=c(-4,4), x.range=NULL, y.max=0.2, col.f="chartreuse1", col.b="chartreuse4", cex=1, h.title="", add=F)
#hist.plot.color(rnorm(1000,0,0.5), bin.num=20, lty=3, lwd=1, xlab="", ylab="", x.at=c(0,1), y.at=NULL, bar.width=0.7, offset=0, digits=2,
#		norm=T, xlim=c(-4,4), x.range=NULL, y.max=0.2, col.f="chartreuse1", col.b="chartreuse4", cex=1, h.title="", add=F)
#hist.plot.color(rnorm(1000,0,1), bin.num=20, lty=3, lwd=1, xlab="", ylab="", x.at=c(0,1), y.at=NULL, bar.width=0.7, offset=0, digits=2,
#		norm=T, xlim=c(-4,4), x.range=NULL, y.max=0.2, col.f="chartreuse1", col.b="chartreuse4", cex=1, h.title="", add=F)
#hist.plot.color(rnorm(1000,0,1.5), bin.num=20, lty=3, lwd=1, xlab="", ylab="", x.at=c(0,1), y.at=NULL, bar.width=0.7, offset=0, digits=2,
#		norm=T, xlim=c(-4,4), x.range=NULL, y.max=0.2, col.f="chartreuse1", col.b="chartreuse4", cex=1, h.title="", add=F)
#par(mfrow=c(1,1));



# 3)
rnd.d <- function(x, min=0.5, max=1.5, sd.one=T){
	# random data set : by a sample of a real data set
	#   uniform random numbers by computer based generator
	#   finally normal distributed noise on the template sample
	# x: real count matrix (if present 22*gene: 1370468, 22*transcript: 5386150)
	# rnorm: mean from x, sd from runif
	# min,max: sd range for rnd number generator
	# sd.one: T: one sd value for all, F: for each another sd value
	# rnd: sampled version of x
	nr <- nrow(x)
	nc <- ncol(x)
	x1 <- as.vector(x)	# unfold by column
	x1len <- length(x1)
	
	if(sd.one){
		r.sd <- runif(1, min, max)
	}else{
		r.sd <- runif(x1len, min, max)
	}
	y <- rnorm(x1len, mean=x1, sd=r.sd)
	
	rnd <- matrix(y, nrow=nr, ncol=nc, byrow=F)	# fold again by column
	dimnames(rnd)[[1]] <- dimnames(x)[[1]]
	dimnames(rnd)[[2]] <- dimnames(x)[[2]]
	# check&correct
	a <- sum(rnd<0)
	a1 <- nr*nc
	cat("% <0 :",a/a1," numbers: ",a," of ",a1,"\n")
	rnd <- abs(rnd)
	return(rnd)
}
#rnd.d(matrix(c(1,8,7,5,3,9,2,5,4),3,3), min=0.5, max=1.5, sd.one=T)
#rnd.d(matrix(c(1,8,7,5,3,9,2,5,4),3,3), min=0.5, max=1.5, sd.one=F)

rnorm(5, mean=c(1,2,3,4,5), sd=c(0.3,0.6,1,1.3,1.6))
rnorm(5, mean=c(1,2,3,4,5), sd=0.3)


test.rnd.d <- function(x, i=1000){
	# test
	print(x)
	z <- x
	for(j in 1:i){
		y <- rnd.d(x)
		z <- z+y
	}
	z <- z/(i+1)
	print(z)
}
#test.rnd.d(matrix(c(1,8,7,5,3,9,2,5,4),3,3), i=1000)



# R  ?runif  ?.Random.seed

rnd.e <- function(x, nat.unif, min=0.5, max=1.5, sd.one=T){
	# random data set : by a sample of a real data set
	#   uniform random numbers from atmospheric noise/quantum noise
	#   (.Machine$integer.max)
	# x: real data/count matrix (if present 22*gene: 1370468, 22*transcript: 5386150)
	#  x1 needs to be even -> bmt
	# nat.unif: vector length: sd:event  1 or x1len + x1len
	# min, max: sd range
	# sd.one: T: one sd value for all, F: for each another sd value
	# bmt: mean from x, sd and event from external rnd source
	# rnd: sampled version of x
	nr <- nrow(x)
	nc <- ncol(x)
	x1 <- as.vector(x)	# unfold by column
	x1len <- length(x1)
	if((x1len %% 2)!=0){ # should be even
		tmp <- T
		x1 <- c(0,x1)
		x1len <- x1len + 1
	}
	nat.u.len <- length(nat.unif)
	if(sd.one){
		needed <- 1 + x1len
		if(needed>nat.u.len){ stop("nat.unif too small - needed: ",needed," input: ",nat.u.len) }
		nu.sd <- nat.unif[1]
		nu.event <- nat.unif[2:(x1len+1)]
	}else{
		needed <- x1len + x1len
		if(needed>nat.u.len){ stop("nat.unif too small - needed: ",needed," input: ",nat.u.len) }
		nu.sd <- nat.unif[1:x1len]
		nu.event <- nat.unif[(x1len+1):(2*x1len)]
	}
	nu.sd <- adaptScale.E(nu.sd, minT=min, maxT=max)
	nu.event <- adaptScale.E(nu.event, minT=0, maxT=1)
	# *2 back take the firsts half-bmt cos
	y <- bmt.ext(nu.event, mu=x1, sd=nu.sd)[1:x1len]
	if(tmp){ y <- y[-1] }
	rnd <- matrix(y, nrow=nr, ncol=nc, byrow=F)	# fold again by column
	dimnames(rnd)[[1]] <- dimnames(x)[[1]]
	dimnames(rnd)[[2]] <- dimnames(x)[[2]]
	# check&correct
	a <- sum(rnd<0)
	a1 <- nr*nc
	cat("% <0 :",a/a1," numbers: ",a," of ",a1,"\n")
	rnd <- abs(rnd)
	return(rnd)
}
#rnd.e(matrix(c(1,8,7,5,3,9,2,5,4),3,3), nat.unif=a[1:11], min=0.5, max=1.5, sd.one=T)		# +1
#rnd.e(matrix(c(1,8,7,5,3,9,2,5,4),3,3), nat.unif=a[1:20], min=0.5, max=1.5, sd.one=F)		# +2

save.image()



## diagnostic
test.rnd.R <- function(n=2, g=1, gs=1, ti.txt="", te.fn="", sdata=F, type="p", file="results03rndtests/01_"){
	# diagnostic for R random numbers & functions
	# n genes, g group(s), gs group size
	# type: "l","p","spl", {3,46, ..}
	# te.fn: one name of a R function accessible in the present session
	fnO1 <- get(te.fn)
	tmp <- fnO1(n, g, gs)
	if(sdata){ assign(paste("aa_",ti.txt,sep=""), value=tmp, envir=.GlobalEnv) }	# save tmp
	
	nc <- ncol(tmp)
	pdf(paste(file,ti.txt,".pdf",sep=""),width=7,height=7)
	par(mfrow=c(2,2))
	for(i in 1:nc){
		hist.plot.p.l.curve(x=tmp[,i], bin.num=40, norm=F, type=type,
				xlab="scale", ylog=F, digits=2, col.f="green", cex=0.7,
				h.title=paste("col# ",i," bin# 40","\n ",ti.txt,sep=""))
	}
	par(mfrow=c(1,1))
	dev.off()
	return()
}
# test.rnd.R(n=60000, g=1, gs=20, ti.txt="rnorm_mean_sd_runif_n60000_g1_gs20", te.fn="rnd.a1", sdata=F)


test.rnd.N <- function(x, rnd=NULL, ti.txt="", te.fn="", sdata=F, type="p", file="results03rndtests/01_"){
	# diagnostic on Natural random numbers & functions
	# x: a real data set (df, matrix)
	# type: "l","p","spl", {3,46, ..}
	# te.fn: one name of a R function accessible in the present session
	# rnd: natural random numbers
	fnO1 <- get(te.fn)
	if(is.null(rnd)){ tmp <- fnO1(x) }else{ tmp <- fnO1(x, rnd) }
	if(sdata){ assign(paste("aa_",ti.txt,sep=""), value=tmp, envir=.GlobalEnv) }	# save tmp
	
	nc <- ncol(tmp)
	pdf(paste(file,ti.txt,".pdf",sep=""),width=7,height=7)
	par(mfrow=c(2,2))
	for(i in 1:nc){
		hist.plot.p.l.curve(x=tmp[,i], bin.num=40, norm=F, type=type,
				xlab="scale", ylog=T, digits=2, col.f="green", cex=0.7,
				h.title=paste("col# ",i," bin# 40","\n ",ti.txt,sep=""))
	}
	par(mfrow=c(1,1))
	dev.off()
	return()
}
# test.rnd.N(x, rnd=??, ti.txt="rnorm_mean_sd_runif_n60000_g1_gs20", te.fn="rnd.a2", sdata=F)


#range(aa_rnorm_mean_sd_runif_n60000_g1_gs20[,1])
#boxplot(aa_rnorm_mean_sd_runif_n60000_g1_gs20)
#hist(aa_rnorm_mean_sd_runif_n60000_g1_gs20[,1], breaks=40)





# further steps -> functions02.R  # sd mad Mad - test approach