# TODO: Add comment
# 
# Author: E.Korsching 9.12.2011, 10-2021
###############################################################################



plot.sample.search.dependencies <- function(x, raw=NULL, ref.size.range=c(1,4), ylim=NULL, y.at=NULL, ylab="correlation", text=T, description="",
		refnames=NULL, testnames=NULL, cex.t=0.8, cex.rt=1, mfrow=NULL, mar.b.f=0, ssqflag=1)
{
	# draw regression plots (1 list element)    [or  superimpose two regression plots (2 list element)]
	#  with the same marker panel in the reference and test set (test might be differ in size)
	#  (avoid special characters in col names)
	# using results from  enumeration.search.dep.boots.f90  to fill cor table, order vec  (example at the end)
	# input : x: list containing a list with: 1. cor table, 2. order vec, 3. plot color :
	#         list(list(cor 1, order 1, color 1), list(cor 2, order 2, color 2))
	#  raw: list with one or two raw data matrix (with column names identical to the cor table names) for superimposing average expression values
	#  ref.size.range: min max values of expression range -- check if appropriate --
	#  ylim: y axis range: c(min,max) given or generated from the data,  y.at: tick positions,  ylab: y label
	#  text: print data name and global ssq (T/F) and if given a description in next line
	#  refnames, testnames: according to the order in given data provide different names
	#         list(list(refnames), list(refnames)), list(list(testnames), list(testnames))
	#  cex.t: rest size, cex.rt: ref,test labels,   mfrow: layout of plot page -OR- NULL if external layouts are given
	#  mar.b.f: extend bottom inner margin (add something)
	
	#ini
	if(missing(x)){ cat("\n correlation matrix, order vector and color must be provided"); break }
	len.list <- length(x)
	if(len.list<1 | len.list>2){ cat("\n list needs one or two list objects"); break }
	
	# create meta data list structure
	metadata <- vector("list", len.list)
	for(j in 1:len.list){
		metadata[[j]] <- vector("list", 6)
	}
	
	tmp.r <- c(0,0)		# start value for range
	for(j in 1:len.list){
		tmp <- x[[j]]
		len.element.list <- length(tmp)
		if(len.element.list!=3){ cat("\n internal list element ",j," must contain on the\n first position a cor table, on the second a order object and on the third a color"); break }
		#create max range
		tmp.r1 <- range(x[[j]][[1]])[1]
		tmp.r2 <- range(x[[j]][[1]])[2]
		if( tmp.r[1]>tmp.r1 ){ tmp.r[1] <- tmp.r1 } 
		if( tmp.r[2]<tmp.r2 ){ tmp.r[2] <- tmp.r2 }
		#meta data from all list elements
		metadata[[j]][[1]] <- nrow(x[[j]][[1]])	# nr
		metadata[[j]][[2]] <- ncol(x[[j]][[1]])	# nc
		metadata[[j]][[3]] <- if(is.null(testnames)){ dimnames(x[[j]][[1]])[[1]] }else{ testnames[[j]][[1]] }		# rnames - get row names
		metadata[[j]][[4]] <- if(is.null(refnames)){ dimnames(x[[j]][[1]])[[2]] }else{ refnames[[j]][[1]] }		# cnames - get col names
		if(!is.null(raw)){	# create mean expression values
			metadata[[j]][[5]] <- colMeans( raw[[j]][ , metadata[[j]][[3]] ] )		# mean raw values for rnames
			metadata[[j]][[6]] <- colMeans( raw[[j]][ , metadata[[j]][[4]] ] )		# mean raw values for cnames
		}
		#check for order
		cat("\n order ",j," - ",metadata[[j]][[3]][ x[[j]][[2]] ] )
	}
	#show range values
	cat("\n calculated range : ", tmp.r, " given range : ", if(is.null(ylim)){ "--" }else{ ylim } )
	if(is.null(ylim)){			# if not given use calculated range values
		ylim <- round((tmp.r+c(-0.1,0.1)), 1)		# adjust calculated y range for design
	}
	# adjust y range  for spare space on top (text and expr graph)
	ylim.a <- round( (ylim+c(0, abs( (ylim[2]-ylim[1])/3 ) ) ), 1)
	
	# x axes offset in y direction
	ylim.off <- abs( (ylim[2]-ylim[1])/30 )
	
	
	#plot
	if(is.null(mfrow)){
		par(mar=c(3+mar.b.f,4,3,2))	#default mar , lines: 5,4,4,2 +0.1
	}else{
		par(mfrow=mfrow, mar=c(3+mar.b.f,4,3,2))	#default mar , lines: 5,4,4,2 +0.1
	}
	
	tmpV <- matrix(0,len.list,metadata[[1]][[2]])		# recalculate min ssq for each sub situation to recalculate min global ssq
	ssqg <- vector(mode="numeric",length=len.list)			# store one or two calculated global ssq
	
	for(i in 1:metadata[[1]][[2]]){		#create nc plots per list object
		for(j in 1:len.list){				#create one or two regression line(s) / x axis(axes) per nc
			# point symbol normalization factor
			sym.norm <- (metadata[[j]][[1]] / 10)*0.083		# 1/12=0.0833
	
			#ascending vector of length nr
			x1 <- (1:metadata[[j]][[1]])
			if(ssqflag==1){
				tmpV[j,i] <- linfit.ssq( x[[j]][[1]][ x[[j]][[2]], i ], x1 )		#####fit for ssq calculation (in x )
			}
			if(ssqflag==2){
				tmpV[j,i] <- linfit.ssq( x1, x[[j]][[1]][ x[[j]][[2]], i ] )		#####fit for ssq calculation (in y )
			}
			if(j==1){
				if(is.null(raw)){	# create points according to mean expression values or not
					plot(x=x1, y=x[[1]][[1]][ x[[1]][[2]], i ], ylab="", tck=F, xlab="", ylim=ylim.a, pch=1, axes=F, cex=cex.t, col=x[[1]][[3]])
				}else{
					plot(x=x1, y=x[[1]][[1]][ x[[1]][[2]], i ], type="n", ylab="", tck=F, xlab="", ylim=ylim.a, pch=1, axes=F, cex=cex.t, col=x[[1]][[3]])
					symbols(x=x1, y=x[[1]][[1]][ x[[1]][[2]], i ], circles=metadata[[1]][[5]] *sym.norm,		# /: raw size is 1..4
							inches=F, bg=x[[1]][[3]], fg=NULL, add=T)
				}
				# axis and labels
				axis(side=1, at=x1, pos=ylim[1]-ylim.off, labels=F, col=x[[j]][[3]])
				text(x=c(1:length(metadata[[1]][[3]][ x[[1]][[2]] ])), y=ylim[1]-ylim.off*2.7, labels=metadata[[1]][[3]][ x[[1]][[2]] ], srt=45, adj=1, xpd=TRUE, cex=cex.rt, col=x[[j]][[3]])
				# ref color first experiment
#				text(x=par()$usr[1]+0.1, y=ylim.a[2], labels=metadata[[1]][[4]][i], srt=0, adj=0, cex=cex.rt+0.2, col="black")
				if(!is.null(raw)){
					# create ref expr color min and max
					tmp.r <- par()$usr
					tmp.r.u <- (tmp.r[2]-tmp.r[1])/10
					symbols(x=c(tmp.r.u*8.4,tmp.r.u*9), y=c(ylim.a[2],ylim.a[2]), circles=ref.size.range *sym.norm,		# /: norm factor
							inches=F, bg="grey", fg=NULL, add=T)
					# create ref expr color one or not
					symbols(x=tmp.r.u*5, y=ylim.a[2], circles=metadata[[1]][[6]][ i ] *sym.norm,
							inches=F, bg=x[[1]][[3]], fg=NULL, add=T)
				}
			}else{
				if(is.null(raw)){	# create points according to mean expression values or not
					points(x=x1, y=x[[j]][[1]][ x[[j]][[2]], i ], pch=1, col=x[[j]][[3]])
				}else{
					symbols(x=x1, y=x[[j]][[1]][ x[[j]][[2]], i ], circles=metadata[[j]][[5]] *sym.norm,
							inches=F, bg=x[[j]][[3]], fg=NULL, add=T)
				}
				# axis and labels
				axis(side=3, at=x1, pos=ylim[2]+ylim.off, labels=F, col=x[[j]][[3]])
				text(x=c(1:length(metadata[[j]][[3]][ x[[j]][[2]] ])), y=ylim[2]+ylim.off*2.7, labels=metadata[[j]][[3]][ x[[j]][[2]] ], srt=45, adj=0, xpd=TRUE, cex=cex.rt, col=x[[j]][[3]])
				# ref color second experiment
				if(!is.null(raw)){
					symbols(x=tmp.r.u*6, y=ylim.a[2], circles=metadata[[j]][[6]][ i ] *sym.norm,
							inches=F, bg=x[[j]][[3]], fg=NULL, add=T)
				}
			}
			loess <- lm(x[[j]][[1]][ x[[j]][[2]], i ]~x1)			# create fit line
			abline(coef=loess$coefficients, col=x[[j]][[3]])		# plot fit line
			sasgc <- round(abs(loess$coefficients[2]),3)			# sasgc : slope as strength of global correlation
			# reference label & sasgc
			text(x=par()$usr[1]+0.1, y=ylim.a[2], labels=paste(metadata[[1]][[4]][i],sep=""), srt=0, adj=0, cex=cex.rt+0.2, col="black")
			text(x=par()$usr[1]+0.1, y=ylim.a[2]-((ylim.a[2]-ylim.a[1])*0.1), labels=paste("dependency coef. ",sasgc,sep=""), srt=0, adj=0, cex=cex.rt, col="black")
			}
		
		# y axis scale/tick marks
		if(is.null(y.at)){ axis(side=2, yaxp=c(ylim[1],ylim[2], round((ylim[2]-ylim[1]),1)*10 ), col="black") } else { axis(side=2, at=y.at, col="black") }
		# y axis label
		mtext(text=ylab, side=2, line=3.0, outer=F, srt=90, adj=0.5, cex=cex.t, col="black")	
	}
	
	# calculate global ssq (ssqg)
	for(j in 1:len.list){
		if(ssqflag==1){ssqg[j] <- round(sum(tmpV[j, 1:metadata[[1]][[2]] ]), digits=1) }	# ssqg - x
		if(ssqflag==2){ssqg[j] <- round(sum(tmpV[j, 1:metadata[[1]][[2]] ]), digits=4) }	# ssqg - y
	}
	
	if(text){		#text below the last graphic tile
		tmp <- paste( deparse(substitute(x[[]])), collapse="" )	#extract data name
		tmp.len <- nchar(tmp)
		if(tmp.len>35){ tmp <- paste(substr(tmp, start=1, stop=35), ".. ") }	#restrict length of name to 35
		if(description!=""){
			mtext(text=description, side=1, line=3.2, adj=0, cex=cex.t-0.2)
		}else{
			mtext(text=tmp, side=1, line=3.2, adj=0, cex=cex.t-0.2)
		}
		if(len.list==1){
			mtext(text=paste("ssqg=",ssqg[1]), side=1, line=3.8, col=x[[1]][[3]], adj=0, cex=cex.t-0.2)
		}else{
			mtext(text=paste("ssqg=",ssqg[1]), side=1, line=3.8, col=x[[1]][[3]], adj=0, cex=cex.t-0.2)
			mtext(text=paste(ssqg[2]), side=1, line=3.8, col=x[[2]][[3]], adj=0.35, cex=cex.t-0.2)		# adj  might be different
		}
	}
	
	cat("\n  min. global ssq = ", paste(ssqg, sep=", "), "\n")
	return()
}


