# TODO: Add comment
# 
# Author: E.Korsching Sep 27, 2013, 2024
###############################################################################



#### boxplots different to R standard boxplot
#    -boxplot.flex  simple + flexible, 
#    -boxplot.ek   more features and can be added to a graph area


boxplot.flex <- function(x, iqr=c(0.25,0.75), flags=c(1,1,1,1,1,1), at=NULL, text=NULL, ylab="",
		width=1, density=-1, border="black", col="lightblue", col.sp=c("blue","green"), lwd=1, cex=0.7, tcex=1, adj=0.5, srt=0)
{
	# plot boxplots
	# x: a data.frame or matrix, and samples in columns
	#    afterwards a data structure is created with 8 data characteristics:
	#    min, min-wisker, min-box, median, mean, max-box, max-whisker, max  (matrix row 1-8, per column sample)
	# flags: decision for (1) or against (0) graph element, triggers also graph ylim : min, min-wisker, median, mean, max-whisker, max
	# iqr: inter quantil range,   flags: 1: yes 0: no  order see above
	# at: at which x values should each boxplot be placed (NULL or as many as nc - positive values)
	# text: NULL  or  as many as columns given
	# border: box border: col or NA,  col: NULL  or  one color  or  as long as x subobjects,  col.sp: c(median color, mean color)
	# ylab: y label,  adj: label offset position,  srt: label rotation
	# width: box width,  density: >=1 shade or -1 fill box
	# lwd: line size,  cex: size of symbols, tcex: text label size
	# normally whisker 1.5*IQR rsp. the last data value inside this range (standard definition); beyond: single points (outliers)
	nr <- nrow(x)
	nc <- ncol(x)
	# calculate boxplot parameter 1-8
	mi <- apply(x,2, min, rm.na=T)	
	ma <- apply(x,2, max, rm.na=T)
	me <- apply(x,2, mean, rm.na=T)
	i1 <- apply(x,2, quantile, probs=iqr[1])
	i2 <- apply(x,2, quantile, probs=iqr[2])
	id <- i2-i1
	md <- apply(x,2, quantile, probs=0.5)
	wa <- md+id
	wi <- md-id
	cat("\n wisker <0 ? ",sum(wi<=0)," lower whisker set to 0")
	wi[wi<=0] <- 0
	x1 <- rbind(min=mi,minw=wi,minb=i1,md=md,me=me,maxb=i2,maxw=wa,max=ma)
	#print(x1)
	
	# some trigger values
	g.max <- max(x1[8,])
	w.max <- max(x1[7,], na.rm=T)
	# define x positions if not given
	if(is.null(at)){
		at <- seq(1,nc)
		x.max <- nc + 0.5
	}else{
		x.max <- max(at) + 0.5
	}
	# customize color values
	if(is.null(col)){
		col.len <- 0
	}else{
		col.len <- length(col)
	}
	if(col.len>1 & col.len<nc){ stop("\n error: you need a color for each boxplot  or  one color  or  set color to NULL") }
	if(col.len==1 & nc>1){ col <- rep(col[1],nc) }
	if(col.len==0){ col <- rep("white",nc) }
	# customize text labels
	if(is.null(text)){
		text.len <- 0
	}else{
		text.len <- length(text)
	}
	if(text.len!=0 & text.len<nc){ stop("\n error: you need a label for each boxplot or set text to NULL") }
	if(text.len==0){ text <- as.character(seq(1,nc)) }
	
	# plot per column
	par(xpd=T)		# mgp=c(3,1,0) if axis lables in m-line
	if(flags[6]==0){ y.max <- w.max }else{ y.max <- g.max }
	
	y.min <- y.max/12
	plot(x=0, y=0, type="n", ylab=ylab, xlab="", xlim=c(0.6,x.max), ylim=c(-y.min, y.max), axes=F, cex=cex)
	
	for(i in 1:nc){
		x.min <- x1[1,i]	# min
		x.min.w <- x1[2,i]	# min-wisker
		x.min.b <- x1[3,i]	# min-box
		x.median <- x1[4,i]	# median
		x.mean <- x1[5,i]	# mean
		x.max.b <- x1[6,i]	# max-box
		x.max.w <- x1[7,i]	# max-whisker
		x.max <- x1[8,i]	# max
		
		# vertical box
		xb <- c(at[i]-width/2, at[i]+width/2, at[i]+width/2, at[i]-width/2)			# x1,x2,x2,x1  box
		yb <- c(x.min.b, x.min.b, x.max.b, x.max.b)									# y1,y1,y2,y2  box
		
		top.h <- c(at[i]-width/4, x.max.w, at[i]+width/4, x.max.w)					# x1,y1,x2,y2  whisker line
		top.v <- c(at[i], x.max.b, at[i], x.max.w)									# x1,y1,x2,y2  connector line
		bottom.v <- c(at[i], x.min.w, at[i], x.min.b)								# x1,y1,x2,y2  connector line
		bottom.h <- c(at[i]-width/4, x.min.w, at[i]+width/4, x.min.w)				# x1,y1,x2,y2  whisker line
		median.l <- c(at[i]-(width/2.1), x.median, at[i]+(width/2.1), x.median)		# x1,y1,x2,y2  median line
		mean.l <- c(at[i]-(width/2.1), x.mean, at[i]+(width/2.1), x.mean)			# x1,y1,x2,y2  mean line
		
		bottom.out <- c(at[i], x.min)					# min outlier
		top.out <- c(at[i], x.max)						# max outlier
		# plot
		if(flags[2]==1){
			segments(x0=bottom.v[1], y0=bottom.v[2], x1=bottom.v[3], y1=bottom.v[4], lty=2, lwd=lwd, col=border)		# whisker down
			segments(x0=bottom.h[1], y0=bottom.h[2], x1=bottom.h[3], y1=bottom.h[4], lty=1, lwd=lwd, col=border)
		}
		if(flags[5]==1){
			segments(x0=top.h[1], y0=top.h[2], x1=top.h[3], y1=top.h[4], lty=1, lwd=lwd, col=border)					# whisker up
			segments(x0=top.v[1], y0=top.v[2], x1=top.v[3], y1=top.v[4], lty=2, lwd=lwd, col=border)
		}
		polygon(xb, yb, density=density, col=col[i], border=border, lwd=lwd)				# box
		if(flags[3]==1){ segments(x0=median.l[1], y0=median.l[2], x1=median.l[3], y1=median.l[4], lwd=(lwd+0.5), lty=1, col=col.sp[1]) }	# median line
		if(flags[4]==1){ segments(x0=mean.l[1], y0=mean.l[2], x1=mean.l[3], y1=mean.l[4], lwd=(lwd+0.5), lty="11", col=col.sp[2]) }		# mean line
		if(flags[1]==1){ points(x=bottom.out[1], y=bottom.out[2], pch=2, cex=cex) }		# min
		if(flags[6]==1){ points(x=top.out[1], y=top.out[2], pch=6, cex=cex) }		# max
		text(x=i, y=-(y.min/2), labels=text[i], cex=tcex, adj=adj, srt=srt)			# x axis labels
	}
	axis( side=2, las=0, cex.axis=1.1, cex.lab=1.1 )
	
	# reset
	par(xpd=F)
	return()
}

