# TODO: Add comment
# 
# Author: E.Korsching 03-2013, 11-2019, 01-2020, 06-2022, 2024
###############################################################################

# heatmap.ek		feature rich + with clustering
# heatmap.ek.s		simple version for par(mfrow) and less features (e.g. no clustering)




## heatmap  with layout() design

heatmap.ek <- function(x, clev=1, values=T, sign.num=2, ctext=NULL, cex.lab=1, move.cv=0, move.ch=0, rtext=NULL, move.rh=0.3, move.rv=0,
		asp=T, title.m=NULL, cex=1.5, cex.mat=1, mat.srt=0, grid="black",
		parmar=c(0.01,0.01,0.01,0.01), layoutW=c(2,1.5,7,1.6), layoutH=c(1,2,7,1), lWe=0.15, lHe=0.05,
		color.g=c("green","red"), col.miss.values="grey", center.col=NULL, special.col=NULL,
		key=T, kxc=0.3, key.v.offset=5, full.key=T, invers=FALSE, k.levels=10, k.sign=2, custom.scale=NULL, cex.key=1,
		cluster=0, cl.method="average", dissmat="cor", show.dgram=T){
	# create a color matrix with or w/o text labels
	#  missing value for clustering: col-wise replaced by median of col  (otherwise it will not work)
	# x: data matrix
	# clev: color levels (internally y) the basis is: the whole matrix: 1 --or special cases -- basis row wise: 2, column wise: 3
	# values: print values in the plot matrix
	# ctext: col label str vector,  rtext: row label str vector,  cex.lab: char size axis text
	#  move.cv, move.ch: column (c) offset vertical (v) or horizontal (r)
	#  move.rh move.rv: row (r) offset horizontal (h) or vertical (v)
	# title.m: Titel,  also for color key label: cex: character size
	## if NULL :  Error in plot.new()  :  figure margins too large
	## ... had to increase the image size, or decrease the resolution png(filename="myfile.png", res=150, width = 1000, height = 1000)
	# asp: NULL: no aspect ratio, T: aspect ratio 
	# cex.mat: text size in matrix, mat.srt: text rotation in matrix, sign.num: number of decimals
	# grid: color or NULL: gobal fg or NA: omit borders
	# parmar: margin variable for all layout elements,
	#  layoutW: layout columns width 4 relative values, layoutH: layout lines height 4 relative values
	#  lWe/lHe: layout values for empty areas - dependencies!
	# color.g: provide 2 or 3 colors for a color gradient [code is depending on this feature]
	#  for linear scales with center : e.g. [-..0..+] : c("red","white","green")
	#  for linear scales w/o a center :scale=NULL,  e.g. [0..max] : c("red","green")
	# col.miss.values: color of missing values: choose a color
	# center.col: value from the data range; will be placed in the center of the color gradient
	# special.col: if not NULL: one or more special colors and positions which override the positions in the color palette
	#  paired: [position(s),color(s)], position starts with 1..n : dependend on k.levels+1 !
	#  e.g. list(c(1),c("lightblue"))  or  list(c(1,5,6),c("lightblue","pink","black"))
	# key: T: show color key,  kxc: correct proportions of key squares
	#  key.v.offset: put space above and below key to adjust proportions
	#  full.key: T: all labels will be shown, F: only extreme labels will be shown
	#  k.levels: define number of color levels  [means n-1 cuts]
	#  custom.scale: set a custom range independent from the data : c(min,max)
	#  invers: TRUE: color gradient will be inverted
	#  k.sign: significant numbers in the key legend
	# cluster: agnes based, 0: no, 1: rows, 2: cols, 3: both,
	#   4: both, but two different distance methods, dissmat needs: c("eu","cor") [row,column]
	#   cl.method: see 'agnes'
	#   dissmat: cor, eu  (cor is converted: -1->1 1->0  dissimilarity, agnes.cl)
	#   show.dgram: T: plot dendrogram
	
	# source other functions
	#require(dendextend)
	source("../0functions/cluster/agnes.cl.R")
	
	if(k.levels<2){ stop("\n k.levels >=2") }
	
	# sub-functions
	is.even <- function(x){ return(x %% 2 == 0) }		#for integers only, returns T or F  or a vector of logicals
#	is.odd <- function(x){ return(x %% 2 != 0) }
	
	f.step <- function(x){
		# calculate grid of step.size mid points - linear : 1-2-3-4  <-1.5-2.5-3.5->
		if(!is.null(custom.scale)){
			start <- custom.scale[1]
			end <- custom.scale[2]
		}else{
			start <- min(x, na.rm=TRUE)
			end <- max(x, na.rm=TRUE)
		}
		cat("\n Data min: ",start," max: ",end,"\n")
		step.size <- (end-start)/(k.levels-1)		# because step parts & step+1 colors
		start.o <- start		#save
		
		# first + last bin = one full bin -- check if appropriate 
		start <- start+step.size/2		#initialization
		
		# key parameter
		key.p <- c(start, step.size, start.o, end)
		
		# create color gradient
		pal <- f.pal.forming(key.p)
		return(list(pal, key.p))
	}
	
	f.pal.forming <- function(key.p){
		# create a color gradient
		if(!is.null(center.col)){
			len.color.g <- length(color.g)
			if(len.color.g<=2){
				stop("\n choose an odd number of colors >= 3 -or- set center.col=NULL")
			}
			even.colors <- is.even(len.color.g)
			if(even.colors){
				stop("\n choose an odd number of colors : e.g. color.g=c('red','white','blue') -or- set center.col=NULL")
			}
			len.color.mult <- len.color.g %/% 2
			f.pal1 <- colorRampPalette(color.g[1:(len.color.mult+1)], space="rgb")		#create function 1
			f.pal2 <- colorRampPalette(color.g[(len.color.mult+1):len.color.g], space="rgb")		#create function 2
			
			# which bin is color center
			i <- T
			j <- 1		# counter for center
			while(i){
				if(center.col<key.p[1]){
					i <- F
				}else{
					j <- j+1
					key.p[1] <- key.p[1] + key.p[2]
				}
			}
			
			pal1 <- f.pal1(j)		#generate color plalette 1
			# create now one color more (+2 because k.levels+1) which will be deleted in the next line (to start the gradient in a correct way)
			pal2 <- f.pal2(k.levels-j+2)		#generate color plalette 2
			pal <- c(pal1,pal2[2:(k.levels-j+2)])	#join both
			
		}else{
			f.pal <- colorRampPalette(color.g, space="rgb")		#create function
			pal <- f.pal(k.levels+1)		#generate color plalette
		}
		
		if(!is.null(special.col)){
			for(i in 1:length(special.col[[1]])){
				pal[special.col[[1]][i]] <- special.col[[2]][i]
			}
		}
		
		# invert col order
		if(invers){ pal <- pal[length(pal):1] }
		
		return(pal)
	}
	
	f.val.to.col <- function(x,y,pal.key.p){
		# translate matrix values to color values,  centered
		start <- pal.key.p[[2]][3]+pal.key.p[[2]][2]/2		#re-initialization
		if(sum(x<start, na.rm=T)!=0){		#this approach is implicitely dealing with rounding errors at the start and end
			tmp <- x<start
			tmp[is.na(tmp)] <- F	#correct NA to FALSE
			y[tmp] <- pal.key.p[[1]][1]
		}
		for(i in 2:k.levels){
			old <- start
			start <- start + pal.key.p[[2]][2]
			if(i<k.levels){
				if(sum(x>=old & x<start, na.rm=T)!=0){
					tmp <- x>=old & x<start
					tmp[is.na(tmp)] <- F	#correct NA to FALSE
					y[tmp] <- pal.key.p[[1]][i]
				}
			}else{
				if(sum(x>=old, na.rm=T)!=0){		#last step
					tmp <- x>=old
					tmp[is.na(tmp)] <- F	#correct NA to FALSE
					y[tmp] <- pal.key.p[[1]][i]
				}
			}
		}
		return(y)
	}
	
	f.key <- function(pal.key.p){
		# create color key legend
		plot(0, 0, xlim=c(0, 4), ylim=c(k.levels+key.v.offset, 0-key.v.offset), type="n", axes=F, xlab="", ylab="")
		
		delta.step <- 1
		dyn.pos <- 1
		kx0 <- 0.5+kxc
		kxe <- 1-kxc
		t.a <- 1.7
		
		old <- pal.key.p[[2]][3]	# for label description
		start <- pal.key.p[[2]][3]+pal.key.p[[2]][2]/2		#re-initialization
		
		# top
		xb <- c(kx0, kx0+kxe, kx0+kxe, kx0)
		yb <- c(dyn.pos-delta.step/2, dyn.pos-delta.step/2, dyn.pos+delta.step/2, dyn.pos+delta.step/2)	
		polygon(xb, yb, density=-1, col=pal.key.p[[1]][1], border=NULL, lwd=0.1)
		text(x=t.a, y=dyn.pos, labels=paste(" >= ",round(old,k.sign)), cex=cex.key, srt=0, adj=0)
		
		old <- start
		start <- start+pal.key.p[[2]][2]
		dyn.pos <- dyn.pos+delta.step
		for(i in 2:(k.levels-1)){
			xb <- c(kx0, kx0+kxe, kx0+kxe, kx0)
			yb <- c(dyn.pos-delta.step/2, dyn.pos-delta.step/2, dyn.pos+delta.step/2, dyn.pos+delta.step/2)	
			polygon(xb, yb, density=-1, col=pal.key.p[[1]][i], border=NULL, lwd=0.1)
			if(full.key){
				text(x=t.a, y=dyn.pos, labels=paste(round(old,k.sign),", <",round(start,k.sign)), cex=cex.key, srt=0, adj=0)
			}
			old <- start
			start <- start+pal.key.p[[2]][2]
			dyn.pos <- dyn.pos+delta.step
		}
		
		# bottom
		xb <- c(kx0, kx0+kxe, kx0+kxe, kx0)
		yb <- c(dyn.pos-delta.step/2, dyn.pos-delta.step/2, dyn.pos+delta.step/2, dyn.pos+delta.step/2)	
		polygon(xb, yb, density=-1, col=pal.key.p[[1]][(k.levels+1)], border=NULL, lwd=0.1)
		text(x=t.a, y=dyn.pos, labels=paste(" <= ", round(pal.key.p[[2]][4],k.sign)), cex=cex.key, srt=0, adj=0)
		
		# missing values
		if(na.sum>0){
			dyn.pos <- dyn.pos+delta.step+1
			xb <- c(kx0, kx0+kxe, kx0+kxe, kx0)
			yb <- c(dyn.pos-delta.step/2, dyn.pos-delta.step/2, dyn.pos+delta.step/2, dyn.pos+delta.step/2)	
			polygon(xb, yb, density=-1, col=col.miss.values, border=NULL, lwd=0.1)
			text(x=t.a, y=dyn.pos, labels=paste("NA"), cex=cex.key, srt=0, adj=0)
		}
		return()
	}
	
	## main function
	nr <- nrow(x)
	nc <- ncol(x)
	
	# color: create y template for color values
	y <- x
	
	# color levels: adjust k.levels if even number of step
	even.k.levels <- is.even(k.levels)
	if(even.k.levels){
		k.levels <- k.levels+1		# +1 because centered , odd : no action
	}
	
	# assign color to y from values in x
	if(clev==1){		# whole matrix
		pal.key.p <- f.step(x)
		y <- f.val.to.col(x,y,pal.key.p)
	}
	if(clev==2){		# row wise
		for(i in 1:nr){
			pal.key.p <- f.step(x[i,,drop=F])
			y[i,] <- f.val.to.col(x[i,,drop=F],y[i,,drop=F],pal.key.p)
		}
	}
	if(clev==3){		# col wise
		for(i in 1:nc){
			pal.key.p <- f.step(x[,i,drop=F])
			y[,i] <- f.val.to.col(x[,i,drop=F],y[,i,drop=F],pal.key.p)
		}
	}
	
	# color: special handling of missing values
	na.s <- is.na(x)
	na.sum <- sum(na.s)
	if(na.sum>0){
		y[na.s] <- col.miss.values
		warning(paste("\n",na.sum," NA(s) in the data set - color set to: ",col.miss.values,"\n", sep=""))
		
		# replace missing values for clustering
		#  by col median assuming that cols are variables - check if appropriate!
		if(cluster!=0){
			x <- replace.NA.rc(x, rc="col", fix=0, with="median")
			warning(paste("\nreplaced missing values for clustering by column median -- check if appropriate!\n", sep=""))
		}
	}
	
	# clustering on x
	if(cluster!=0){		# nothing
		if(nr<2 & nc<2){ stop("\n clustering needs at least matrix with 2 rows/columns") }
		if(cluster==1){		# per row
			cl.rows <- agnes.cl(t(x), label=NULL, extension=T, ttext="", dissmat=dissmat, cl.method=cl.method, scl="n01", cex=0.6, plot=F, return=T)
			# shift to greater zero - xlm.r
			if(min(cl.rows$height) < 0){
				cl.rows$height <- cl.rows$height - min(cl.rows$height)
			}
			xlm.x <- (max(cl.rows$height) - min(cl.rows$height)) / 10
			xlm.r <- c(max(cl.rows$height) + xlm.x, min(cl.rows$height) - xlm.x)
			cl.rows <- as.dendrogram(cl.rows)
			cl.rows.o <- order.dendrogram(cl.rows)
			cat("\n row order : ",cl.rows.o)
			# change order
			x <- x[cl.rows.o,]
			y <- y[cl.rows.o,]
		}
		if(cluster==2){		# per column
			cl.cols <- agnes.cl(x, label=NULL, extension=T, ttext="", dissmat=dissmat, cl.method=cl.method, scl="n01", cex=0.6, plot=F, return=T)
			# shift to greater zero - ylm.r
			if(min(cl.cols$height) < 0){
				cl.cols$height <- cl.cols$height - min(cl.cols$height)
			}
			ylm.x <- (max(cl.cols$height) - min(cl.cols$height)) / 10
			ylm.r <- c(max(cl.cols$height) + ylm.x, min(cl.cols$height) - ylm.x)
			cl.cols <- as.dendrogram(cl.cols)
			cl.cols.o <- order.dendrogram(cl.cols)
			cat("\n col order : ",cl.cols.o)
			# change order
			x <- x[,cl.cols.o]
			y <- y[,cl.cols.o]
		}
		if(cluster==3){		# both, row (first) and column (second)
			cl.rows <- agnes.cl(t(x), label=NULL, extension=T, ttext="", dissmat=dissmat, cl.method=cl.method, scl="n01", cex=0.6, plot=F, return=T)
			# shift to greater zero - xlm.r
			if(min(cl.rows$height) < 0){
				cl.rows$height <- cl.rows$height - min(cl.rows$height)
			}
			xlm.x <- (max(cl.rows$height) - min(cl.rows$height)) / 10
			xlm.r <- c(max(cl.rows$height) + xlm.x, min(cl.rows$height) - xlm.x)
			cl.rows <- as.dendrogram(cl.rows)
			cl.rows.o <- order.dendrogram(cl.rows)
			cat("\n row order : ",cl.rows.o)
			# change order
			x <- x[cl.rows.o, ]
			y <- y[cl.rows.o, ]
			
			cl.cols <- agnes.cl(x, label=NULL, extension=T, ttext="", dissmat=dissmat, cl.method=cl.method, scl="n01", cex=0.6, plot=F, return=T)
			# shift to greater zero - ylm.r
			if(min(cl.cols$height) < 0){
				cl.cols$height <- cl.cols$height - min(cl.cols$height)
			}
			ylm.x <- (max(cl.cols$height) - min(cl.cols$height)) / 10
			ylm.r <- c(max(cl.cols$height) + ylm.x, min(cl.cols$height) - ylm.x)
			cl.cols <- as.dendrogram(cl.cols)
			cl.cols.o <- order.dendrogram(cl.cols)
			cat("\n col order : ",cl.cols.o)
			# change order
			x <- x[, cl.cols.o]
			y <- y[, cl.cols.o]
			# change order
			#x <- x[cl.rows.o, cl.cols.o]
			#y <- y[cl.rows.o, cl.cols.o]
			
		}
		if(cluster==4){		# both (like above), but two different distance measures
			cl.rows <- agnes.cl(t(x), label=NULL, extension=T, ttext="", dissmat=dissmat[1], cl.method=cl.method, scl="n01", cex=0.6, plot=F, return=T)
			# shift to greater zero - xlm.r
			if(min(cl.rows$height) < 0){
				cl.rows$height <- cl.rows$height - min(cl.rows$height)
			}
			xlm.x <- (max(cl.rows$height) - min(cl.rows$height)) / 10
			xlm.r <- c(max(cl.rows$height) + xlm.x, min(cl.rows$height) - xlm.x)
			cl.rows <- as.dendrogram(cl.rows)
			cl.rows.o <- order.dendrogram(cl.rows)
			cat("\n row order : ",cl.rows.o)
			
			cl.cols <- agnes.cl(x, label=NULL, extension=T, ttext="", dissmat=dissmat[2], cl.method=cl.method, scl="n01", cex=0.6, plot=F, return=T)
			# shift to greater zero - ylm.r
			if(min(cl.cols$height) < 0){
				cl.cols$height <- cl.cols$height - min(cl.cols$height)
			}
			ylm.x <- (max(cl.cols$height) - min(cl.cols$height)) / 10
			ylm.r <- c(max(cl.cols$height) + ylm.x, min(cl.cols$height) - ylm.x)
			cl.cols <- as.dendrogram(cl.cols)
			cl.cols.o <- order.dendrogram(cl.cols)
			cat("\n col order : ",cl.cols.o)
			# change order
			x <- x[cl.rows.o, cl.cols.o]
			y <- y[cl.rows.o, cl.cols.o]
		}
	}
	
	# generate the layout matrix
	# t  t  t  t	t title
	# 0  0  cn 0	cn column names
	# rn cl H  k	rn row names, cl cluster tree rows, H heatmap, k legend
	# 0  0  cl 0	cl cluster tree columns
	lmat <- rbind(  c(1,1,1,1),
					c(2,3,4,5),
					c(6,7,8,9),
					c(10,11,12,13) )
	# adjust the layout width/height values according to visibility
	# title.m
	if(is.null(title.m)){
		layoutH[1] <- lHe		# value for empty region
	}
	# ctext
	if(is.null(ctext)){
		layoutH[2] <- lHe
	}
	# rtext
	if(is.null(rtext)){
		layoutW[1] <- lWe
	}
	# cluster !rows
	if((cluster==0 | cluster==2) & !show.dgram){
		layoutW[2] <- lWe
	}
	# cluster !cols
	if((cluster==0 | cluster==1 ) & !show.dgram){
		layoutH[4] <- lHe
	}
	# key
	if(key==F){
		layoutW[4] <- lWe
	}
	# show the outcome
	cat("\n layout width ",layoutW," layout height ",layoutH,"\n")
	
	## start plotting
	layout(lmat, width=layoutW, height=layoutH)	
#	layout.show(13); return();
	
	# adjust aspect ratio
	if(is.null(asp)){
		aarx <- 0;  aary <- 0
	}else{
		if(nc>nr){ aarx <- 0;  aary <- nc-nr }
		if(nc<nr){ aarx <- nr-nc;  aary <- 0 }
		if(nc==nr){ aarx <- 0;  aary <- 0 }
	}
	
	# set margins
	par(mar=parmar)
	
	# title
	if(!is.null(title.m)){
		plot(0, 0, xlim=c(0, 10), ylim=c(0, 1), type="n", axes=F, xlab="", ylab="" )
		text(x=0.5, y=0.5, labels=title.m, cex=cex+0.2, adj=0)
	}else{
		plot.new()
	}
	
	# advance
	plot.new()
	plot.new()
	
	# column description
	if(!is.null(ctext)){
		plot(0, 0, xlim=c(0.1, nc+0.6+aarx), ylim=c(0, 1), type="n", axes=F, xlab="", ylab="" )
		if(cluster==2 | cluster==3 | cluster==4){ ctext <- ctext[cl.cols.o] }
		for(i in 1:nc){		# top
			text(x=i+move.ch, y=move.cv, ctext[i], cex=cex.lab, srt=65, adj=0)
		}
	}else{
		plot.new()
	}
	
	# advance
	plot.new()
	
	# row description
	if(!is.null(rtext)){
		plot(0, 0, xlim=c(0, 10), ylim=c(nr+0.6, -(0.1+aary)), type="n", axes=F, xlab="", ylab="" )
		if(cluster==1 | cluster==3 | cluster==4){ rtext <- rtext[cl.rows.o] }
		for(i in 1:nr){		# left
			text(x=move.rh, y=i+move.rv, rtext[i], cex=cex.lab, srt=0, adj=0)
		}
	}else{
		plot.new()
	}
	
	# cluster rows
	if((cluster==1 | cluster==3 | cluster==4) & show.dgram){
		plot(cl.rows, leaflab="none", yaxt="n", horiz=T, xlim=xlm.r, ylim=c(nr+0.6, -(0.1+aary)))
	}else{
		plot.new()
	}
	
	# heatmap
	plot(0, 0, xlim=c(0.1, nc+0.6+aarx), ylim=c(nr+0.6, -(0.1+aary)), type="n", axes=F, ann=F)		# ann: no xlab ylab
	for(i in 1:nc){
		for(j in 1:nr){
			xb <- c(i-0.5, i+0.5, i+0.5, i-0.5)
			yb <- c(j-0.5, j-0.5, j+0.5, j+0.5)	
			polygon(xb, yb, density=-1, col=y[j,i], border=grid, lwd=0.1)
			if(values){
				text(x=i, y=j, round(x[j,i],sign.num), cex=cex.mat, srt=mat.srt, adj=c(0.5,0.5))
			}
		}
	}
#	testplotarea()
	
	# color key (optional)
	if(key){
		f.key(pal.key.p)
	}else{
		plot.new()
	}
	
	# advance
	plot.new()
	plot.new()
	
	# cluster columns
	if((cluster==2 | cluster==3 | cluster==4) & show.dgram){
		plot(cl.cols, leaflab="none", yaxt="n", xlim=c(0.1, nc+0.6+aarx), ylim=ylm.r)
		#test
		
	}else{
		plot.new()
	}
	
	# advance
	plot.new()
	# fin plot
	
	# clean up
	par(mfrow=c(1,1))	# reset layout mode
	par(mar=c(5.1,4.1,4.1,2.1))
}



