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



#### functions


# fit with Bezier curve
bezierCurve <- function(x, y, n=10)
{
	# x, y: the x and y coordinates of the hull points
	# n: the number of points in the curve.
	outx <- NULL
	outy <- NULL
	
	i <- 1
	for (t in seq(0, 1, length.out=n))
	{
		b <- bez(x, y, t)
		outx[i] <- b$x
		outy[i] <- b$y
		
		i <- i+1
	}
	
	return (list(x=outx, y=outy))
}

bez <- function(x, y, t)
{
	outx <- 0
	outy <- 0
	n <- length(x)-1
	for (i in 0:n)
	{
		outx <- outx + choose(n, i)*((1-t)^(n-i))*t^i*x[i+1]
		outy <- outy + choose(n, i)*((1-t)^(n-i))*t^i*y[i+1]
	}
	
	return (list(x=outx, y=outy))
}

# Bezier usage and all others
test.fit <- function(ad){
	plot(ad$ax, ad$ay, "b", pch=20, col="gray")
	lines(bezierCurve(ad$ax, ad$ay, 40), col="red")
	
	loessFit <- loess(ay~ax, ad, span=0.6)
	lines(predict(loessFit), col="green")
	
	asS <- smooth.spline(x=ad$ax, y=ad$ay, spar=0.2)
	lines(asS, col="black", lwd=2, lty="14")
	approxFit <- approx(ad, n=15)
	lines(approxFit, col="purple3")
	lowessFit <- data.frame(lowess(ad, f=.6, iter=1))
	lines(lowessFit, col="orange3", lwd=2, lty="14")
}

#test.fit(ad)



# plot - boxplots and density curves
multi.box.densi.panel <- function(x, type="b", iqr=c(0.25,0.75), bin.num, logscale="", flags=c(1,1,1,1,1,1),
		cex=0.6, adj=1, col.b="black", col.f="blue", col.m=c("gray","green"), line=0, prefix=""){
	# plot a panel of box plots or density plots
	# x: data.frame,   type: "b" boxplot, "d" density plot,   logscale: "y" y is log scale,   spar=NULL possible
	source("../0functions/0general/bin.R")
	source("../0functions/0general/boxplot.ek.R")
	ch.df <- is.data.frame(x)
	if(!ch.df){ stop("\nx needs to be a data frame") }
	nc <- ncol(x)
	cnames <- names(x)
	x.name=deparse(substitute(x))
	mar1 <- c(8,5,1,1)
	if(type=="b"){
		ta <- T; tb <- F
	}else if(type=="d"){
		ta <- F; tb <- T
	}else if(type=="bd"){
		ta <- T; tb <- T
	}else{ stop("\ncheck type") }
	
	#layout
	if(ta==T & tb==T){
		layout( matrix(c(1,2), nrow=2, ncol=1, byrow=T), widths=c(34), heights=c(13,13), respect=T)		# relative proportions
	}else{
		layout( matrix(1, nrow=1, ncol=1, byrow=T), widths=c(34), heights=c(13), respect=T)
	}
	
	cex.sav <- par()$cex
	adj.sav <- par()$adj
	par(cex=cex,adj=adj)
	par(mar=mar1)
	
	if(logscale=="y"){
		xx <- x+1			# +1: to get the log scale running if '0'
	}else{
		xx <- x
	}
	if(ta){
		#bs <- boxplot(xx,plot=F)$stats
		boxplot.flex(x=xx,
				iqr=iqr,
				flags=flags,
				at=NULL,
				text=cnames,
				ylab=if(logscale=="y"){"intensity (+1)"}else{"intensity"},
				width=0.7,
				border=col.b,
				col=col.f,
				col.sp=col.m,
				lwd=1, cex=0.7, tcex=1, adj=1, srt=90)
		
		#axis( side=1, at=c(1:nc), labels=cnames, tick=T, las=2 )
		tmp <- paste(prefix," - ",x.name)
		mtext(text=tmp, side=3, line=line, outer=F, adj=1)
	}
	if(tb){
		n <- bin.num
		mat <- matrix(0,n,nc)
		for(i in 1:nc){
			a <- bin(x[,i], bin.num=n)
			mat[,i] <- a$erg[,1]
			#cat("\ncount",range(mat[,i]))
		}
		mmat <- max(mat)
		mat <- mat / mmat
		#return(mat)
		plot(x=0, y=0, type="n", ylab="density [50 bin]", xlab="", xlim=c(0.5,(nc+0.5)), ylim=c(n,-1), axes=F, cex=cex)
		for(i in 1:nc){
			lines(bezierCurve(x=i+(-mat[,i]), y=(1:n), 40), col="red")
		}
		axis( side=1, at=c(1:nc), labels=cnames, tick=T, las=2 )
		axis( side=2, at=c(1,bin.num), labels=c("high  expr.","low expr."), tick=T, las=0, cex.axis=1.1, cex.lab=1.1 )
		tmp <- paste(prefix," - ",x.name)
		mtext(text=tmp, side=3, line=line, outer=F, adj=1)
	}
	# restore
	par(mfrow=c(1,1), mar=c(5,4,4,2), cex=cex.sav, adj=adj.sav)
}



#cat("\nIQR before",apply(mat,2,IQR))
## to avoid smooth.spline()-IQR =0 error, add some small noise in 0 cells
## two col sums before: 1: 1.764706, 15: 1.764706 -- after: 1: 1.764737 15: 1.764746
#ll <- sum(mat==0)
#llj <- jitter(rep(0.000001,ll))
#mat[mat==0] <- llj
#cat("\nIQR  after",apply(mat,2,IQR))

#asS <- smooth.spline(x=i+(-mat[,i]), y=(1:n), spar=spar)		# spar 0.35
#lines(asS, col=col.f[1], lwd=2)

#smooth.spline(x=15+(-a[,15]), y=(1:50), spar=0.2)
#IQR(a[,15])	# 0
#IQR(a[,14])	# 0.05882353
#quantile(a[,15], probs = c(0,0.25,0.5,0.75,1))
#quantile(a[,14], probs = c(0,0.25,0.5,0.75,1))


rm.zero.row <- function(x, row.zero=5){
	# remove all rows <= row.zero counts
	a <- apply(x,1,sum)
	x <- x[!(a<=row.zero), ]
	return(x)
}


# test function normalization
test.iqr.m.norm <- function(x, row.zero=0, out="", iqr=c(0.25,0.75), pdf.n="", width=7, height=11){
	# test function for IQR & median normalization
	# x: data set, columns: samples, rows: genes or else ,   row.zero: NA or del rows with entries smaller than
	# iprange: IQR: 0.25,0.75, 5r95: 0.05,0.95
	# out: symbol for norm output in workspace,   pdf.n: e.g. "out.pdf"
	# every gene will be analyzed on bandwidth across all samples until all samples are encluded
	#  the statistics table will be returned
	source("../0functions/0general/heatmap.ek.R")
	iqrange.fn <- function(x, iqr, na.rm=F, type=7){
		diff(quantile(as.numeric(x), probs=iqr, na.rm=na.rm, names=F, type=type))
	}
	nc <- ncol(x)
	# remove all rows <= row.zero counts
	if(!is.na(row.zero)){ x <- rm.zero.row(x, row.zero=5) }
	nr <- nrow(x)
	#return(x)
	if(pdf.n!=""){
		pdf(pdf.n, width=width, height=height)
		multi.box.densi.panel(x, type="bd", iqr=iqr, flags=c(0,0,1,1,1,1), bin.num=50, prefix=paste("IN raw IQR",iqr[1],"-",iqr[2]))
	}
	# estimate all iqr & median
	colIqr <- apply(x, 2, iqrange.fn, iqr=iqr, na.rm=T, type=7)
	colMea <- apply(x, 2, mean, na.rm=T)	#colMed <- apply(x, 2, median, na.rm=T)
	cat("\ncolIQR",round(colIqr,2),"\ncolMea",round(colMea,2))
	# set target iqr & mean/median
	t.mip <- mean(colIqr)
	t.mea <- mean(colMea)
	cat("\nt.iqr",round(t.mip,2),", t.mea",round(t.mea,2),"\n")
	
	x1 <- data.frame(matrix(0,nr,nc))
	dimnames(x1)[[1]] <- dimnames(x)[[1]]
	dimnames(x1)[[2]] <- dimnames(x)[[2]]
	for(i in 1:nc){
		div.iqr <- colIqr[i] / t.mip
		#dif.mea <- colMea[i] - t.mea
		#x1[,i] <- (x[,i] / div.iqr) - dif.mea		# mean or median scaling
		x1[,i] <- x[,i] / div.iqr
	}
	colIqr <- apply(x1, 2, iqrange.fn, iqr=iqr, na.rm=T, type=7)
	colMea <- apply(x1, 2, mean, na.rm=T)
	cat("\ncolIqr",round(colIqr,2),"\ncolMea",round(colMea,2),"\n")
	if(pdf.n!=""){
		multi.box.densi.panel(x1, type="bd", iqr=iqr, flags=c(0,0,1,1,1,1), bin.num=50, prefix=paste("NORM IQR",iqr[1],"-",iqr[2]))
	}
	if(out!=""){ assign(out, x1, envir=.GlobalEnv) }
	# create result table
	stat <- data.frame(matrix(0,nr,nc))
	dimnames(stat)[[1]] <- dimnames(x)[[1]]
	dimnames(stat)[[2]] <- dimnames(x)[[2]]
	# bandwidth is starting by +-1 around the mean of the gene row and will be increased by 1 in each direction
	#  until every cell in the row is included - so a 'times increasing count' until every row member is included
	for(j in 1:nr){	# gene rows
		#cat("\nj",j)
		a <- T
		k <- 1	# bandwidth and counter
		l <- 0	# stop indicator
		logi.old <- rep(F,nc)
		x.r <- unlist(x1[j,])
		x.m <- mean(x.r)
		#cat("\nmean",x.m,"\n")
		while(a){
			logi <- (x.r<=(x.m+k) & x.r>=(x.m-k))	# check
			#cat("\nlogi",logi)
			l.t <- sum(logi)
			if(l<l.t){ # one more?
				l.p <- xor(logi.old, logi)   		#xor(c(F,F,F),c(F,T,F))
				stat[j, l.p] <- k
				l <- l.t
				logi.old <- logi
				if(l==nc){ a <- F }
			}
			k <- k+1
		}
		#print(stat)	# development table
	}
	if(pdf.n!=""){
		par(mfrow=c(1,1))
		heatmap.ek.s(x=x1, pc=F, pr=F, values=F, sign.num=2, ctext=NULL, cex.lab=0.5, move.cv=0, move.ch=0,
				rtext=apply(round(t(apply(x1,1,range)),1),1,paste,collapse=":"), move.rh=0, move.rv=0,
				title.m="bandwidth expression per row - after norm IQR", cex=0.7, cex.mat=0.7, mat.srt=0, grid="gray100", diagC=F,
				parmar=c(0.1,8,2,2), color.g=c("white","black"), col.miss.values="grey", center.col=NULL, special.col=NULL,
				k.levels=50, custom.scale=NULL)		# pr color scaling per row
		heatmap.ek.s(x=stat, pc=F, pr=F, values=F, sign.num=2, ctext=NULL, cex.lab=0.5, move.cv=0, move.ch=0,
				rtext=apply(round(t(apply(stat,1,range)),1),1,paste,collapse=":"), move.rh=0, move.rv=0,
				title.m="bandwidth cycle per row - number range", cex=0.7, cex.mat=0.7, mat.srt=0, grid="gray100", diagC=F,
				parmar=c(0.1,8,2,2), color.g=c("white","black"), col.miss.values="grey", center.col=NULL, special.col=NULL,
				k.levels=50, custom.scale=NULL)		# pr color scaling per row
		dev.off()
	}
	return(stat)
}

