# TODO: Add comment
# 
# Author: E.Korsching 10.9.2009
###############################################################################



scale.it <- function(
	x,					# data.frame  or  list of vectors
	mask=NULL,			# logical data.frame with T : mask for normalization but keep in data set
	method="IQRm",		# normalization method - see below
	trim=1,				# trim fn: values beyond a certain quantil value are set to that value (is changing data!) 0..1, 1: no limit
	max.y=1,			# several fn: max value to which all data are scaled
	center=F,			# scale fn: for scale(x, center=T, scale=T)
	divided=F,			# scale fn: scale parameter
	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",
	before=F			# if the normalization fails, the before plot gets column names
	)
{
	# this normalization is only implemented for columns
	
	# ini
	method <- match.arg(method,c("shiftscale","scale","median","iqr75","IQRm","Bol","mad","sd","range","sum"))
	
	ch.df <- is.data.frame(x)
	if(ch.df){
		nc <- dim(x)[2]
	}else{
		nc <- length(x)
	}
	cnames <- names(x)
	x.name=deparse(substitute(x))
#	rnames <- row.names(x)
	
	# trim
	if(trim < 1){
		trim.fn <- function(x,trim,xmax){
			xmax <- max(x)
			x.up.lim <- quantile(x[x>0],trim)	# -not- appropriate for every situation !
			x[x > x.up.lim] <- x.up.lim
			nxmax <- max(x)
			cat("\nold max",xmax," new ",round(trim*100,2), "% trimmed max",round(nxmax,3))
			return(x)
		}
		if(ch.df){
			x <- apply(x,2,trim.fn,trim)
		}else{
			x <- lapply(x,trim.fn,trim)
		}
		cat("\n")
	}
	
	# mask data for normalization		TO BE TESTED
	if(!is.null(mask)){
		z <- x
		x <- vector("list",0)
		for(i in 1:nc){
			x[[i]] <- z[mask[,i],i]
		}
	}
	
	# 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 )
		if(before){
			axis( side = 1 , at=c(1:length(names(x))) , labels=names(x), tick=TRUE, las=2, srt=srt )
			mtext(text=paste(x.name," - raw distributions\n"), side=3, line=line, outer=F, adj=1)
			return()	# manually decided - normalization not working
		}else{
			if(!is.null(mask)){
				tmp <- paste(x.name," - ",method," normalization\n >> some data masked")
			}else{
				tmp <- paste(x.name," - ",method," normalization\n")
			}
			mtext(text=tmp, side=3, line=line, outer=F, adj=1)
		}
	}
	
	if(method=="shiftscale"){
		shift.scale <- function(x){ (x-min(x))/(max(x)-min(x)) }
		
		if(is.null(mask) & ch.df){
			x <- apply(x,2,shift.scale)*max.y		# subtract min divide by max-min	
		}else{
			x <- lapply(x,shift.scale)*max.y
		}
		xlabel <- paste("scaling: (Xi-min(Xi))/(max(Xi)-min(Xi))")
	}
	if(method=="mad"){
		mad.scale <- function(x){ (x-median(x))/(mad(x)) }
		
		if(is.null(mask) & ch.df){
			x <- apply(x,2,mad.scale)*max.y		# subtract median divide by MAD
		}else{
			x <- lapply(x,mad.scale)*max.y
		}
		xlabel <- paste("scaling:(Xi-median(Xi))/mad(Xi)")
	}
	if(method=="scale"){
		if(is.null(mask) & ch.df){
			x <- apply(x,2,scale,center,divided)*max.y		# subtract col means divide by sd
		}else{
			x <- lapply(x,scale,center,divided)*max.y
		}
		xlabel <- paste("scaling:(Xi-mean(Xi)) / std.dev.(Xi)")
	}
	if(method=="median"){
		if(is.null(mask) & ch.df){
			av.median <- mean(apply(x,2,median))
		}else{
			av.median <- lapply(lapply(x,median),mean)
		}
		if(max.y==1){ max.y <- av.median }	# -not- appropriate for every situation !
		cat("\ndata will be scaled to a new average median for all variables:",max.y)
		if(is.null(mask) & ch.df){
			x <- apply(x,2,function(x,max.y){(x/median(x))*max.y}, max.y)			# divide by median , multiply by new max
		}else{
			x <- lapply(x,function(x,max.y){(x/median(x))*max.y}, max.y)
		}
		xlabel <- paste("scaling:",round(max.y,2)," / median(Xi)  ( av. median=",round(av.median),")")
	}
	if(method=="iqr75"){
		#	IQR.by.basic.fn <- function(x, na.rm = FALSE){ diff(quantile(as.numeric(x), c(0.25, 0.75), na.rm = na.rm)) }
		
		if(is.null(mask) & ch.df){
			# interquantile range normalization and 0.75 percentil & max adjustment
			nr <- nrow(x)
			# scale the IRQs to the max IQR value  (stretch column wise)
			colIQR <- apply(x, 2, IQR, na.rm=T)		# use standard fn -> see IQR help
			divisor <- matrix(rep(colIQR/max(colIQR), nr), nrow=nr, byrow=T)	# same dim as x - in each col the appropriate constant
			x.adj <- x/divisor	# col ranges will be adjusted by the denominator
			# move the whole columns by their 0.75 percentil values to the max 0.75 percentil value  (column wise)
			coluIQR <- apply(x.adj, 2, quantile, probs=0.75)		# upper IQR value
			adjustment <- matrix(rep(max(coluIQR) - coluIQR, nr), nrow=nr, byrow=T)
			x <- x.adj + adjustment
		}else{
			# list data
			colIQR <- lapply(x, IQR, na.rm=T)
			divisor <- unlist(colIQR)/max(unlist(colIQR))
			for(i in 1:nc){
				x[[i]] <- x[[i]]/divisor[i]
			}
			coluIQR <- lapply(x, quantile, probs=0.75)		# upper IQR value
			adjustment <- max(coluIQR) - coluIQR
			for(i in 1:nc){
				x[[i]] <- x[[i]]+adjustment[i]
			}
		}
		xlabel <- paste("interquantile range normalization and 0.75 percentil & max adjustment")
	}
	if(method=="IQRm"){
		if(is.null(mask) & ch.df){
			# interquantile range normalization and median adjustment
			nr <- nrow(x)
			# scale the IRQs to the max IQR value  (stretch column wise)
			colIQR <- apply(x, 2, IQR, na.rm=T)		# use standard fn -> see IQR help
			divisor <- matrix(rep(colIQR/max(colIQR), nr), nrow=nr, byrow=T)	# same dim as x - in each col the appropriate constant
			x.adj <- x/divisor	#col ranges will be adjusted by the denominator
			# move the whole column by their median values to the max of the medians  (column wise)
			colMed <- apply(x.adj, 2, median, na.rm=T)
			adjustment <- matrix(rep(max(colMed) - colMed, nr), nrow=nr, byrow=T)
			x <- x.adj + adjustment
		}else{
			# list data
			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("interquantile range normalization and median adjustment")
	}
	if(method=="sd"){
		if(is.null(mask) & ch.df){
			x <- apply(x,2,function(x,max.y){(x/sd(x))*max.y}, max.y)		# divide by sd multiply by max
		}else{
			x <- lapply(x,function(x,max.y){(x/sd(x))*max.y}, max.y)
		}
		xlabel <- paste("scaling: Xi / std.dev.(Xi)")
	}
	if(method=="range"){
		range.fn <- function(x){ max(x)-min(x) }
		
		if(is.null(mask) & ch.df){
			x <- apply(x,2,function(x,max.y){(x/range.fn(x))*max.y}, max.y)		# divide by spread multiply by max
		}else{
			x <- lapply(x,function(x,max.y){(x/range.fn(x))*max.y}, max.y)
		}
		xlabel <- paste("scaling: 1 / (max(Xi)-min(Xi))")
	}
	if(method=="sum"){
		if(is.null(mask) & ch.df){
			sum.xi <- colSums(x)
			max.sum <- max(sum.xi)
			x <- apply(x,2,function(x,max.sum){x/(sum(x)/max.sum)}, max.sum)		# divide by max based scale factor
		}else{
			sum.xi <- vector("numeric",nc)
			for(i in 1:nc){
				sum.xi[i] <- sum(unlist(x[[i]]))
			}
			max.sum <- max(sum.xi)
			x <- lapply(x,function(x,max.sum){x/(sum(x)/max.sum)}, max.sum)
		}
		xlabel <- paste("scaling: sum(Xi)/max(sum(Xi))")
	}
	
	cat("\n Cols are normalized\n")
	if(is.null(mask) & 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")
		sum.na <- sum(is.na(x))
		cat("\n NA values ? : ",sum.na," [0:none,>0:yes] \n")
		if(sum.na>0){ cat("\n Prevent NA values in return variable - stop\n"); return(x) }
	}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")
	}
#	x <- data.frame(x,row.names=rnames, check.rows=FALSE)
#	names(x) <- cnames
	
	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(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)
	}
}


#scale.it(
#		x=data.frame(a=c(1,2,6,8,9,2,4,9,5,4,9),b=c(3,6,5,9,9,2,7,9,2,7,9)),					#data.frame
#		method="iqr",		#normalization method "shiftscale","scale","median","IQRm","mad","sd","range","sum"
#		plot.it=T,			#plot normalization graph?
#		save.graph=F,				#save normalization graph in file
#		file.name="normalization",	#name of file and/or data.frame
#		save.data=F,				#save data in workspace
#		cex=0.6,
#		col.f=c("blue","green"),
#		col.b="black"
#)