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



plot.sample.search.dependencies.2 <- function(x, raw=NULL, ref.size.range=c(1,4), ylim=NULL, y.at=NULL, ylab="correlation", text=T, description="",
			cex=0.8, mfrow=c(2,3), sym.norm=0.125, test.meta=F, ssqflag=1, regrl=T)
{
	# draw :  superimpose two regression plots (2 list element) and show expression strength
	#  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
	#  cex: text size,  mfrow: layout like mfrow philosophy (rows,cols)
	#  test.meta: T: only give meta data structure back
	#  sym.norm: needs adjustment if small test set, circle symbol size for expression values: scaling factor
	#   (! both experiments need to have the same scale)
	#   radius is given: 0.25 is maximum radius,  raw size is 1..4,   because area is quadratic, sqrt on expr. values 1..2 : -> /8
	#  regrl:T: plot regression line
	
	
	#ini
	if(missing(x)){ stop("\n correlation matrix, order vector and color must be provided") }
	len.list <- length(x)
	if(len.list<1 | len.list>2){ stop("\n list needs one or two list objects") }
	
	# source
	source(file="../0functions/0general/barplot.ek.R")
	
	# cex scaling
	cex.s <- 0.3
	
	# create meta data list structure (not ordered - ordering will take place in the sub functions)
	metadata <- vector("list", len.list)
	for(j in 1:len.list){
		metadata[[j]] <- vector("list", 6)
		names(metadata[[j]]) <- c("nr.test","nc.ref","no.rnames","no.cnames","mean.test","mean.ref")
	}
	# prefix 'no': not ordered (all)
	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]] <- dimnames(x[[j]][[1]])[[1]]		# rnames - get row names
		metadata[[j]][[4]] <- dimnames(x[[j]][[1]])[[2]]		# cnames - get col names
		if(!is.null(raw)){	# create mean expression values from 'raw' data
			metadata[[j]][[5]] <- colMeans( raw[[j]][ , metadata[[j]][[3]] ] )		# mean raw values for test , order of rnames
			metadata[[j]][[6]] <- colMeans( raw[[j]][ , metadata[[j]][[4]] ] )		# mean raw values for ref , order of cnames
		}
		# check
		cat("\n ref          ",j," - ",metadata[[j]][[4]] )
		cat("\n test         ",j," - ",metadata[[j]][[3]] )
		# check for test order (apply given order on given data)
		cat("\n test ordered ",j," - ",metadata[[j]][[3]][ x[[j]][[2]] ] )
	}
	# metadata out
	if(test.meta){ return(metadata) }
	
	# 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
	}
	
	# x axes offset in y direction
	ylim.off <- abs( (ylim[2]-ylim[1])/30 )
	
	# xlim
	xlim <- c(0.5,length(x[[1]][[2]]+0.5))
	
	
	# plot layout
	# height
	h.pl <- c(1,7,2)	# h.pl <- c(1,7,1.1,0.9)
	h.pl.l1 <- length(h.pl)
	h.pl <- rep(h.pl, mfrow[1])
	h.pl.l2 <- length(h.pl)
	# matrix
	m1 <- matrix(0, nrow=h.pl.l2, ncol=mfrow[2])
	ii <- 1
	for(i in 1:mfrow[1]){	# order to fill layout
		for(k in 1:mfrow[2]){
			for(j in ((1:h.pl.l1)+(i-1)*h.pl.l1)){
				m1[j,k] <- ii
				ii <- ii+1
			}
		}
	}
	# activate plot and fill according to imprinted order
	layout(m1, heights=h.pl)
	#layout.show(ii-1); return()
	
	par.tmp <- par(no.readonly=T)
	par(mar=c(3,4,3,2))		#default mar - lines: 5,4,4,2 +0.1 , ek 3,4,3,2
	
	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 all calculated global ssq
	
	for(i in 1:metadata[[1]][[2]]){		#create nc plots per list object
		
		## subplot 1 - name, strength, range - for every block
		par(mar=c(0.1,0.1,0.1,0.1))		#default mar - lines: 5,4,4,2 +0.1
		plot(x=1, y=1, type="n", ylab="", xlab="", xlim=xlim, ylim=c(0,2), tck=F, pch=1, axes=F, cex=cex, col=x[[1]][[3]])
		tmp.l <- 0.55	# y, work line
		# get graph area size
		tmp.r <- par()$usr
		tmp.r.u <- tmp.r[2]-tmp.r[1]
		# title text
		text(x=tmp.r.u*0.07, y=tmp.l, labels=metadata[[1]][[4]][i], srt=0, adj=0, cex=cex+4*cex.s, col="black")		#x[[2]][[3]]
		# expression strength (circle size)
		if(!is.null(raw)){
			# create mean expression input symbols
			#text(x=tmp.r.u*0.62, y=tmp.l, labels=paste("mean input:"), srt=0, adj=0, cex=cex+cex.s, col="black")
			#symbols(x=c(tmp.r.u*0.85,tmp.r.u*0.9), y=c(tmp.l,tmp.l), circles=sqrt(ref.size.range) *sym.norm, inches=F, bg="grey", fg=NULL, add=T)
			# create one or two expression reference labels
			for(j in 1:len.list){
				symbols(x=tmp.r.u*(0.3+(0.06*(j-1))), y=tmp.l, circles=sqrt(metadata[[j]][[6]][ i ]) *sym.norm, inches=F, bg=NULL, fg=x[[j]][[3]], add=T)	# switch bg,fg hollow,filled
			}
		}
		
		## subplot 2 - graph - for every block
		par(mar=c(3,4,3,2))
		for(j in 1:len.list){				#create one or two regression line(s) / x axis(axes) per nc
			#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 )
			}
			if(ssqflag==2){
				tmpV[j,i] <- linfit.ssq( x1, x[[j]][[1]][ x[[j]][[2]], i ] )
			}
			if(j==1){	# first regression
				if(is.null(raw)){	# create points according to mean expression values..
					plot(x=x1, y=x[[j]][[1]][ x[[j]][[2]], i ], ylab="", xlab="", xlim=xlim, ylim=ylim, tck=F, pch=1, axes=F, cex=cex, col=x[[j]][[3]])
				}else{	# ..or not
					plot(x=1, y=1, type="n", ylab="", xlab="", xlim=xlim, ylim=ylim, tck=F, pch=1, axes=F, cex=cex, col=x[[1]][[3]])
					symbols(x=x1, y=x[[j]][[1]][ x[[j]][[2]], i ], circles=sqrt(metadata[[j]][[5]][ x[[j]][[2]] ]) *sym.norm, inches=F, bg=NULL, fg=x[[j]][[3]], add=T)
				}# circles not*2
				# 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[[j]][[3]][ x[[j]][[2]] ])), y=ylim[1]-ylim.off*2.7, labels=metadata[[j]][[3]][ x[[j]][[2]] ], srt=45, adj=1, xpd=TRUE, cex=cex, col=x[[j]][[3]])
			}else{		# second regression
				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, cex=cex, col=x[[j]][[3]])
				}else{
					symbols(x=x1, y=x[[j]][[1]][ x[[j]][[2]], i ], circles=sqrt(metadata[[j]][[5]][ x[[j]][[2]] ]) *sym.norm, inches=F, bg=NULL, fg=x[[j]][[3]], 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, col=x[[j]][[3]])
			}
			loess <- lm(x[[j]][[1]][ x[[j]][[2]], i ]~x1)			# create regression line
			if(regrl){ abline(coef=loess$coefficients, col=x[[j]][[3]]) }		# plot regression line
			sasgc <- round(abs(loess$coefficients[2]),3)			# sasgc : lin.model coef. (slope): differential dependency
			# sasgc
			text(x=((xlim[2]-xlim[1])/3), y=if(j==1){ ylim[1]+((ylim[2]-ylim[1])*0.05) }else{ ylim[2]-((ylim[2]-ylim[1])*0.05) },
					labels=paste("lin.model coef. ",sasgc,sep=""), srt=0, adj=0, cex=cex+cex.s, col=x[[j]][[3]])
		}
		
		# 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=x[[2]][[3]]) }
		# y axis label
		mtext(text=ylab, side=2, line=3.0, outer=F, srt=90, adj=0.5, cex=cex, col="black")
		
		
		## subplot 3 - ? - for every block #par(mar=c(0.2,4,0.1,2))	# area merged with subplot 4
		
		
		## subplot 4 - ssq global - only for last block
		par(mar=c(0.3,1,0.3,1))
		plot(x=1, y=1, type="n", ylab="", xlab="", xlim=xlim, ylim=c(0,2), tck=F, pch=1, axes=F, cex=cex, col=x[[1]][[3]])
		
		if(i == metadata[[1]][[2]]){
			# 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) }
				if(ssqflag==2){ssqg[j] <- round(sum(tmpV[j, 1:metadata[[1]][[2]] ]), digits=4) }
			}
			
			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!=""){
					text(x=1, y=1.3, labels=description, adj=0, cex=cex+cex.s)
				}else{
					text(x=1, y=1.3, labels=tmp, adj=0, cex=cex+cex.s)
				}
				if(len.list==1){
					text(x=1, y=0.3, labels=paste("ssqg=",ssqg[1]), col=x[[1]][[3]], adj=0, cex=cex+cex.s)
				}else{
					text(x=1.1, y=0.24, labels=paste("ssqg=",ssqg[1]), col=x[[1]][[3]], adj=0, cex=cex+cex.s)
					text(x=3, y=0.3, labels=paste("/  ",ssqg[2]), col=x[[2]][[3]], adj=0, cex=cex+cex.s)
				}
			}
		}
		
	}
	cat("\n  min. global ssq = ", paste(ssqg, sep=", "), "\n")
	par(mar=c(5,4,4,2))		#default mar
	return()
}



##two
#pdf(file=paste(getwd(),"/results/robert/test",".",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.2(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.2(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=0.8, mfrow=c(2,3), ssqflag=1 )
#dev.off()