#a <- test.iqr.m.norm(x=d.exp.test01, row.zero=5, out="aa", iqr=c(0.25,0.75), pdf.n="test_iqr_m_norm__control_data_1.pdf", width=7, height=11)




## function for  per factor optimization
#    expand / shrink until most are inside a defined bandwidth

optimize.per.sample <- function(x, row.zero=5, binn=50){
	# aka normalize the expression data
	# 1) bandwidth raster 2) groups
	# x: matrix data.frame samples in columns
	# binn: >0: fixed bin number =0: static metric 1:-10 10:-100 ... 
	source("../0functions/0general/bin.R")
	nc <- ncol(x)
	# remove all rows <= row.zero counts
	if(!is.na(row.zero)){ x1 <- rm.zero.row(x, row.zero=5) }
	nr <- nrow(x1)
	# create bandwidth	
	if(binn>0){
		kat <- seq(1,2*50,2)
		b.r <- apply(x1,2, bin, bin.num=binn)
	}else{
		kat <- c(1,10,100)
		b.r <- apply(x1,2, bin, bin.fix=c(0,15,400,1000000), fix.open=F)
	}
	print(b.r[[1]]$erg)
	idx <- which((x1[2,1]>=b.r[[1]]$erg[,2] & x1[2,1]<b.r[[1]]$erg[,3]),arr.ind=T)
	kat[idx]
	return(idx)
	
	## todo ##
	
}
#a <- optimize.per.sample(d.exp.test01, row.zero=5)




# mean of row w/o zero elements
# calc factor for becoming mean value w/o zero elements
# analyse factors
test.row.factors.for.mean <- function(x){
	nr <- nrow(x)
	nc <- ncol(x)
	x.m <- vector("numeric",nr)
	for(i in 1:nr){
		x.0 <- x[i,]==0
		r.idx <- which(x.0!=T, arr.ind=T)
		#cat("\n r.idx",r.idx)
		x.zero <- sum(x.0)
		x.m[i] <- sum(x[i,]) / (nc-x.zero)
		y <- rep(0,nc)
		for(j in r.idx){
			y[j] <- x[i,j]/x.m[i]
		}
		x[i,] <- y
	}
	return(list(factor.mat=x, x.mean=x.m))
}
#a <- d.exp.test01[1:3,]
#a[1,6] <- 0
#a[2,16] <- 0
#a <- test.row.factors.for.mean(a)


#d.exp.test01[1:10,]
#aa <- rm.zero.row(d.exp.test01, row.zero=5)
#a <- test.row.factors.for.mean(aa)
#a$x.mean
#a$factor.mat[1:10,]
#hist(a$factor.mat[,1], 16)
#hist(unlist(a$factor.mat[1,]), 16)




# next step: start around x.mean with diff.row /10 +step and return vector with peaks  or direct peaks

# find multi modal
# install.packages("multimode")		also installing the dependencies ‘FNN’, ‘kernlab’, ‘mclust’, ‘multicool’, ‘pracma’, ‘diptest’, ‘ks’, ‘rootSolve’
vector.modes <- function(x, prob=0.95, steps=100){
	require(multimode)
	xlen <- length(x)
	x <- unlist(x)
	if(sum(x==0)>0){ x <- x+1; cat("\nnote: x+1") }	# bw needs to be positive
	x.r <- range(x)
	bw <- seq(x.r[1], x.r[2], length.out=steps)
	s.fn <- sapply(bw, function(i) multimode::nmodes(x, i) )
	par(mfrow=c(2,1)); plot(s.fn, log='y', type='s', yaxt='n'); axis(2, 1:1e5)
	n_mod <- quantile(s.fn, prob)	# estimate - not always proper
	mm <- multimode::locmodes(x, mod0=n_mod)
	hist(x, xlen); abline(v=mm$locations[seq(1,length(a.mm$locations),2)], col="blue"); par(mfrow=c(1,1))
	return(mm$locations[seq(1,length(a.mm$locations),2)])
}
#set.seed(123)
#a <- c(rnorm(250, 0.125, 0.03), rnorm(25, 0.85, 0.05), rnorm(200, 1.24, 0.02))
#vector.modes(x=d.exp.test01[1,], prob=0.96, steps=100)



# install.packages("pracma")
#pracma::Mode(unlist(d.exp.test01[1,]))




## relations between sample values and between the same gene/row values
# in a vector all versus all
all.ver.all <- function(x){
	# relations (quotient)
	xlen <- length(x)
	vlen <- (xlen^2-xlen)/2
	vec <- vector("numeric",vlen)
	k <- 1
	for(i in 1:(xlen-1)){#start, horizontal
		for(j in (i+1):xlen){#end
			vec[k] <- x[i]/x[j]
			k <- k+1
		}
	}
	return(vec)
}
#all.ver.all(x=c(1,2,3,4,5))

all.col.all.row.relations <- function(x){
	# x: data.frame or matrix
	if(is.data.frame(x)){ x <- as.matrix(x) }
	nr <- nrow(x)
	nc <- ncol(x)
	rlen <- (nc^2-nc)/2
	m.rrel <- matrix(0,nr,rlen)
	for(i in 1:nr){
		m.rrel[i,] <- all.ver.all(x[i,])#horizontal
	}
	clen <- (nr^2-nr)/2
	m.crel <- matrix(0,clen,nc)
	for(i in 1:nc){
		m.crel[,i] <- all.ver.all(x[,i])#vertical
	}
	m.rrel[is.infinite(m.rrel)] <- 0
	m.rrel[is.nan(m.rrel)] <- 0
	m.crel[is.infinite(m.crel)] <- 0
	m.crel[is.nan(m.crel)] <- 0
	return(list(rrel=m.rrel,crel=m.crel))
}
#a <- all.col.all.row.relations(x=matrix(1:25,5,5,byrow=T))

#a <- aa[1:4,1:4]
#a <- all.col.all.row.relations(x=a)
#aa[1:4,1:4]
#round(a$rrel,1)
#round(a$crel,1)