#a <- boxplot.flex(x=aaa,
#		flags=c(1,1,1,1,1,1), at=NULL, text=NULL, ylab="freq", width=0.9, density=-1, border="black", col=NULL, col.sp=c("black","orange"), lwd=1, cex=0.7, tcex=1, adj=0.5, srt=0)




boxplot.ek <- function(x, z=c(0,1), range.outer.plot=NULL, at=NULL, text=NULL, text.pos="top", col=NULL,
		skip.tails=F, prob.marks=c(0.025, 0.25, 0.50, 0.75, 0.975),
		label=T, width=1, density=-1, border=T, col.b="black", lwd=2, cex=0.7, tcex=1, adj=0.5, srt=0, horiz=F)
{
	# plot boxplots in existing plot
	# x: list of vector(s), each vector gives a boxplot -or- a data.frame, and each column gives a boxplot
	# z: pair: c(shift factor, scaling factor)
	# range.outer.plot: c(min,max) of surrounding plot
	# at: at which x values should each boxplot be placed (NULL or as many as in x)
	# text: NULL  or  as long as x subobjects,   text.pos: "top" "bottom"
	# col: NULL  or  one color  or  as long as x subobjects
	# skip.tails: T: no outliers, F: show outliers,  prob.marks: definition of top whisker, 3rd quantil, median, 2cd quantil, bottom whisker
	# label: T: box labels,  adj: label offset position,  srt: label rotation
	# width: box width,  density: shade or fill box,  border: box border,  col.b: box border color
	# lwd: line size,  cex: size of symbols, tcex: text label size,  horiz: T: horizontal barplots, F: vertical barplots
	
	# check if whisker 1.5*IQR rsp. the last data value inside this range (standard definition); beyond: single points (outliers)
	
	# ini
	if(is.null(range.outer.plot)){ stop("\n surrounding plot range has to be given ") }
#	present.coordinates <- par("usr")
	
	# adjust differences in the data type
	if(is.data.frame(x)){
		x.len <- ncol(x)	# number of list elements
		
		x.list <- vector("list",0)
		for(i in 1:x.len){
			x.list[[i]] <- x[,i]
		}
		x <- x.list
	}else{
		x.len <- length(x)	# number of list elements
	}
	
	# define x positions if not given
	if(is.null(at)){ at <- seq(1,x.len,1) }
	
	# customize color values
	if(is.null(col)){
		col.len <- 0
	}else{
		col.len <- length(col)
	}
	if(col.len>1 & col.len<x.len){ cat("\n error: you need a color for each boxplot  or  one color  or  set color to NULL"); return() }
	if(col.len==1 & x.len>1){ col <- rep(col[1],x.len) }		# rep 'no text'
	if(col.len==0){ col <- rep("white",x.len) }		# rep 'no text'
	# customize text labels
	if(is.null(text)){
		text.len <- 0
	}else{
		text.len <- length(text)
	}
	if(text.len!=0 & text.len<x.len){ cat("\n error: you need a label for each boxplot or set text to NULL"); return() }
	if(text.len==0){ text <- rep("",x.len) }		# rep 'no text'
	
	# text offset for top/bottom label
#	if(skip.tails){ x.r <- c(x.quantile[1], x.quantile[5]) }else{ x.r <- range(tmp.x) }
	x.r <- range.outer.plot		#range(x)
	text.pos.off <- (x.r[2]-x.r[1])/8
	if(text.pos=="top"){
		text.pos.set <- x.r[2]+text.pos.off
	}else{
		text.pos.set <- x.r[1]-text.pos.off
	}
	# text offset for median value
	median.pos.off <- width/6
	cat("\n range x ",x.r," text.pos.off ",text.pos.off," text.pos.set ",text.pos.set)
	
	# plot per subobject
	par(adj=adj)
	for(i in 1:x.len){
		tmp.x <- x[[i]]
		x.mean <- mean(tmp.x)							# calculation of mean
		x.quantile <- quantile(tmp.x, prob.marks)		# calculation of limits according to defined probabilities
		cat("\n percentil values ",x.quantile,"\n")
		cat("\n z[1] ",z[1],"\n")
		top.out <- tmp.x[tmp.x>x.quantile[5]]			# top outlier
		bottom.out <- tmp.x[tmp.x<x.quantile[1]]		# bottom outlier
		top.out.len <- length(top.out)
		bottom.out.len <- length(bottom.out)
		
		if(!horiz){		# vertical
			xb <- c(at[i]-width/2, at[i]+width/2, at[i]+width/2, at[i]-width/2)										# x1,x2,x2,x1  box
			yb <- c(z[1]+x.quantile[2], z[1]+x.quantile[2], z[1]+x.quantile[4], z[1]+x.quantile[4]) *z[2]			# y1,y1,y2,y2  box
			#
			top.h <- c(at[i]-width/4, (z[1]+x.quantile[5]) *z[2], at[i]+width/4, (z[1]+x.quantile[5]) *z[2])		# x1,y1,x2,y2  whisker
			top.v <- c(at[i], (z[1]+x.quantile[4]) *z[2], at[i], (z[1]+x.quantile[5]) *z[2])						# x1,y1,x2,y2  line
			bottom.v <- c(at[i], (z[1]+x.quantile[1]) *z[2], at[i], (z[1]+x.quantile[2]) *z[2])						# x1,y1,x2,y2  line
			bottom.h <- c(at[i]-width/4, (z[1]+x.quantile[1]) *z[2], at[i]+width/4, (z[1]+x.quantile[1]) *z[2])		# x1,y1,x2,y2  whisker
			median.l <- c(at[i]-width/2, (z[1]+x.quantile[3]) *z[2], at[i]+width/2, (z[1]+x.quantile[3]) *z[2])		# x1,y1,x2,y2  median line
			#
			out.points.x <- rep(at[i], (top.out.len+bottom.out.len))	# x, top-bottom
			out.points.y <- (c(top.out, bottom.out)+z[1]) *z[2]			# y, top-bottom
			#
			text.x <- at[i]
			text.y <- (text.pos.set+z[1]) *z[2]
			if(text.pos=="top"){ text.adj <- 0.5 }else{ text.adj <- 0.5 }
			median.x <- at[i]+width/2 + median.pos.off
			median.y <- (z[1]+x.quantile[3]) *z[2]
			median.value <- round((z[1]+x.quantile[3]) *z[2],2)
			mean.x <- at[i]
			mean.y <- (z[1]+x.mean) *z[2]
		}else{		# horizontal
			yb <- c(at[i]-width/2, at[i]-width/2, at[i]+width/2, at[i]+width/2)										# y1,y1,y2,y2  box
			xb <- c(z[1]+x.quantile[2], z[1]+x.quantile[4], z[1]+x.quantile[4], z[1]+x.quantile[2]) *z[2]			# x1,x2,x2,x1  box
			#
			top.h <- c((z[1]+x.quantile[5]) *z[2], at[i]-width/4, (z[1]+x.quantile[5]) *z[2], at[i]+width/4)		# x1,y1,x2,y2  whisker
			top.v <- c((z[1]+x.quantile[4]) *z[2], at[i], (z[1]+x.quantile[5]) *z[2], at[i])						# x1,y1,x2,y2  line
			bottom.v <- c((z[1]+x.quantile[1]) *z[2], at[i], (z[1]+x.quantile[2]) *z[2], at[i])						# x1,y1,x2,y2  line
			bottom.h <- c((z[1]+x.quantile[1]) *z[2], at[i]-width/4, (z[1]+x.quantile[1]) *z[2], at[i]+width/4)		# x1,y1,x2,y2  whisker
			median.l <- c((z[1]+x.quantile[3]) *z[2], at[i]-width/2, (z[1]+x.quantile[3]) *z[2], at[i]+width/2)		# x1,y1,x2,y2  median line
			#
			out.points.y <- rep(at[i], (top.out.len+bottom.out.len))	# y, top-bottom
			out.points.x <- (c(top.out, bottom.out)+z[1]) *z[2]			# x, top-bottom
			#
			text.x <- (text.pos.set+z[1]) *z[2]
			text.y <- at[i]
			if(text.pos=="top"){ text.adj <- 0 }else{ text.adj <- 1 }
			median.x <- (z[1]+x.quantile[3]) *z[2]
			median.y <- at[i]+ width/2 + median.pos.off
			median.value <- round((z[1]+x.quantile[3]) *z[2],2)
			mean.x <- (z[1]+x.mean) *z[2]
			mean.y <- at[i]
		}
		# plot
		polygon(xb, yb, density=density, col=col[i], border=if(border){col.b}else{NA}, lwd=lwd)						# box
		segments(x0=top.h[1], y0=top.h[2], x1=top.h[3], y1=top.h[4], lty=1, lwd=lwd)								# whisker
		segments(x0=top.v[1], y0=top.v[2], x1=top.v[3], y1=top.v[4], lty=2, lwd=lwd)
		segments(x0=bottom.v[1], y0=bottom.v[2], x1=bottom.v[3], y1=bottom.v[4], lty=2, lwd=lwd)					# whisker
		segments(x0=bottom.h[1], y0=bottom.h[2], x1=bottom.h[3], y1=bottom.h[4], lty=1, lwd=lwd)
		segments(x0=median.l[1], y0=median.l[2], x1=median.l[3], y1=median.l[4], lwd=(lwd+2))						# median line
		if(label){ text(x=median.x, y=median.y, labels=paste(median.value,sep=""), cex=tcex, adj=0, srt=srt) }		# box labels
		if(!skip.tails){ points(x=out.points.x, y=out.points.y, pch=16, cex=cex) }
		points(x=mean.x, y=mean.y, pch=3, cex=cex*2)																# add 'mean' marker
		text(x=text.x, y=text.y, labels=text[i], cex=tcex, adj=text.adj, srt=srt)									# x axis labels
	}
	
	par(adj=0.5)
	return()
}