## check given refnames, testnames
#rk.Rmale.2.d.1h[[3]]	# given order
#rk.Rmale.2.d.1h[[1]]	# given order

## one - plus give names
#pdf(file=paste(getwd(),"/results/test02",".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""),
#		width=11, height=7, onefile=T, title="MBC.perm.lumi2", pointsize=12)
#plot.sample.search.dependencies(x=list(
#				list(rk.Rmale.2.d.1h[[3]][1:9,,drop=F], rk.Rmale.2.d.1h[[1]][1,2:10], "black")
#		),
#		raw=NULL, ref.size.range=c(1,4),
#		ylim=NULL, y.at=NULL, ylab="correlation", text=T, description="",
#		refnames=list(list(c("R-Bl","R-Bg","R-Bz"))), testnames=list(list(c("T Bl","T Bg","T Bz","T Bk","T Be","T Bu","T Bd","T Br","T Bm"))),
#		cex.t=0.8, cex.rt=1, mfrow=c(2,3), ssqflag=1 )
#dev.off()

## one
#pdf(file=paste(getwd(),"/results/test01",".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""),
#		width=11, height=7, onefile=T, title="MBC.perm.lumi2", pointsize=12)
#plot.sample.search.dependencies(x=list(
#				list(rk.Rmale.2.d.1h[[3]][1:9,,drop=F], rk.Rmale.2.d.1h[[1]][1,2:10], "black")
#		),
#		raw=NULL, ref.size.range=c(1,4),
#		ylim=NULL, y.at=NULL, ylab="correlation", text=T, description="", cex.t=0.8, cex.rt=1, mfrow=c(2,3), ssqflag=1 )
#dev.off()