#a <- all.col.all.row.relations(x=aa)
#pdf("test_norm_factor_in_x_y_01.pdf", width=7, height=11)
#heatmap.ek.s(as.matrix(a$rrel), pc=F, pr=F, values=F, sign.num=2, ctext=NULL, cex.lab=1, move.cv=0, move.ch=0, rtext=NULL, move.rh=0, move.rv=0,
#		title.m=NULL, cex=1, cex.mat=0.7, mat.srt=0, grid="black", diagC=F,
#		parmar=c(0.01,1,5,0.01), color.g=c("green","red"), col.miss.values="grey", center.col=NULL, special.col=NULL,
#		k.levels=20, custom.scale=NULL)
#
#heatmap.ek.s(as.matrix(a$crel[1:108,]), pc=F, pr=F, values=F, sign.num=2, ctext=NULL, cex.lab=1, move.cv=0, move.ch=0, rtext=NULL, move.rh=0, move.rv=0,
#		title.m=NULL, cex=1, cex.mat=0.7, mat.srt=0, grid="black", diagC=F,
#		parmar=c(0.01,1,5,0.01), color.g=c("green","red"), col.miss.values="grey", center.col=NULL, special.col=NULL,
#		k.levels=20, custom.scale=NULL)
#dev.off()




## simple 'matrix' distances between sample values and between the same gene/row values
# vector-1
naive.m.dist <- function(x){
	xlen <- length(x)
	vlen <- xlen-1
	vec <- vector("numeric",vlen)
	k <- 1
	for(i in 1:(xlen-1)){#start, horizontal
		vec[k] <- abs(x[i]-x[(i+1)])
		k <- k+1
	}
	return(vec)
}
#naive.m.dist(x=c(1,2,3,4,5))

all.col.all.row.nadi <- function(x){
	# x: data.frame or matrix
	if(is.data.frame(x)){ x <- as.matrix(x) }
	nr <- nrow(x)
	nc <- ncol(x)
	clen <- nc-1
	mr.nd <- matrix(0,nr,clen)
	for(i in 1:nr){
		mr.nd[i,] <- naive.m.dist(x[i,])#horizontal
	}
	rlen <- nr-1
	mc.nd <- matrix(0,rlen,nc)
	for(i in 1:nc){
		mc.nd[,i] <- naive.m.dist(x[,i])#vertical
	}
	mr.nd[is.infinite(mr.nd)] <- 0
	mr.nd[is.nan(mr.nd)] <- 0
	mc.nd[is.infinite(mc.nd)] <- 0
	mc.nd[is.nan(mc.nd)] <- 0
	return(list(rnd=mr.nd,cnd=mc.nd))
}
#a <- all.col.all.row.nadi(x=matrix(1:25,5,5,byrow=T))

#a <- aa[1:4,1:4]
#a <- all.col.all.row.nadi(x=a)
#aa[1:4,1:4]
#round(a$rnd,1)
#round(a$cnd,1)

#a <- all.col.all.row.nadi(x=aa)
## x lattice distance
#round(range(a$rnd,1),1)		# 0 1565
#hist(a$rnd)
#a1 <- bin(a$rnd, bin.num=50)
#a2 <- find.peak(a1$erg[,1], lowlim=3)
#plot(1:a1$para[1],a1$erg[,1]); abline(v=a2$peak.pos, col="blue")
## y lattice distance
#round(range(a$cnd,1),1)		# 0 1623
#hist(a$cnd)
#a1 <- bin(a$cnd, bin.num=50)
#a2 <- find.peak(a1$erg[,1], lowlim=3)
#plot(1:a1$para[1],a1$erg[,1]); abline(v=a2$peak.pos, col="blue")



## peak in row
#a1 <- bin(aa[1,], bin.num=16)
#a2 <- find.peak(a1$erg[,1], lowlim=0)
#plot(1:a1$para[1],a1$erg[,1]); abline(v=a2$peak.pos, col="blue")

all.peaks.all.rows <- function(x){
	# x: data.frame or matrix
	f1 <- function(x, lenrow){
		vec <- rep(0,lenrow)
		vec[x] <- 1
		return(vec)
	}
	if(is.data.frame(x)){ x <- as.matrix(x) }
	nr <- nrow(x)
	nc <- ncol(x)
	pnum <- vector("integer",nr)
	groups <- matrix(0,nr,nc)
	peaks <- matrix(0,nr,nc)
	for(i in 1:nr){
		ll <- find.peak(x[i,], lowlim=0)
		pnum[i] <- ll$gnumber
		groups[i,] <- ll$groups
		peaks[i,] <- f1(ll$peak.pos, nc)
	}
	return(list(pnum=pnum, groups=groups, peaks=peaks))
}

#a <- all.peaks.all.rows(x=aa)




## interpolation fitting polynomial model again
# https://stackoverflow.com/questions/3822535/fitting-polynomial-model-to-data-in-r
# Which model is the "best fitting model" depends on what you mean by "best".
# R has tools to help, but you need to provide the definition for "best" to choose between them.
		
#ax <- 1:10
#ay <- ax + c(-0.5,0.5)
#
#plot(ax,ay, xlim=c(0,11), ylim=c(-1,12))
#
#fit1 <- lm( ay~offset(ax) -1 )
#fit2 <- lm( ay~ax )
#fit3 <- lm( ay~poly(ax,3) )
#fit4 <- lm( ay~poly(ax,9) )
#library(splines)
#fit5 <- lm( ay~ns(ax, 3) )
#fit6 <- lm( ay~ns(ax, 9) )
#
#fit7 <- lm( ay ~ ax + cos(ax*pi) )
#
##axx <- seq(0,11, length.out=250)
#lines(ax, predict(fit1, data.frame(x=ax)), col='blue')
#lines(ax, predict(fit2, data.frame(x=ax)), col='green')
#lines(ax, predict(fit3, data.frame(x=ax)), col='red')
#lines(ax, predict(fit4, data.frame(x=ax)), col='purple')
#lines(ax, predict(fit5, data.frame(x=ax)), col='orange')
#lines(ax, predict(fit6, data.frame(x=ax)), col='grey')
#lines(ax, predict(fit7, data.frame(x=ax)), col='black')





## basic functions again    <->  functions02.R
mean.n.dev <- function(x){
	# mean norm deviation
	n <- length(x)
	xm <- mean(x, na.rm=T)
	d <- abs( (x-xm) ) / n	# absolut deviation
	nd <- d/xm	# norm by mean
	return(nd)
}

m.aad <- function(x){
	# aad mean norm
	n <- length(x)
	xm <- mean(x, na.rm=T)
	aad <- sum( abs( (x-xm) ), na.rm=T ) / n	# average absolute deviation
	naad <- aad/xm	# norm by mean => (coefficient of) dispersion
	return(naad)
}
# with SD/mean -> coefficient of variation
#m.aad(c(1.1,0.95,0.99,1.2, 1.05, 1,1.1,0.88,0.95,1.1,0.98))	# 0.073  small fluctuation
#m.aad(c(1.5,0.85,0.9,1.4,1.15,1,1.3,0.78,0.9,1.2,0.88))		# 0.195  large fluctuation
#m.aad(c(1.1,0.95,0.99,1.2, 1.8,  1,1.1,0.88,0.95,1.1,0.98))	# 0.136  small fluctuation with one large peak


max.n.dev <- function(x, basicX=T){
	# max norm deviation
	n <- length(x)
	xm <- mean(x, na.rm=T)
	d <- abs( (x-xm) ) / n	# absolut deviation
	if(basicX){ xmax <- max(x) }else{ xmax <- max(d) }
	nd <- d/xmax	# norm by max
	return(nd)
}

max.aad <- function(x){
	# aad max norm
	n <- length(x)
	xm <- mean(x, na.rm=T)
	aad <- sum( abs( (x-xm) ), na.rm=T ) / n	# average absolute deviation
	xmax <- max(x)
	naad <- aad/xmax	# norm by max
	return(naad)
}


n.sum <- function(x){
	# sum by n norm
	n <- length(x)
	nsum <- sum(x)/n	# absolute sum := mean
	return(nsum)
}




## peak detection
# see 0functions/prox01/functions00.R

# ?? normalize by  x/(max-min)  or  x/max



## grid like I
# intra vector : neighbouring positions 1-2, 2-3, 3-4, ...
neighbour.dist <- function(x){
	# distance between the neighbouring vector positions (so not all possible ones)
	rx <- range(x)
	if(rx[1]<0 | rx[2]<0){ stop("x has negative numbers") }
	if(rx[2]==rx[1]){ stop("range is zero") }
	xlen <- length(x)
	dlen <- xlen-1
	vec <- vector("numeric",dlen)
	k <- 1
	for(i in 1:dlen){
		vec[k] <- x[i]-x[(i+1)]
		k <- k+1
	}
	return(vec)
}
#neighbour.dist(x=c(1,2,3,4,5))		# minus is up, plus is down

sample.neighbour.dist <- function(x, samp=0, seedn=NULL){
	# x: dist vector,  samp: number of samplings per column
	# seedn: if not NULL needs to be an integer e.g. 174
	if(samp==0){ stop("no sampling") }
	xlen <- length(x)
	if(!is.null(seedn)){ set.seed(seedn) }
	mats <- matrix(0, xlen, samp)
	for(i in 1:samp){
		mats[,i] <- sample(x)		# random permutation
	}
	return(mats)
}
#sample.neighbour.dist(x=c(1,2,3,4,5), samp=2, seedn=174)