##
#plot(c(1,20),c(1,20),xlim=c(-1,4))
#boxplot.ek(x=list(c(1,1,1,2,3,4,4,4,4),c(3,4,5,6,4,19)),
#		z=c(1,2), range.outer.plot=c(0,1), at=NULL, text=c("value1","value2"), text.pos="top", col=c("blue","green"),
#		skip.tails=T, prob.marks=c(0.025, 0.25, 0.50, 0.75, 0.975),
#		label=T, width=0.7, density=-1, border=F, col.b="black", lwd=2, cex=0.7, adj=0.5, srt=45, horiz=F)

##plot(1,1, type="n", xlim=c(0,5), ylim=c(0,1300), xlab="", ylab="", axes=T)	#vertical
#plot(1,1, type="n", xlim=c(0,1300), ylim=c(0,5), xlab="", ylab="", axes=T)	#horizontal
#boxplot.ek(x=b.dh1[,3:6],
#		z=c(0,1), range.outer.plot=c(0,1), at=NULL, text=c("aaaa","bbbb","cc","dd"), text.pos="top", col=NULL,
#		skip.tails=F, prob.marks=c(0.025, 0.25, 0.50, 0.75, 0.975),
#		label=T, width=0.7, density=-1, border=T, col.b="black", lwd=2, cex=0.7, adj=0.5, srt=45, horiz=T)
##points(x=rep(3,8),y=b.dh1[,5], cex=1, pch=16)