## two
#pdf(file=paste(getwd(),"/results/test00",".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""),
#		width=11, height=7, onefile=T, title="MBC.perm.lumi2", pointsize=12)
#plot.sample.search.dependencies(x=list(
#				list(rk.Rmale.2.d.1h[[3]][1:9,,drop=F], rk.Rmale.2.d.1h[[1]][1,2:10], "black"),
#				list(rk.Hfemale.c2.d.1j[[3]][1:9,,drop=F], o=rk.Hfemale.c2.d.1j[[1]][1,2:10], "blue")
#		),
#		raw=list(Rmale.2.d, Hfemale.c2.d), ref.size.range=c(1,4),
#		ylim=NULL, y.at=NULL, ylab="correlation", text=T, description="1h - MBC, 1j - FBC", cex=0.8, mfrow=c(2,3), ssqflag=1 )
#
#plot.sample.search.dependencies(x=list(
#				list(rk.Rmale.2.d.1i[[3]][1:7,,drop=F], o=rk.Rmale.2.d.1i[[1]][1,2:8], "black"),
#				list(rk.Hfemale.c2.d.1k[[3]][1:7,,drop=F], o=rk.Hfemale.c2.d.1k[[1]][1,2:8], "red")
#		),
#		raw=NULL, ref.size.range=c(1,4),
#		ylim=c(-0.4,0.4), y.at=NULL, ylab="correlation", text=T, description="1i - MBC, 1k - FBC", cex.t=0.8, cex.rt=1, mfrow=c(2,3), ssqflag=1 )
#dev.off()