## grid like II
# search for extremes and normals
row.with.wo.extremes <- function(x, minf=2, maxf=1.2, naad=0.4, su=10){
	# search for rows with and without extremes
	# x: matrix with expression values
	# minf: multiplication smaller,   maxf: multiplication larger
	# naad: mean normalized average absolut deviation value
	# su: sum of counts
	source("../0functions/prox01/test_norm.R")
	source("../0functions/sets/intersect_e.R")
	nr <- nrow(x)
	mat <- data.frame(matrix(0,nr,14))
	dimnames(mat)[[1]] <- dimnames(x)[[1]]
	dimnames(mat)[[2]] <- c("row","pos_min","pos_max","median","n.aad","min","1rest_min","1rest_max","1diff_rest","2diff_rest","2rest_min","2rest_max","max","sum")
	for(i in 1:nr){
		mat[i,1] <- i
		xx <- unlist(x[i,])
		mat[i,2] <- which(min(xx)==xx, arr.ind=T)[1]	# might be more equal ones, take the first
		mat[i,3] <- which(max(xx)==xx, arr.ind=T)[1]
		mat[i,4] <- median(xx)
		mat[i,5] <- m.aad(xx)
		mat[i,6] <- min(xx)
		rr <- range(xx[(min(xx)!=xx)])
		mat[i,7] <- rr[1]
		mat[i,8] <- rr[2]
		mat[i,9] <- rr[2]-rr[1]
		rr <- range(xx[(max(xx)!=xx)])
		mat[i,10] <- rr[2]-rr[1]
		mat[i,11] <- rr[1]
		mat[i,12] <- rr[2]
		mat[i,13] <- max(xx)
		mat[i,14] <- sum(xx)
	}
	# filter - max peak
	maxp <- mat[(mat[,"max"]>(mat[,"2rest_max"]*maxf))&(mat[,"sum"]>su),]
	# filter - min peak
	minp <- mat[(mat[,"min"]<(mat[,"1rest_min"]/minf))&(mat[,"sum"]>su),]
	# intersect
	i.e <- intersect.e(x=maxp, y=minp)
	cat("\n maxp-minp intersection size:", nrow(i.e$xy),"(>0 to mmp) / maxp",nrow(maxp),"/ minp",nrow(minp),"\n")
	maxp <- i.e$xs
	minp <- i.e$ys
	mmp <- i.e$xy
	# order
	maxp <- maxp[order(maxp[,"n.aad"],decreasing=T),]
	minp <- minp[order(minp[,"n.aad"],decreasing=T),]
	# filter - very uniform
	loo <- mat[(mat[,"n.aad"]<naad),]
	# intersect
	i.e <- intersect.e(x=loo, y=maxp)
	cat("\n loo-maxp intersection size:", nrow(i.e$xy),"(>0 to mmp) / loo",nrow(loo),"\n")
	loo <- i.e$xs
	mmp <- rbind(mmp,i.e$xy)
	i.e <- intersect.e(x=loo, y=minp)
	cat("\n loo-minp intersection size:", nrow(i.e$xy),"(>0 to mmp)\n")
	loo <- i.e$xs
	mmp <- rbind(mmp,i.e$xy)
	# order
	loo <- loo[order(loo[,"n.aad"],decreasing=F),]
	
	return(list(maxp=maxp,minp=minp,mmp=mmp,loo=loo))
}
#aa[1:5,];row.with.wo.extremes(aa[1:5,], minf=2, maxf=1.2, naad=0.6, su=10)


# distance
all.neighbour.dist <- function(x){
	# all versus all (no redundancy) -- one sample vector
	rx <- range(x)
	if(rx[1]<0 | rx[2]<0){ stop("x has negative numbers") }
	if(rx[2]==rx[1]){ warning("range is zero") }
	xlen <- length(x)
	num.res <- all.neighbour.dist.sub1(xlen)
	dlen <- xlen-1
	vec <- vector("numeric", num.res)
	k <- 1
	for(i in 1:dlen){
		for(j in (i+1):xlen){
			vec[k] <- x[i]-x[j]
			k <- k+1
		}
	}
	return(vec)
}
#all.neighbour.dist(x=c(1,2,3,4,5))		# no diag, two coordinates, both increasing, plus-minus numbers
#all.neighbour.dist(x=c(0,0,0,0,0))

all.neighbour.dist.sub1 <- function(xlen){
	# result length
	num.res <- ((xlen^2) - xlen) / 2
	return(num.res)
}
all.neighbour.dist.sub2 <- function(xlen){
	# combination label
	num.res <- all.neighbour.dist.sub1(xlen)
	dlen <- xlen-1
	mat <- matrix(0, num.res, 2)
	k <- 1
	for(i in 1:dlen){
		for(j in (i+1):xlen){
			mat[k,] <- c(i,j)
			k <- k+1
		}
	}
	return(mat)
}


#
lines.neighbour.dist <- function(x){
	# all versus all others (redundancy, but all lines are separately) -- matrix of lines
	rx <- range(x)
	if(rx[1]<0 | rx[2]<0){ stop("x has negative numbers") }
	if(rx[2]==rx[1]){ stop("range is zero") }
	xlen <- length(x)
	num.res <- xlen-1
	mat <- matrix(0, num.res, xlen)
	for(i in 1:xlen){
		x1 <- x[i]
		x2 <- x[-i]
		for(j in 1:num.res){
			mat[j,i] <- x1-x2[j]
		}
	}
	return(mat)
}
#lines.neighbour.dist(x=c(1,2,3,4,5))		# no diag, two coordinates, both increasing, plus-minus numbers

line.neighbour.dist <- function(x, sel){
	# only selected versus rest -- one sample vector
	# sel: position number - part of x
	rx <- range(x)
	if(rx[1]<0 | rx[2]<0){ stop("x has negative numbers") }
	if(rx[2]==rx[1]){ stop("range is zero") }
	xlen <- length(x)
	if(sel>xlen|sel<=0){ stop("sel out of range") }
	num.res <- xlen-1
	vec <- vector("numeric", num.res)
	x1 <- x[sel]
	x2 <- x[-sel]
	for(i in 1:num.res){
		vec[i] <- x1-x2[i]
	}
	return(vec)
}
#line.neighbour.dist(x=c(1,2,3,4,5),sel=2)		# sel versus others, order outside sel preserved, plus-minus numbers

line.neighbour.dist.sub1 <- function(xlen){
	# result length
	num.res <- xlen-1
	return(num.res)
}
line.neighbour.dist.sub2 <- function(xlen, sel){
	# combination label
	num.res <- line.neighbour.dist.sub1(xlen)
	mat <- matrix(0, num.res, 2)
	pos <- 1:xlen
	pos.m.sel <- pos[-sel]
	for(i in 1:num.res){
		mat[i, ] <- c(sel,pos.m.sel[i])
	}
	return(mat)
}
#line.neighbour.dist.sub2(xlen=9, sel=2)


grid.distribution <- function(x, y){
	# difference
	xlen <- length(x)
	ylen <- length(y)
	if(xlen!=ylen){ stop("x and y need to have the same length") }
	z <- x-y	# might be positive or negative
	plot(1:xlen,z,type="h",xlab="",ylab="",col="blue")
	return()
}
# grid.distribution(x=c(1,3,8,9,2), y=c(2,3,7,11,4))





## algo -1-  use distances to normalize data

one.group.norm <- function(x, ref=1, ssize=200, perc=c(25,75), show=F, sel.factor=F){
	# assumed one group for norm of expression data
	# x: matrix of expression values, samples in columns
	# ref: reference column number  , ssize: sampling size is starting at > ssize
	# perc: difference band d as percent range
	# show: T: all, F: minimum
	if(is.data.frame(x)){ x <- as.matrix(x) }
	nr <- nrow(x)
	nc <- ncol(x)
	if(nc<2){ stop("2 and more columns") }
	if(ssize<3){ stop("3 and more features") }
	# create distance grid
	if(nr>ssize){	#index selection
		idx.sel <- sample(nr,ssize)
		#idx.sel <- 1:200
		idxS <- T
	}else{
		idx.sel <- 1:nr
		idxS <- F
	}
	if(show){ cat("\nidx.sel",idx.sel,"\n") }
	xlen <- length(idx.sel)
	rn <- dimnames(x)[[1]]
	cn <- dimnames(x)[[2]]
	num.res <- all.neighbour.dist.sub1(xlen)
	idx.p <- all.neighbour.dist.sub2(xlen)	# matrix, all possible index pairs of the ssize
	di <- matrix(0,num.res,nc)
	for(i in 1:nc){	# distance
		di[,i] <- all.neighbour.dist(x[idx.sel,i])	# signed!
	}
	# all sample pairings
	ncc <- 1:nc
	nsel <- ncc[!(ref==ncc)]
	m.sub <- vector("list",(nc-1))
	k <- 1
	for(i in nsel){
		# difference pairs between distances, based on ref, filtered on range (and NaN,Inf ...)
		m.sub[[k]] <- one.group.norm.sub1(di[,ref], di[,i], idx.p, perc, num.res, show)
		k <- k+1
	}
	# mark samples where no indices remain
	m.sub.note <- one.group.norm.sub3(m.sub, show)
	#m.sub.note <- c(2,3)	# test line
	# ranking bandwidth overlap
	idx.ov <- one.group.norm.sub4(m.sub, m.sub.note)	# matrix, rows: distance element
	# return selected factors for checks
	if(sel.factor){ return(one.group.norm.sub5b(idx.ov, idx.sel, idxS)) }
	# create scaling factors
	m.factors <- one.group.norm.sub5(x, idx.ov, idx.sel, idxS, m.sub.note, rn, show)
	# adjust all or subset of columns (missing means peculiarities or problems)
	x <- one.group.norm.sub6(x, m.factors, m.sub.note, ncc)
	return(x)
}
#a <- one.group.norm(x=aa[1:4,1:3], ref=1, ssize=3, perc=c(25,75), show=T)

