# 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(1,1), sym.norm=0.125, test.meta=F, ssqflag=1)
{
	# 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 of plot page (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
	
	
	#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 }
	
	# source
	source(file="../0functions/0general/barplot.ek.R")
	
	# cex scaling
	cex.s <- 0.3
	
	# create meta data list structure (not ordered)
	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")
	}
	# 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]] )
		# check for test order (apply given order on given data)
#		cat("\n test         ",j," - ",metadata[[j]][[3]] )
		cat("\n test ordered ",j," - ",metadata[[j]][[3]][ x[[j]][[2]] ] )
	}
	
	# create an assignment between second(!) and first test data set (if two are given) based on the not ordered (!) versions
	#  and order at the end according to first
	# both test data sets are assumed to have the same length and identical identifier (! throughout)
	if(len.list==2){
		diff.index.first.sec <- vector("numeric",length=metadata[[1]][[1]])
		for(i in 1:metadata[[1]][[1]]){
			diff.index.first.sec[i] <- match(x=metadata[[2]][[3]][i], table=metadata[[1]][[3]])
		}
		#show assignment index between first and second
#		cat("\n assignment (n.o.)  ",diff.index.first.sec)
		
		# difference (correcting order for the second, keeping first raw order)
		diff.first.sec <- matrix(0, nrow=metadata[[1]][[1]], ncol=metadata[[1]][[2]])
		for(i in 1:metadata[[1]][[2]]){
			diff.first.sec[,i] <- abs( x[[1]][[1]][,i] - x[[2]][[1]][diff.index.first.sec,i] )
		}
		# order according to first
		diff.first.sec <- diff.first.sec[ x[[1]][[2]], ]
#		cat("\n difference (o.)\n")
#		print(diff.first.sec)
		tmp.d <- max(diff.first.sec)
	}
	
	# 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
	}
	
	# test out
	if(test.meta){
		return(metadata)
	}
	
	# 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 ( like an advanced par(mfrow) )
	# height
	h.pl <- c(1,7,1.1,0.9)
	h.pl <- rep(h.pl,mfrow[1])
	# matrix
	m1 <- matrix(0,nrow=mfrow[1]*4,ncol=mfrow[2])
	ii <- 1
	for(i in 1:mfrow[1]){
		for(k in 1:mfrow[2]){
			for(j in ((1:4)+(i-1)*4)){
				m1[j,k] <- ii
				ii <- ii+1
			}
		}
	}
	
	# activate plot and fill according to imprinted order
	layout(m1,heights=h.pl)
	
	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
		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.18, y=tmp.l, labels=metadata[[1]][[4]][i], srt=0, adj=0, cex=cex+4*cex.s, col=x[[2]][[3]])
		# expression size or not
		if(!is.null(raw)){
			# create ref expr color min and max
			text(x=tmp.r.u*0.51, y=tmp.l, labels=paste("range of mean expression:"), srt=0, adj=0, cex=cex+cex.s, col=x[[2]][[3]])
			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)
			for(j in 1:len.list){				#create one or two expr ref labels and min-max label
				# create ref expr color
				symbols(x=tmp.r.u*(0.07+(0.06*(j-1))), y=tmp.l, circles=sqrt(metadata[[j]][[6]][ i ]) *sym.norm, inches=F, bg=x[[j]][[3]], fg=NULL, add=T)
			}
		}
		
		## subplot 2 - graph
		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 )		#####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){	# first regression
				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="", xlab="", xlim=xlim, ylim=ylim, tck=F, pch=1, axes=F, cex=cex, col=x[[1]][[3]])
				}else{
					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[[1]][[1]][ x[[1]][[2]], i ], circles=sqrt(metadata[[1]][[5]][ x[[1]][[2]] ]) *sym.norm, 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, 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=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, col=x[[j]][[3]])
			}
			loess <- lm(x[[j]][[1]][ x[[j]][[2]], i ]~x1)			# create regression line
			abline(coef=loess$coefficients, col=x[[j]][[3]])		# plot regression line
			sasgc <- round(abs(loess$coefficients[2]),3)			# sasgc : slope as strength of global correlation
			# sasgc
			text(x=((xlim[2]-xlim[1])/2), y=if(j==1){ ylim[2]-((ylim[2]-ylim[1])*0.05) }else{ ylim[1]+((ylim[2]-ylim[1])*0.05) },
					labels=paste("dependency 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=x[[2]][[3]]) } 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=x[[2]][[3]])
		
		
		## subplot 3 - difference (perspective from j=1)
		par(mar=c(0.2,4,0.1,2))
		plot(x=1, y=1, type="n", ylab="", xlab="", xlim=xlim, ylim=c(-0.35,tmp.d), tck=F, pch=1, axes=F, cex=cex, col=x[[2]][[3]])
		# plot and labels
		barplot.ek(x=x1, y=diff.first.sec[,i], col=x[[1]][[3]], col.b=NA, width=0.3, density=-1, lwd=0.5, adj=0.5, srt=0, horiz=F)#xlab="", tick=F,
		text(x=1, y=-0.35, labels="difference (cor.)", srt=0, adj=0, xpd=TRUE, cex=cex+cex.s, col=x[[1]][[3]])
		
		
		## subplot 4 - only for last graph
		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) }	# 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!=""){
					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()