#test.heat <- function(outfile="test"){
#	require(mdendro)
#	# heatmap2: cor,cor,eu, heatmap3:eu,eu,eu
#	m <- matrix(c(1,2,2,4, 3,4,1,1, 2,3,3,3),4,3)
#	dimnames(m)[[1]] <- c("R11","R22","R33","R44")
#	dimnames(m)[[2]] <- c("Cxx","Cyy","Czz")
#	m
#	t(m)
#	
##dev.off()
#	
#	cor(m)[lower.tri(cor(m), diag=F)]
#	
#	
#	library(pheatmap)
#	
#	pdf(outfile, width=11, height=7)		# heatmap2: cor,cor,eu, heatmap3:eu,eu,eu
#	par(mfrow=c(2,2))
## agnes works on cols
## pur agnes
##a <- agnes(m, diss=, metric=, method=, )
#	
## agnes.cl ek		cor:column wise - working differnt in pheatmap()  (!)
## on dim 2
#	a <- agnes.cl(m, label=NULL, extension=T, ttext="", dissmat="eu", cl.method="average", scl="n01", cex=0.6, plot=F, return=T)
## eu cor
#	cat("\n ",a$height)
#	if(min(a$height) < 0){ a$height <- a$height - min(a$height) }
#	cat("\n ",a$height)
#	plot(as.dendrogram(a), yaxt="n", horiz=F, xlim=c(3+0.6, -0.1), ylim=c(-0.1, max(a$height)) )		#  leaflab="none",
#	
#	a <- agnes.cl(m, label=NULL, extension=T, ttext="", dissmat="eu", cl.method="average", scl="n01", cex=0.6, plot=F, return=T)
#	cat("\n ",a$height)
#	if(min(a$height) < 0){ a$height <- a$height - min(a$height) }
#	cat("\n ",a$height)
#	plot(as.dendrogram(a), yaxt="n", horiz=T, xlim=c(-0.1, max(a$height)), ylim=c(3+0.6, -0.1) )
#	
## on dim 1
#	a <- agnes.cl(t(m), label=NULL, extension=T, ttext="", dissmat="eu", cl.method="average", scl="n01", cex=0.6, plot=F, return=T)
#	cat("\n ",a$height)
#	if(min(a$height) < 0){ a$height <- a$height - min(a$height) }
#	cat("\n ",a$height)
#	plot(as.dendrogram(a), yaxt="n", horiz=F, xlim=c(4+0.6, -0.1), ylim=c(-0.1, max(a$height)) )		#  leaflab="none",
#	
#	a <- agnes.cl(t(m), label=NULL, extension=T, ttext="", dissmat="eu", cl.method="average", scl="n01", cex=0.6, plot=F, return=T)
#	cat("\n ",a$height)
#	if(min(a$height) < 0){ a$height <- a$height - min(a$height) }
#	cat("\n ",a$height)
#	plot(as.dendrogram(a), yaxt="n", horiz=T, xlim=c(-0.1, max(a$height)), ylim=c(4+0.6, -0.1) )
#	
#	par(mfrow=c(1,1))
#	
#	heatmap.ek(x=m, ctext=c("Cxx","Cyy","Czz"), rtext=c("R11","R22","R33","R44"), title.m="TTT", key=T, cluster=3, dissmat="eu", show.dgram=T, cex=2, cex.lab=2, cex.mat=1)
## eu cor
#	
#	pheatmap(m, clustering_distance_rows="euclidean", clustering_distance_cols="euclidean", clustering_method="average", color=colorRampPalette(c("green", "red"))(50))
## euclidean correlation
#	
#	heatmap(m, hclustfun=linkage, scale="none" )
## only euclidean
#	
#	dev.off()
#}
#
#test.heat("test_5.pdf")