one.group.norm.sub1 <- function(x, y, idx.p, perc, num.res, show){
	# create weight
	mz <- one.group.norm.sub2(x, y, fn="max")		# "mean"
	# difference
	d <- x - y
	dz <- mz / d
	dzr <- rank(dz)
	dzr.ra <- range(dzr)
	dzr.ra.thres.lo <- (dzr.ra[2] * (perc[1])) / 100
	dzr.ra.thres.up <- (dzr.ra[2] * (perc[2])) / 100
	if(show){ cat("\ndzr.ra",dzr.ra," unique level ",length(unique(dzr))," of ",num.res,"positions, thres",dzr.ra.thres.lo," ",dzr.ra.thres.up) }
	# filter on certain percent range
	dzr.sel <- dzr>=dzr.ra.thres.lo & dzr<=dzr.ra.thres.up
	sub <- data.frame(idxp1=idx.p[dzr.sel,1], idxp2=idx.p[dzr.sel,2], dzS=dz[dzr.sel], rS=dzr[dzr.sel])	# might not be a data frame ->> nrow is NULL not only 0,1,..,n
	# filter additionally on rows which are not NaN or Inf in dz (column dzS)
	sub <- sub[is.finite(sub[,"dzS"]), ]
	# there could be the case of no distance remaining (discarded in .sub3)  or  distances = 0
	return(sub)
}

one.group.norm.sub2 <- function(x, y, fn){
	# weight
	x1 <- cbind(x,y)
	fno <- get(fn, envir=.GlobalEnv)
	z <- apply(x1,1,fno)
	return(z)
}

one.group.norm.sub3 <- function(m.sub, show){
	# mark non informative samples
	#  (not now: and low row number samples)
	xlen <- length(m.sub)	# list
	ov <- vector("integer",xlen)
	note <- NULL
	for(i in 1:xlen){
		out <- nrow(m.sub[[i]])
		if(out==0){ note <- c(note,i) }
		ov[i] <- out
	}
	# plot
	if(show){
		cat("\nsmallest # of rows",min(ov[ov>0]),"\n")
		str.x <- paste("pool - number of remaining rows per sample (ordered)",sep="")
		str.y <- paste("min-max",min(ov),max(ov),sep=" ")
		str.t <- paste("# of samples with 0 rows :",length(note),"/",length(ov),"\nsample# ",paste(note,collapse=","),sep=" ")
		plot(1:length(ov),ov[order(ov)], type="h", xlab=str.x, ylab=str.y, xlim=c(1,length(ov)), ylim=c(0,max(ov)), col="blue", las=2)
		title(main=str.t)
	}
	return(note)
}

one.group.norm.sub4 <- function(m.sub, note){
	# overlap
	source("../0functions/sets/iu.set.R")
	source("../0functions/0general/list2df.R")
	if(length(note)>=1){ x <- m.sub[-note] }else{ x <- m.sub }	# m.sub NULL -> length 0
	xlen <- length(x)
	jj <- NULL
	for(i in 1:xlen){
		jj[[i]] <- paste(x[[i]][,1],x[[i]][,2],sep="-")
	}
	uniD <- iu.set(x=jj, idv=rep(0,xlen), setfn="ic")
	ssp <- strsplit(uniD, split="-", fixed=T)
	uniM <- list2mat( ssp, rows=T )
	mode(uniM) <- "integer"		# character to integer
	return(uniM)
}

one.group.norm.sub5 <- function(x, idx.ov, idx.sel, idxS, note, rn, show){
	# adjustment factors
	#  mean over band selected features
	# adjust x with note
	if(length(note)>=1){ x1 <- x[ , -note] }else{ x1 <- x }
	nc <- ncol(x1)
	#
	io.len <- nrow(idx.ov)
	if(idxS){	# translate
		for(i in 1:io.len){
			idx.ov[i,1] <- idx.sel[idx.ov[i,1]]
			idx.ov[i,2] <- idx.sel[idx.ov[i,2]]
		}
	}
	if(show){ cat("after translate, rows 1:15\n"); print(idx.ov[1:15,]) }
	# expr. tag mean
	uni.idx.ov <- unique(as.integer(idx.ov))	# remove idx duplicates
	uni.idx.ov <- uni.idx.ov[order(uni.idx.ov)]
	if(show){ cat("\nuni.idx.ov ",uni.idx.ov) }
	uni.idx.ov.l <- length(uni.idx.ov)
	m.expr.tag <- vector("integer",uni.idx.ov.l)
	for(i in 1:uni.idx.ov.l){
		m.expr.tag[i] <- mean(x1[uni.idx.ov[i],])
	}
	names(m.expr.tag) <- rn[uni.idx.ov]
	if(show){ cat("\nm.expr.tag ",m.expr.tag,"\nnames: ",rn[uni.idx.ov]) }
	# mean of tag deviation factors
	features.factors <- matrix(0,uni.idx.ov.l,nc)
	for(i in 1:uni.idx.ov.l){
		features.factors[i,] <- x1[uni.idx.ov[i], ] / m.expr.tag[i]	# per sample c(1,2,4,0,1)/2
	}
	# mean of factors per sample
	m.factors <- colMeans(features.factors)
	if(show){ cat("\nm.factors ",m.factors) }
	return(m.factors)
}

one.group.norm.sub5b <- function(idx.ov, idx.sel, idxS){
	# return the real index numbers for export
	io.len <- nrow(idx.ov)
	if(idxS){	# translate
		for(i in 1:io.len){
			idx.ov[i,1] <- idx.sel[idx.ov[i,1]]
			idx.ov[i,2] <- idx.sel[idx.ov[i,2]]
		}
	}
	return(idx.ov)
}

one.group.norm.sub6 <- function(x, m.factors, note, ncc){
	# adjust columns with functional m.factor and leave the others unchanged
	if(length(note)>=1){ nca <- ncc[-note] }else{ nca <- ncc }
	unch <- NULL
	chan <- NULL
	k <- 1
	for(i in nca){
		if(m.factors[k]>0){		# m.factor=0 -> unchanged
			x[,i] <- x[,i] / m.factors[k]
			chan <- c(chan,i)
		}else{
			unch <- c(unch,i)
		}
		k <- k+1
	}
	# feedback
	lnote <- length(note)
	lunch <- length(unch)
	lchan <- length(chan)
	lall <- ncol(x)
	cat("\nadjusted #",lchan," - ",chan,"\nnot adj. diff. calc. #",lnote," - ",note,"\nnot adj. m.factor 0  #",lunch," - ",unch,"\nof  ",lall,"\n")
	# color vector for levelplot()
	gord <- c(chan,note,unch)
	gcol <- c(rep("green",lchan),rep("black",lnote),rep("black",lunch))
	gcol <- gcol[order(gord)]
	assign("acol",gcol,envir=.GlobalEnv)
	return(x)
}

#aa[1:4,1:3]
#p1water ps_U343 v1water
#piR-hsa-119247      117      72      68
#piR-hsa-2781178      78      20       1
#piR-hsa-220202       70      19       2
#piR-hsa-90490        66      18       2
#
#idx.sel 2 1 3
#
#d 13 7 -6
#dz -3 1.142857 -8.833333
#dzr 2 3 1
#dzr.ra 1 3  unique level  3  of  3 positions, thres 0.75   2.25 
#dzr.sel TRUE FALSE TRUE
#
#d 28 9 -19
#dz -1.392857 0.8888889 -3.473684
#dzr 2 3 1
#dzr.ra 1 3  unique level  3  of  3 positions, thres 0.75   2.25 
#dzr.sel TRUE FALSE TRUE
#
#Names : 
#		1 2 
#
#idx before
#[,1] [,2]
#1    1    2
#2    2    3
#after
#[,1] [,2]
#1    2    1
#2    1    3
#
#uni.idx.ov  1 2 3
#m.expr.tag  85.66667 33 30.33333  names:  piR-hsa-119247 piR-hsa-2781178 piR-hsa-220202
#m.factors  2.012362 0.6909671 0.2966705
#
#p1water   ps_U343   v1water
#piR-hsa-119247     58.14    104.20    229.21
#piR-hsa-2781178    38.76     28.94      3.37
#piR-hsa-220202     34.78     27.49      6.74
#piR-hsa-90490      32.79     26.05      6.74



## algo -2-  different approach  with  optimization