# see boew1.R
#f3v <- function(aa, ann, nc, col, cex=1, tcex=0.7, text.pos="bottom"){
#	an <- deparse(substitute(aa))
#	
#	pdf(file=paste(getwd(),"/results/23102015/",an,".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""),
#			width=11, height=7, onefile=T, title=paste(an, sep=""), pointsize=12)
#	
#	plot(1,1, type="n", xlim=c(0,nc+1), ylim=c(-0.2,1.1), xlab="", ylab="", cex=cex, axes=F)	#vertical
#	mtext(text="ratio distribution of cells with positive markers", side=2, line=0.5, cex=cex, adj=1)
#	boxplot.ek(x=aa,
#			z=c(0,1), range.outer.plot=c(0,1), at=NULL, text=ann, text.pos=text.pos, col=col,
#			skip.tails=T, prob.marks=c(0.025, 0.25, 0.50, 0.75, 0.975),
#			label=F, width=0.7, density=-1, border=T, col.b="black", lwd=2, cex=0.7, tcex=tcex, adj=0.5, srt=0, horiz=F)
#	
#	axis(side=2, at=c(0,0.25,0.5,0.75,1), labels=T, tick=T, pos=0.2, lty="solid", cex.axis=cex, col=NULL)
#	dev.off()
#}
#
#f3v(aa=c.DCIS.ERK56K18, ann=c("ER, K5/6, K8/18", "ER, K8/18"), nc=length(names(c.DCIS.ERK56K18)), col=c("blue"), cex=1.5, tcex=1.5, text.pos="bottom" )