## heatmap simple  --  no cluster, no layout() design, no color key

heatmap.ek.s <- function(x, clev=1, values=T, 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(1,1,5,1), color.g=c("green","red"), col.miss.values="grey", center.col=NULL, special.col=NULL,
		k.levels=10, custom.scale=NULL){
	# create a color matrix with or w/o text labels
	# x: data matrix
	# clev: color levels (internally y):  basis: the whole matrix: 1 --or special cases -- basis row wise: 2, column wise: 3
	# values: print values in the plot matrix
	# ctext: col label str vector,  rtext: row label str vector,  cex.lab: char size axis text
	#  move.cv, move.ch: column (c) offset vertical (v) or horizontal (r)
	#  move.rh move.rv: row (r) offset horizontal (h) or vertical (v)
	# title.m: Titel, cex: character size
	## if NULL : Error in plot.new() : figure margins too large
	## ... had to increase the image size, or decrease the resolution png(filename="myfile.png", res=150, width = 1000, height = 1000)
	# cex.mat: text size in matrix, mat.srt: text rotation in matrix, sign.num: number of decimals
	# grid: color or NULL: gobal fg or NA: omit borders
	#  diagC: F: no action, T: omit top.left-bottom.right diagonal colors and values
	# parmar: margin variable for all layout elements,
	# color.g: provide 2 or 3 colors for a color gradient [code is depending on this feature]
	#  for linear scales with center : e.g. [-..0..+] : c("red","white","green")
	#  for linear scales w/o a center :scale=NULL,  e.g. [0..max] : c("red","green")
	# col.miss.values: color of missing values: choose a color
	# center.col: value from the data range; will be placed in the center of the color gradient
	# special.col: if not NULL: one or more special colors and positions which override the positions in the color palette
	#  paired: [position(s),color(s)], position starts with 1..n : dependend on k.levels+1 !
	#  e.g. list(c(1),c("lightblue"))  or  list(c(1,5,6),c("lightblue","pink","black"))
	# k.levels: define number of color levels  [means n-1 cuts]
	#  custom.scale: set a custom range independent from the data : c(min,max)
	
	if(k.levels<2){ stop("\n k.levels >=2") }
	
	# sub-functions
	is.even <- function(x){ return(x %% 2 == 0) }		#for integers only, returns T or F  or a vector of logicals
	
	f.step <- function(x){
		# calculate grid of step.size mid points - linear : 1-2-3-4  <-1.5-2.5-3.5->
		if(!is.null(custom.scale)){
			start <- custom.scale[1]
			end <- custom.scale[2]
		}else{
			start <- min(x, na.rm=TRUE)
			end <- max(x, na.rm=TRUE)
		}
		#cat("\n Data min: ",start," max: ",end,"\n")
		step.size <- (end-start)/(k.levels-1)		# because step parts & step+1 colors
		start.o <- start		#save
		
		# first + last bin = one full bin -- check if appropriate 
		start <- start+step.size/2		#initialization
		
		# key parameter
		key.p <- c(start, step.size, start.o, end)
		
		# create color gradient
		pal <- f.pal.forming(key.p)
		return(list(pal, key.p))
	}
	
	f.pal.forming <- function(key.p){
		# create a color gradient
		if(!is.null(center.col)){
			len.color.g <- length(color.g)
			if(len.color.g<=2){
				stop("\n choose an odd number of colors >= 3 -or- set center.col=NULL")
			}
			even.colors <- is.even(len.color.g)
			if(even.colors){
				stop("\n choose an odd number of colors : e.g. color.g=c('red','white','blue') -or- set center.col=NULL")
			}
			len.color.mult <- len.color.g %/% 2
			f.pal1 <- colorRampPalette(color.g[1:(len.color.mult+1)], space="rgb")		#create function 1
			f.pal2 <- colorRampPalette(color.g[(len.color.mult+1):len.color.g], space="rgb")		#create function 2
			
			# which bin is color center
			i <- T
			j <- 1		# counter for center
			while(i){
				if(center.col<key.p[1]){
					i <- F
				}else{
					j <- j+1
					key.p[1] <- key.p[1] + key.p[2]
				}
			}
			
			pal1 <- f.pal1(j)		#generate color plalette 1
			# create now one color more (+2 because k.levels+1) which will be deleted in the next line (to start the gradient in a correct way)
			pal2 <- f.pal2(k.levels-j+2)		#generate color plalette 2
			pal <- c(pal1,pal2[2:(k.levels-j+2)])	#join both
			
		}else{
			f.pal <- colorRampPalette(color.g, space="rgb")		#create function
			pal <- f.pal(k.levels+1)		#generate color plalette
		}
		
		if(!is.null(special.col)){
			for(i in 1:length(special.col[[1]])){
				pal[special.col[[1]][i]] <- special.col[[2]][i]
			}
		}
		
		return(pal)
	}
	
	f.val.to.col <- function(x,y,pal.key.p){
		# translate matrix values to color values,  centered
		start <- pal.key.p[[2]][3]+pal.key.p[[2]][2]/2		#re-initialization
		if(sum(x<start, na.rm=T)!=0){		#this approach is implicitely dealing with rounding errors at the start and end
			tmp <- x<start
			tmp[is.na(tmp)] <- F	#correct NA to FALSE
			y[tmp] <- pal.key.p[[1]][1]
		}
		for(i in 2:k.levels){
			old <- start
			start <- start + pal.key.p[[2]][2]
			if(i<k.levels){
				if(sum(x>=old & x<start, na.rm=T)!=0){
					tmp <- x>=old & x<start
					tmp[is.na(tmp)] <- F	#correct NA to FALSE
					y[tmp] <- pal.key.p[[1]][i]
				}
			}else{
				if(sum(x>=old, na.rm=T)!=0){		#last step
					tmp <- x>=old
					tmp[is.na(tmp)] <- F	#correct NA to FALSE
					y[tmp] <- pal.key.p[[1]][i]
				}
			}
		}
		return(y)
	}
	
	## main
	nr <- nrow(x)
	nc <- ncol(x)
	rx <- range(x)
	y <- x				# create template
	
	# adjust k.levels if even number of step
	even.k.levels <- is.even(k.levels)
	if(even.k.levels){
		k.levels <- k.levels+1		# +1 because centered , odd : no action
	}
	
	# assign color to y from values in x
	if(clev==1){		# whole matrix
		pal.key.p <- f.step(x)
		y <- f.val.to.col(x,y,pal.key.p)
	}
	if(clev==2){		# row wise
		for(i in 1:nr){
			pal.key.p <- f.step(x[i,,drop=F])
			y[i,] <- f.val.to.col(x[i,,drop=F],y[i,,drop=F],pal.key.p)
		}
	}
	if(clev==3){		# col wise
		for(i in 1:nc){
			pal.key.p <- f.step(x[,i,drop=F])
			y[,i] <- f.val.to.col(x[,i,drop=F],y[,i,drop=F],pal.key.p)
		}
	}
	
	# color: special handling of missing values
	na.s <- is.na(x)
	na.sum <- sum(na.s)
	if(na.sum>0){
		y[na.s] <- col.miss.values
		warning(paste("\n",na.sum," NA(s) in the data set - color set to: ",col.miss.values,"\n", sep=""))
	}
	
	# special handling of diagonal values
	if(diagC){
		if(nc==nr){
			for(i in 1:nc){
				x[i,i] <- NA
				y[i,i] <- NA
			}
		}else{	# skip
			cat("\n diagonal handling not possible rows!=columns")
		}
	}
	
	# set margins, clip to figure region (margins now open for text)
	parback <- par(no.readonly=T)
	par(mar=parmar,xpd=T)
	
	# heatmap
	# adjust aspect ratio
	diff <- abs(nc-nr)
	if(nc>nr){ sx <- nc; sy <- nr+diff }
	if(nc<nr){ sx <- nc+diff; sy <- nr }
	if(nc==nr){ sx <- nc; sy <- nr }
	plot(0, 0, xlim=c(-1, sx), ylim=c(sy, -1), type="n", axes=F, ann=F)	# blank plot
	mtext(text=title.m, side=3, line=1, cex=cex.lab)
	mtext(text=paste("range : ",round(rx[1],1)," - ",round(rx[2],1),sep=""), side=3, line=0, cex=cex.lab)
	for(i in 1:nc){
		for(j in 1:nr){
			xb <- c(i-0.5, i+0.5, i+0.5, i-0.5)
			yb <- c(j-0.5, j-0.5, j+0.5, j+0.5)	
			polygon(xb, yb, density=-1, col=y[j,i], border=grid, lwd=0.2)
			if(values){
				text(x=i, y=j, round(x[j,i],sign.num), cex=cex.mat, srt=mat.srt, adj=c(0.5,0.5))
			}
		}
	}
	# column description
	if(!is.null(ctext)){
		for(i in 1:nc){		# top
			text(x=i+move.ch, y=move.cv, ctext[i], cex=cex.lab, srt=65, adj=0)
		}
	}
	# row description
	if(!is.null(rtext)){
		for(i in 1:nr){		# left
			text(x=move.rh, y=i+move.rv, rtext[i], cex=cex.lab, srt=0, adj=1)
		}
	}
	
	# clean up
	par(parback)
	return()
}


#m <- matrix(c(1,2,2,4, 3,4,1,1, 2,3,3,3),4,3)
#dimnames(m)[[1]] <- c("R11","R22","R33","R44")
#dimnames(m)[[2]] <- c("Cxx","Cyy","Czz")
#m
#
#heatmap.ek.s(x=t(m), clev=1, values=T, sign.num=2, ctext=dimnames(m)[[1]], cex.lab=1, rtext=dimnames(m)[[2]],
#		title.m="Title", cex=1, cex.mat=0.7, mat.srt=0, grid="black", diagC=F,
#		parmar=c(1,1,5,5), color.g=c("green","red"), col.miss.values="grey", center.col=NULL, special.col=NULL,
#		k.levels=10, custom.scale=NULL)