b.one.group.norm <- function(x, ssize=200, pu=0.95, pd=0.05, p.dist=c(0.25,0.75), lib.norm=NULL, show=3, sel.factor=F){
	# assumed one group for norm of expression data
	# x: matrix of expression values, samples in columns
	# ssize: sampling size is starting at > ssize
	# pu,pd: upper and lower trim p factor
	# lib.norm: normalize library size: NULL: no normalization, or 0: max count column taken, or number: column number
	# show: 1..n
	# sel.factor: return values
	x1 <- b.one.group.norm.sub0(x, ssize)	# check x
	# create distance grid
	if(x1$nr>ssize){	#index selection
		#idx.sel <- sample(x1$nr,ssize)
		idx.sel <- 1:200
		idxS <- T
	}else{
		idx.sel <- 1:x1$nr
		idxS <- F
	}
	if(show==1){ cat("\nidx.sel",idx.sel,"\n") }
	# level total counts per sample (no -> !concern)
	# filter on expression range
	x2 <- b.one.group.norm.sub2(x2, pu, pd)
	if(is.null(x2)){ stop("filter on expression range : nothing left") }
	return(x2)
	
}
#b.one.group.norm(x, ssize=200, p.expr=c(0,0.05,0.95,1), p.dist=c(0.25,0.75), lib.norm=F, show=3, sel.factor=F)


b.one.group.norm.sub0 <- function(x, ssize){
	# check basics of x
	if(is.data.frame(x)){ x <- as.matrix(x) }
	if(anyNA(x)){ stop("NA in x not allowed") }		# alt.: stopifnot all expression T stop at first
	col.sums <- apply(x,2,sum)
	cat("\ncolumn sums",col.sums)
	nr <- nrow(x)
	col.zero <- apply(x,2, function(x){ sum(x==0) })
	cat("\ncolumn zero share",col.zero,"\nof length",nr,"\n")
	if(any(col.sums==0)){ stop("remove columns which have 0 counts") }
	nc <- ncol(x)
	if(nc<2){ stop("x: 2 and more columns") }
	if(ssize<3){ stop("ssize: 3 and more features") }
	return(list(x=x,nc=nc,nr=nr))
}
#b.one.group.norm.sub0(matrix(c(0,0,0,1,1,0,1,1,1),3,3,byrow=F), 3)

b.one.group.norm.sub1 <- function(x, lib.norm=0){				# not better than other
	# adjust all samples to a common count (strength) level
	#  (like DESeq and others with the 'library size' correction)
	#  but here the approach tries to tries to distribute broader
	#  by not overemphasizing with global multiplication factors-
	#  in any case, the basic idea contains several problems
	# x: expression matrix or part of,   lib.norm: column number or 0 (max count column)
	nr <- nrow(x)
	nc <- ncol(x)
	#cat("\nnr",nr," nc",nc," colSums",colSums(x),"\n")
	pos.0 <- apply(x,2,function(s){ s==0 })
	pos.1 <- !pos.0
	#print(pos.0); print(pos.1)
	z <- apply(x,2,sum)		# column sum
	if(lib.norm==0){
		zsel <- max(z)
		zselI <- which(zsel==z, arr.ind=T)		# could be >1 positions
	}else{
		zsel <- z[lib.norm]
		zselI <- lib.norm
	}
	ncc <- 1:nc
	zI <- ncc[!(zselI==ncc)]
	zd <- zsel - z			# column differences
	zmf <- vector("numeric",nc)
	for(i in 1:nc){			# mean of filled positions
		zmf[i] <- z[i] / sum(pos.1[,i])
	}
	sum.p1m <- sum(pos.1[,zselI])	# selected count sample: filled positions (>0)
	#cat("z",z," zsel",zsel," zI",zI," zselI",zselI," zd",zd," zmf",zmf," sum.p1m",sum.p1m)
	for(i in zI){
		tmpI0 <- logic.dir(pos.0[,zselI], pos.0[,i], 1)		# true pos.0 positions cf. reference column
		if(sum(tmpI0)!=0){
			tmpS <- sum(tmpI0)			# number of true pos.0
			z.s <- tmpS / sum.p1m		# share which goes from sel-z.1 to z.0 part, equally distributed
			z.sc <- zd[i] * z.s			# share in counts
			z.pos <- z.sc / tmpS		# share in counts per every true pos.0
			mx <- min(x[x[,i]>0,i])
			if(z.pos>mx){ z.pos <- mx }			# limit to min >0 in sample
			zd[i] <- zd[i] - (z.pos * tmpS)		# correct zd
			x[tmpI0,i] <- z.pos					# assign to all true pos.0
			#cat("\nz.s",z.s," z.sc",z.sc," z.pos",z.pos," mx",mx," i",i," tmpI0",tmpI0)
		}
	}
	#cat("\n"); print(x); cat("colSums",colSums(x))
	for(i in zI){
		tmpI0 <- logic.dir(pos.0[,zselI], pos.0[,i], 1)		# true pos.0 positions cf. reference column
		tmpI1 <- logic.dir(pos.1[,zselI], pos.1[,i], 2)		# true pos.1 positions cf. reference column
		#cat("\nsum(x[tmpI0,i])",sum(x[tmpI0,i])," sum(x[tmpI1,i])",sum(x[tmpI1,i]))
		rs <- z[zselI] - sum(x[tmpI0,i])	# down corrected sum of sel-sample-sum
		zf <- 1 / (sum(x[tmpI1,i]) / rs)	# multiplication factor for all true pos.1
		x[tmpI1,i] <- x[tmpI1,i] * zf		# application on true pos.1
		#cat("\nrs",rs," zf",zf," z[zselI]",z[zselI]," i",i," tmpI0",tmpI0," tmpI1",tmpI1)
	}
	cat("\n#",zselI," column sums",colSums(x),"\n")
	return(x)
}
#matrix(c(1,0,0,0,2,2,0,0,3,3,0,3),4,3,byrow=F); b.one.group.norm.sub1(matrix(c(1,0,0,0,2,2,0,0,3,3,0,3),4,3,byrow=F), 2)	# 0 2
#a <- b.one.group.norm.sub1(aa, 14)

logic.dir <- function(ref,x, r){
	# a directed decision scheme ref -> x for logical operations
	# ref,x: two vectors, ref: reference column, x: target column
	# r: rule: row number of rule
	rlen <- length(ref)
	xlen <- length(x)
	if(rlen!=xlen){ stop("logic.dir() x,y needs to be the same length") }
	rule <- matrix(c(	# ref,x,outcome
					F,T,T,
					T,T,T),
			2,3,byrow=T)
	f1 <- function(x, z){
		if(z[1]==x[1] & z[2]==x[2]){ return(T) }else{ return(F) }
	}
	mat <- cbind(ref,x)
	out <- apply(mat,1, f1, z=rule[r,])
	return(out)
}
#logic.dir(ref=T,x=T, r=1); logic.dir(ref=c(F,T,T,F),x=c(T,F,T,F), r=1); logic.dir(ref=c(F,T,T,F),x=c(T,F,T,F), r=2)

b.one.group.norm.sub1.1 <- function(x, lib.norm=0, pu=0.95, pd=0.05){				# slightly better than sub1, but still not the right direction
	# adjust all samples to a common count (strength) level
	# by simple multiplication factors (spread) with 10% trim before
	#  (like DESeq and others with the 'library size' correction)
	#  in any case, the basic idea contains several problems
	# x: expression matrix or part of,   lib.norm: column number or 0 (max count column)
	# pu,pd: upper and lower trim proportion 0..1
	nc <- ncol(x)
	trim.xi <- function(x, pu, pd){
		pup <- quantile(x, pu)
		pdown <- quantile(x, pd)
		return(c(pup,pdown))
	}
	lim <- apply(x,2, trim.xi, pu=pu, pd=pd)	# pu first row, pd second row
	zt <- vector("numeric",nc)
	for(i in 1:nc){
		trunc <- x[x[i] <= lim[1,i] & x[i] >= lim[2,i], i]
		zt[i] <- sum(trunc)
	}
	if(lib.norm==0){
		zsel <- max(zt)
		zselI <- which(zsel==zt, arr.ind=T)		# could be >1 positions
	}else{
		zsel <- zt[lib.norm]
		zselI <- lib.norm
	}
	zf <- zt[zselI[1]] / zt
	for(i in 1:nc){
		x[,i] <- x[,i] * zf[i]
	}
	cat("\n#",zselI," column sums",colSums(x),"\n")
	return(x)
}
#matrix(c(1,0,0,0,2,2,0,0,3,3,0,3),4,3,byrow=F); b.one.group.norm.sub1.1(matrix(c(1,0,0,0,2,2,0,0,3,3,0,3),4,3,byrow=F), 0, pu=0.95, pd=0.05)	# 0 2
#a <- b.one.group.norm.sub1.1(aa, 0, pu=0.975, pd=0.025)  # 0 14


b.one.group.norm.sub2 <- function(x, pu=0.95, pd=0.05){
	# filter expression range
	# x: expression matrix
	# pu,pd: 2 cut values
	nc <- ncol(x)
	trim.xi <- function(x, pu, pd){
		pup <- quantile(x, pu)
		pdown <- quantile(x, pd)
		return(c(pup,pdown))
	}
	lim <- apply(x,2, trim.xi, pu=pu, pd=pd)	# pu first row, pd second row
	
	
	return(x)
}
#a <- b.one.group.norm.sub1(x=aaa, prob=c(0, 0.02, 0.98, 1))	# aaa[,1,drop=F]




