# TODO: Add comment
# 
# Author: E.Korsching  pre 2009
###############################################################################



scale.iqrm <- function(
	x,					# data.frame  or  list of vectors
	plot.it=T,			# plot normalization graph?
	logscale="",			# y axis with log scale: "y"
	save.graph=F,				# save normalization graph in pdf file
	file.name="normalization",	# name of file and/or data.frame
	save.data=F,				# save data in workspace
	outlier=F,			# show outlier
	mfrow=c(2,1),
	mar=c(7,5,4,4),
	srt=90,
	cex=0.6,
	line=0,
	adj=1,
	col.f=c("blue","green"),
	col.b="black"
	)
{
	# this normalization is only implemented for columns
	# IQRm
	
	# ini
	
	ch.df <- is.data.frame(x)
	if(ch.df){
		nc <- ncol(x)
	}else{
		nc <- length(x)		# vectors in a list
	}
	cnames <- names(x)
	x.name=deparse(substitute(x))
#	rnames <- row.names(x)
	
	# ini2
	mar1 <- c(1,5,4,1)
	mar2 <- c(8,5,1,1)
	
	if(plot.it){		# before scaling
		#layout
		layout(
				matrix(c(1,2), nrow=2, ncol=1, byrow=TRUE),
				widths=c(34),				#relative proportions
				heights=c(9,13),				#relative proportions
				respect=TRUE
		)
#		layout.show(2)	#function to check the layout
		
		if(save.graph){
#			png(filename=paste(getwd(),"/results/",file.name,".png", sep=""),
#				width=800, height=800, units="px", pointsize=12, bg="white")
			pdf(file=paste(getwd(),"/",file.name,".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""),
				width=11, height=7, onefile=T, title="normalization", pointsize=12)
		}
		
		cex.sav <- par()$cex
		adj.sav <- par()$adj
		par(cex=cex,adj=adj)
		
		par(mar=mar1)
		if(logscale=="y"){
			if(ch.df){			# +1: to get the log scale running if '0'
				xx <- x+1
			}else{
				xx <- vector("list",0)
				for(i in 1:nc){
					xx[[i]] <- x[[i]]+1
				}
			}
		}else{
			xx <- x
		}
		
		boxplot(x=xx,				# wisker : 1.5 * IQR range top/bottom
			outline=outlier,
#			outwex= .4,
			notch=F,
			ylim=range(boxplot(xx,plot=F)$stats),
			ylab=if(logscale=="y"){"intensity (+1)"}else{"intensity"},
			cex.lab=1.1,
			boxwex= .7,
			log=logscale,
			plot=T,
			axes=F,
			col=col.f[1],
			border=col.b
		)
		axis( side=2, cex.axis=1.1, cex.lab=1.1 )
		tmp <- paste(x.name," - IQR median normalization\n")
		mtext(text=tmp, side=3, line=line, outer=F, adj=1)
	}
	
	# normalization		new version 2023-01
	if(ch.df){
		nr <- nrow(x)
		# Adjust IQ ranges to be the same as max of IQRs
		colIQR <- apply(x, 2, IQR, na.rm=T)		# use standard fn -> see IQR help
		cat("\n colIQR",colIQR)		# IQR of 0 will not work
		divisor <- colIQR/max(colIQR)		# for each col the appropriate constant
		for(i in 1:nc){
			x[,i] <- x[,i]/divisor[i]		#col ranges werden durch den Teiler angeglichen
		}
		
		# Adjust medians to be the same as max of medians
		colMed <- apply(x, 2, median, na.rm=T)
		adjustment <- max(colMed) - colMed
		for(i in 1:nc){
			x[,i] <- x[,i]+adjustment[i]
		}
	}else{
		# Adjust IQ ranges to be the same as max of IQRs
		colIQR <- lapply(x, IQR, na.rm=T)
		divisor <- unlist(colIQR)/max(unlist(colIQR))
		for(i in 1:nc){
			x[[i]] <- x[[i]]/divisor[i]
		}
		colMed <- vector("numeric",nc)
		for(i in 1:nc){
			colMed[i] <- median(unlist(x[[i]]))
		}
		adjustment <- max(colMed) - colMed
		for(i in 1:nc){
			x[[i]] <- x[[i]]+adjustment[i]
		}
	}
	xlabel <- paste("scaling: to interquartile range and median: X.adj=Xi * max(IQR(Xi))/IQR(Xi) , X.norm=X.adj * max(median(X.adj.i))/median(X.adj.i)")
	
	cat("\n Cols are normalized\n")
	if(ch.df){
		cat("\n the new values are in the range: ", round(range(x, na.rm=T)[1],2), " - ", round(range(x, na.rm=T)[2],2), "\n")
		cat("\n NA values ? : ",sum(is.na(x))," [0:none,>0:yes] \n")
	}else{
		sum.xi <- matrix(0,nc,2)
		for(i in 1:nc){
			sum.xi[i,] <- range(unlist(x[[i]]))
		}
		cat("\n the new values are in the range: ", round(min(sum.xi[,1]),2), " - ", round(max(sum.xi[,2]),2), "\n")
		na.xi <- vector("numeric",nc)
		for(i in 1:nc){
			na.xi[i] <- sum(is.na(unlist(x[[i]])))
		}
		cat("\n NA values ? : ",sum(na.xi)," [0:none,>0:yes] \n")
	}
	
	if(plot.it){		# after scaling
		par(mar=mar2)
		if(logscale=="y"){
			if(ch.df){			# +1: to get the log scale running if '0'
				xx <- x+1
			}else{
				xx <- vector("list",0)
				for(i in 1:nc){
					xx[[i]] <- x[[i]]+1
				}
			}
		}else{
			xx <- x
		}
		
		boxplot(x=xx,
			outline=outlier,
#			outwex= .4,
			notch=F,
			ylim=range(boxplot(xx,plot=F)$stats),
			ylab=if(logscale=="y"){"intensity (+1)"}else{"intensity"},
			cex.lab=1.1,
			boxwex= .7,
			log=logscale,
			plot=T,
			axes=F,
			col=col.f[2],
			border=col.b
		)
		axis( side = 1 , at=c(1:length(names(x))) , labels=names(x), tick=TRUE, las=2, srt=srt )
		axis( side = 2 , cex.axis=1.1, cex.lab=1.1)
		
		if(save.graph){	dev.off() }
	}
	
	# restore
	par(mfrow=c(1,1), mar=c(5,4,4,2), cex=cex.sav, adj=adj.sav)
	# save file and data set in workspace   or   return data set
	if(save.data){
		if(is.null(file.name)|file.name==""){ cat("\n No data name \n"); break }
		assign(file.name, x , where=1)
		cat("\n Normalized data is saved in ",file.name,"\n")
	}else{
	    return(x)
	}
}


#aa <- scale.iqrm(
#	x=data.frame(cbind(c(1,5,7,10,9,100),c(3,4,6,10,8,90))),	# data.frame  or  list of vectors
##	x=cbind(runif(10,1,10),runif(10,1,10)),	# data.frame  or  list of vectors
#	plot.it=T,			# plot normalization graph?
#	logscale="",			# y axis with log scale: "y"
#	save.graph=F,				# save normalization graph in pdf file
#	file.name="normalization",	# name of file and/or data.frame
#	save.data=F,				# save data in workspace
#	outlier=F,			# show outlier
#	mfrow=c(2,1),
#	mar=c(7,5,4,4),
#	srt=90,
#	cex=0.6,
#	line=0,
#	adj=1,
#	col.f=c("blue","green"),
#	col.b="black"
#	)