## algo -3-  different approach  with  optimization - but w/o norm + filter in the beginning

c.one.group.norm <- function(x, ssize=200, sampleN=0){
	# assumed one group for norm of expression data
	# x: matrix of expression values, samples in columns
	# ssize: sampling size is starting at > ssize
	# sampleN: number samplings
	x1 <- c.one.group.norm.sub0(x, ssize)	# check x
	# create distance grid
	if(x1$nr>ssize){	#index selection
		#idx.sel <- sample(x1$nr,ssize)
		idx.sel <- 1:200
		idxS <- T
	}else{
		idx.sel <- 1:x1$nr
		idxS <- F
	}
	#cat("\nidx.sel",idx.sel,"\n")
	# level total counts per sample (no -> !concern)
	# filter on expression range (no -> !concern)
	nr <- x1$nr
	nc <- x1$nc
	x <- x1$x
	xlen <- length(idx.sel)
	rn <- dimnames(x)[[1]]
	cn <- dimnames(x)[[2]]
	# dist
	num.res <- all.neighbour.dist.sub1(xlen)	# length of result vector
	idx.p <- all.neighbour.dist.sub2(xlen)	# matrix, all possible index pairs of the ssize
	di <- matrix(0,num.res,nc)
	for(i in 1:nc){	# distance
		di[,i] <- all.neighbour.dist(x[idx.sel,i])	# signed!
	}
	# start values
	start <- c.one.group.norm.sub1.1(di, num.res, nc)
	#return(start)
	erg <- c.one.group.norm.sub2(di, num.res, nc, start$gsadm, sampleN, start$rz)
	return(list(first=start,sample=erg))
}
#a <- c.one.group.norm(aa, ssize=200, sampleN=10000)		# 10000 -- 1 min

# dist sum sample
#dist.sum.sample <- apply(di,2,sum)
# norm sum sample
#di.n <- matrix(0,num.res,nc)
#for(i in 1:nc){
#	di.n[,i] <- di[,i] / dist.sum.sample[i]
#}

c.one.group.norm.sub0 <- function(x, ssize){
	# check basics of x
	if(is.data.frame(x)){ x <- as.matrix(x) }
	if(anyNA(x)){ stop("NA in x not allowed") }		# alt.: stopifnot all expression T stop at first
	col.sums <- apply(x,2,sum)
	cat("\ncolumn sums",col.sums)
	nr <- nrow(x)
	col.zero <- apply(x,2, function(x){ sum(x==0) })
	cat("\ncolumn zero share",col.zero,"\nof length",nr,"\n")
	if(any(col.sums==0)){ stop("remove columns which have 0 counts") }
	nc <- ncol(x)
	if(nc<2){ stop("x: 2 and more columns") }
	if(ssize<3){ stop("ssize: 3 and more features") }
	return(list(x=x,nc=nc,nr=nr))
}
#c.one.group.norm.sub0(matrix(c(0,0,0,1,1,0,1,1,1),3,3,byrow=F), 3)

#c.one.group.norm.sub1 <- function(x, num.res, nc){
#	# start values
#	dv.mat <- matrix(0,num.res,nc)
#	rm <- rowMeans(x)		# row means over di.n
#	dm <- dv.mat
#	for(i in 1:nc){		# difference to row means
#		dm[,i] <- x[,i] - rm[i]
#	}
#	adm <- dv.mat
#	for(i in 1:nc){		# absolute differences
#		adm[,i] <- abs(dm[,i])
#	}
#	sadm <- apply(adm,1, sum)
#	gsadm <- sum(sadm)
#	madm <- rowMeans(adm)	# mean of absolute differences
#	dmadm <- dv.mat
#	for(i in 1:nc){		# adm difference to madm
#		dmadm[,i] <- adm[,i] - madm[i]
#	}
#	#sdmadm <- apply(dmadm,2, function(x, num.res){ sum(x)/num.res }, num.res=num.res)	# sample sum of all deviations divided by n
#	rdmadm <- apply(dmadm,2, range)
#	return(list(gsadm=gsadm, sadm=sadm, rdmadm=rdmadm))
#}

c.one.group.norm.sub1.1 <- function(x, num.res, nc){
	# start values
	mat <- matrix(0,num.res,nc)
	rm <- rowMeans(x)		# row means over di.n,  num.res
	dm <- mat
	for(i in 1:nc){		# difference to row means
		dm[,i] <- x[,i] - rm
	}
	adm <- mat
	for(i in 1:nc){		# absolute differences
		adm[,i] <- abs(dm[,i])
	}
	sadm <- apply(adm,1, sum)	# sum per row
	gsadm <- sum(sadm)		# global sum
	rz <- apply(dm,2, function(x){ y <- range(x); z <- y[2]-y[1]; return(c(y,z)) })
	dimnames(rz)[[1]] <- c("min","max","diff")
	return(list(gsadm=gsadm, sadm=sadm, rz=rz))
}

c.one.group.norm.sub2 <- function(x, num.res, nc, gsadm.S, sampleN=0, rz){
	# sampling
	mat <- matrix(0,num.res,nc)
	gsadmAA <- gsadm.S
	gsadmA <- vector("numeric",sampleN)	# vector for all gsadmA
	fvec <- NULL	# gsadmA + factor sample vector -- if smaller than last gsadmAA
	for(k in 1:sampleN){
		f <- c.one.group.norm.sub3(rz)
		x1 <- mat
		for(i in 1:nc){
			x1[,i] <- x[,i] * f[i]
		}
		rm <- rowMeans(x1)		# row means over di.n
		dm <- mat
		for(i in 1:nc){		# difference to row means
			dm[,i] <- x1[,i] - rm
		}
		adm <- mat
		for(i in 1:nc){		# absolute differences
			adm[,i] <- abs(dm[,i])
		}
		sadm <- apply(adm,1, sum)	# sum per row
		gsadmA[k] <- sum(sadm)		# global sum
		if(gsadmA[k]<gsadmAA){
			fvec <- rbind( fvec, c(k, gsadmA[k], f) )	# k, global sum value, factors to the end
			gsadmAA <- gsadmA[k]
		}
	}
	#return(list(fvec=fvec, gsadm.S=gsadm.S, gsadmA=gsadmA, sampleN=sampleN))
	return(list(fvec=fvec, gsadm.S=gsadm.S, sampleN=sampleN))
}

c.one.group.norm.sub3 <- function(rz){
	# sample
	rzlen <- ncol(rz)
	f <- vector("numeric",rzlen)
	for(i in 1:rzlen){
		tmp <- runif(1,rz[1,i],rz[2,i])
		f[i] <- tmp/rz[3,i]
	}
	return(f)
}
#c.one.group.norm.sub3(matrix(c(-0.5,1,1.5, -0.4,1.1,1.7, 0.1,0.9,0.8),3,3,byrow=F))


# apply
c.one.group.norm.A <- function(x, factors){
	# apply optimal spread factors on expression data
	# x: matrix of expression values, samples in columns
	# factors: spread factors as many as columns in x and in an appropriate order
	nr <- nrow(x)
	nc <- ncol(x)
	# apply
	for(i in 1:nc){
		x[,i] <- x[,i] * factors[i]
	}
	return(x)
}
#a <- c.one.group.norm.A(aa, factors=a1list[[14]]$sample$fvec[14,3:length(a1list[[14]]$sample$fvec[14,])] )





## algo -4-  different approach  with  optimization - more consistent concerning factors

d.one.group.norm <- function(x, ssize=200, sampleN=0){
	# assumed one group for norm of expression data
	# x: matrix of expression values, samples in columns
	# ssize: sampling size is starting at > ssize
	# sampleN: number of samplings
	require(future)
	plan(multisession)
	
	x1 <- d.one.group.norm.sub0(x, ssize)	# check x
	# create distance grid
	if(x1$nr>ssize){	#index selection
		#idx.sel <- sample(x1$nr,ssize)
		idx.sel <- 1:200
		idxS <- T
	}else{
		idx.sel <- 1:x1$nr
		idxS <- F
	}
	nr <- x1$nr
	nc <- x1$nc
	x <- x1$x
	xlen <- length(idx.sel)
	rn <- dimnames(x)[[1]]
	cn <- dimnames(x)[[2]]
	# dist
	num.res <- all.neighbour.dist.sub1(xlen)	# length of result vector
	idx.p <- all.neighbour.dist.sub2(xlen)	# matrix, all possible index pairs of the ssize
	di <- matrix(0,num.res,nc)
	for(i in 1:nc){	# distance
		di[,i] <- all.neighbour.dist(x[idx.sel,i])	# signed!
	}
	# start values
	start <- d.one.group.norm.sub1(di, num.res, nc)
	#return(start)
	erg <- d.one.group.norm.sub2(di, num.res, nc, start$gsadm, sampleN, start$fdm)
	return(list(first=start,sample=erg))
}
#a <- d.one.group.norm(aa, ssize=200, sampleN=10)

d.one.group.norm.sub0 <- function(x, ssize){
	# check basics of x
	if(is.data.frame(x)){ x <- as.matrix(x) }
	if(anyNA(x)){ stop("NA in x not allowed") }		# alt.: stopifnot all expression T stop at first
	col.sums <- apply(x,2,sum)
	cat("\ncolumn sums",col.sums)
	nr <- nrow(x)
	col.zero <- apply(x,2, function(x){ sum(x==0) })
	cat("\ncolumn zero share",col.zero,"\nof length",nr,"\n")
	if(any(col.sums==0)){ stop("remove columns which have 0 counts") }
	nc <- ncol(x)
	if(nc<2){ stop("x: 2 and more columns") }
	if(ssize<3){ stop("ssize: 3 and more features") }
	return(list(x=x,nc=nc,nr=nr))
}
#d.one.group.norm.sub0(matrix(c(0,0,0,1,1,0,1,1,1),3,3,byrow=F), 3)

d.one.group.norm.sub1 <- function(x, num.res, nc){
	# start values
	mat <- matrix(0,num.res,nc)
	rm <- rowMeans(x)		# row means over di - length num.res
	dm <- mat
	for(i in 1:nc){		# difference to row means
		dm[,i] <- x[,i] - rm
	}
	fdm <- mat
	for(i in 1:num.res){		# create ideal absolute factors with the proper direction
		fdm[i,] <- rm[i] / x[i,]
	}
	for(i in 1:nc){		# correct  Inf -> 0  -Inf -> 0  -NaN -> 0  ##### check if appropriate
		fdm[fdm[,i]==Inf,i] <- 0
		fdm[fdm[,i]==-Inf,i] <- 0
		fdm[is.nan(fdm[,i]),i] <- 0
	}
	adm <- mat
	for(i in 1:nc){		# create global row var sum for start
		adm[,i] <- abs(dm[,i])
	}
	sadm <- apply(adm,1, sum)
	gsadm <- sum(sadm)
	return(list(gsadm=gsadm, sadm=sadm, fdm=fdm))
}
#a <- d.one.group.norm.sub1(matrix(c(1,3,7, 5,6,5, 2,0,0),3,3,byrow=T), 3, 3)

d.one.group.norm.sub2 <- function(x, num.res, nc, sum.dm, sampleN=0, fdm){
	# sampling
	mat <- matrix(0,num.res,nc)
	gsfdmAA <- sum.dm
	gsfdmA <- vector("numeric",sampleN)	# vector for all sampled gsadmA
	fvec <- NULL	# results: gsadmA + factor sample vector -- if smaller than last gsadmAA
	test <- NULL
	loof <- future({
	for(i in 1:sampleN){
		x1 <- mat
		for(k in 1:nc){	# sample factor for each row position and for each column
			f <- d.one.group.norm.sub3.1(fdm[,k])
			x1[,k] <- x[,k] * f
		}
		gsfdmA[i] <- d.one.group.norm.sub22(x1, num.res, nc)
		if(gsfdmA[i]<gsfdmAA){
			fvec <- rbind( fvec, c(i, gsfdmA[i], f) )	# k, global sum value, factors to the end
			gsfdmAA <- gsfdmA[i]
		}
	}
	})
	v <- value(loof)
	return(list(fvec=fvec, sum.dm=sum.dm, gsfdmA=gsfdmA, sampleN=sampleN))
}
#d.one.group.norm.sub2(x, num.res, nc, sum.dm, sampleN=0, fdm)

d.one.group.norm.sub22 <- function(x, num.res, nc){
	# global row var value
	mat <- matrix(0,num.res,nc)
	rm <- rowMeans(x)		# row means over di - length num.res
	dm <- mat
	for(i in 1:nc){		# difference to row means
		dm[,i] <- x[,i] - rm
	}
	# create global dm sum
	adm <- mat
	for(i in 1:nc){
		adm[,i] <- abs(dm[,i])
	}
	sadm <- apply(adm,1, sum)	# sum per row
	gsadm <- sum(sadm)		# global sum
	return(gsadm)
}

#d.one.group.norm.sub3 <- function(fdm.r){	# too many runif() calls...
#	# sample
#	fdmrlen <- length(fdm.r)
#	f <- vector("numeric",fdmrlen)
#	logiT <- fdm.r>1		# no if() runtime 11 min + 1.5 min worse
#	logiF <- fdm.r<1 & fdm.r!=0
#	ff <- NULL
#	for(i in fdm.r[logiT]){
#		ff <- c(ff, runif(1,1,i) )
#	}
#	f[logiT] <- ff
#	ff <- NULL
#	for(i in fdm.r[logiF]){
#		ff <- c(ff, runif(1,i,1) )
#	}
#	f[logiF] <- ff
#	return(f)
#}
#d.one.group.norm.sub3(matrix(c(0.5, 1, 1.5, 0.4, 0, 2),1,6))

d.one.group.norm.sub3.1 <- function(fdm.col){
	# sample
	fdmclen <- length(fdm.col)
	f <- vector("numeric",fdmclen)
	rfdmc <- range(fdm.col)
	rnd <- runif(fdmclen,rfdmc[1],rfdmc[2])		# rnd for overall range
	for(i in 1:fdmclen){	# adjust scale to individual target range  ((x-minS)*(maxT-minT)/(maxS-minS))+minT
		if(fdm.col[i]>1){
			f[i] <- ((fdm.col[i]-rfdmc[1])*(fdm.col[i]-1)/(rfdmc[2]-rfdmc[1]))+1
		}else if(fdm.col[i]<1 & fdm.col[i]>0){
			f[i] <- ((fdm.col[i]-rfdmc[1])*(1-fdm.col[i])/(rfdmc[2]-rfdmc[1]))+fdm.col[i]
		}
	}
	return(f)
}
#d.one.group.norm.sub3.1(matrix(c(0.5, 1, 1.5, 0.4, 0, 2),1,6))


# apply
d.one.group.norm.A <- function(x, factors){
	# apply optimal spread factors on expression data
	# x: matrix of expression values, samples in columns
	# factors: spread factors as many as columns in x and in an appropriate order
	nr <- nrow(x)
	nc <- ncol(x)
	# apply
	for(i in 1:nc){
		x[,i] <- x[,i] * factors[i]
	}
	return(x)
}
#a <- d.one.group.norm.A(aa, factors=a1list[[14]]$sample$fvec[14,3:length(a1list[[14]]$sample$fvec[14,])] )





#### further approach
#    norm by all possible pairwise d values between two groups
#    only good from 3 vs. 3 samples on (9 d [differences])
# estimate per feature line
# range of d-all: rd (min max)
# d-all / min(abs(d-all))		standard difference
# analysis of feature distributions

# over all features


# Two sample standardization over all pairwise difference values
# one feature
standardization.analysis.feature <- function(x){
	# standardization per feature 
	#  pairwise difference (d) based
	# x: a vector of differences per biological feature
	
	# all differences
	and <- all.neighbour.dist(x)		# fn: no redundancy, return: vec
	# order and visualize
	and.l <- length(and)
	and.r <- range(and)
	and.o <- and[order(and, decreasing=F)]
	plot(1:and.l, and.o, type="p", xlab="Ordered expression values per feature", ylab="Distance [Expression unit]")
	# test function creating an outlier limit ?

	return()
}
#standardization.analysis.feature(x=unlist(sal.singleend.unstranded[6,]))


#nrow(sal.singleend.unstranded)					# 244825
#length(unlist(sal.singleend.unstranded[6,]))	# 22
#range(unlist(sal.singleend.unstranded[6,]))		# 0.000 1238.577
#
#all.neighbour.dist.sub1(xlen=22)	# 231 : unique neighbour combinations in a vector of 22
#all.neighbour.dist.sub2(xlen=22)	# matrix : two columns with vector positions of combinations
	

# wrapper for all lines
standardization.analysis.samples <- function(x){
	# x: matrix or df of features in columns and samples in rows
	idf <- is.data.frame(x)
	nc <- ncol(x)
	nr <- nrow(x)
	d.f.len <- all.neighbour.dist.sub1(xlen=nc)
	erg <- matrix(0,nr,d.f.len)
	for(i in 1:nr){
		if(idf){ lin <- unlist(x[i,]) }else{ lin <- x[i,] }
		erg[i,] <- all.neighbour.dist(lin)
	}
	return(erg)
}
#a <- standardization.analysis.samples(x=sal.singleend.unstranded)
# 2 min 244825x231

#a1 <- date()
#pdf("test_all_dist.pdf",width=7,height=20)
#heatmap.ek.s(a[1:2000,], pc=F, pr=F, values=F, sign.num=2, ctext=NULL, cex.lab=1, move.cv=0, move.ch=0, rtext=NULL, move.rh=0, move.rv=0,
#		title.m=NULL, cex=1, cex.mat=0.7, mat.srt=0, grid=NA, diagC=F,
#		parmar=c(0.01,1,5,0.01), color.g=c("green","red"), col.miss.values="grey", center.col=NULL, special.col=NULL,
#		k.levels=10, custom.scale=NULL)
#dev.off()
#a1;date()	# min



# analysis of results
#standardization.analysis.conclusion <- function(x, file=NULL, main="", sub="", height=9){
#	# file: not NULL: string: "path/name.extension" is saving a graph if length.g1 <= 4
#	
#}





# options(width=180)
# save.image()