# TODO: Add comment
# 
# Author: E.Korsching  11-2025
###############################################################################



#### workflow functions & experiments on the way
# sim: similarity,  dis: dissimilarity/distance

#options(width=250)


# vector distance
eu.dis <- function(v1, v2, p=1){
	# Euclidean distance/dissimilarity
	# no NA
	if(length(v1)!=length(v2)){ stop("v1,v2 vector length differs") }
	d <- sqrt( sum( (v1-v2)^2 ) )
	if(p==1){		# dist (not norm)
		return(d)
	}
	if(p==2){		# sim norm
		s <- 1 / exp(d)		# to compare with cosine similarity
		return(s)
	}
	return()
}
#eu.dis(c(1,0,0), c(1,1,1), p=1)	# 1.41

# vector difference
v.diff <- function(v1, v2){
	# difference function: subtract and split in distance and direction
	if(length(v1)!=length(v2)){ stop("v1,v2 vector length differs") }
	z <- v1-v2
	sg <- sign(z)
	z <- abs(z)
	return(list(dist=z,sign=sg))
}
#v.diff(c(1,0,0), c(1,1,1))			#  0 1 1   0 -1 -1

# average vector
v.avg <- function(v1, v2){
	# average of two vectors
	if(length(v1)!=length(v2)){ stop("v1,v2 vector length differs") }
	z <- (v1+v2)/2
	return(z)
}
#v.avg(c(1,0,0), c(1,1,1))			# 1.0 0.5 0.5

# vector magnitude
v.mag <- function(v1){
	# vector magnitude
	z <- sqrt(sum(v1^2))
	return(z)
}
#v.mag(c(1,0,0))					# 1



# enhanced  abs()  function
abs.sign <- function(x){
	# split in value and sign
	sg <- sign(x)
	z <- abs(x)
	return(list(x=z,sign=sg))
}
#abs.sign(c(5,-4,9,0,3,-9))			# 5 4 9 0 3 9   1 -1  1  0  1 -1



cos.sim <- function(v1, v2, p=1){
	# Cosine similarity, angle no magnitude
	# no NA
	if(length(v1)!=length(v2)){ stop("v1,v2 vector length differs") }
	s <- sum(v1*v2) / ( sqrt( sum( v1^2 ) ) * sqrt( sum( v2^2 ) ) )
	if(p==1){
		return(s)		# cos value 0..1
	}
	if(p==2){
		d <- 1 - s
		return(d)
	}
}
#cos.sim(c(1,0,0), c(1,1,1), p=1)	# 0.577


pear.cor <- function(v1, v2, p=1){
	# Pearson correlation similarity
	# no NA
	if(length(v1)!=length(v2)){ stop("v1,v2 vector length differs") }
	if(p==3){	# analytical
		v1m <- v1-mean(v1)
		v2m <- v2-mean(v2)
		c <- sum(v1m*v2m)
		v1sqrt <- sqrt(sum((v1-mean(v1))^2))
		v2sqrt <- sqrt(sum((v2-mean(v2))^2))
		f <- v1sqrt*v2sqrt
		g <- c / f
		cat("\n v1m ",v1m," v2m ",v2m," sum(v1m*v2m) ",c," -- v1sqrt ",v1sqrt," v2sqrt ",v2sqrt," v1sqrt*v2sqrt ",f," c/f ",g,"\n")
		return()
	}
	ps <- sum((v1-mean(v1))*(v2-mean(v2))) / ( sqrt( sum((v1-mean(v1))^2) ) * sqrt( sum((v2-mean(v2))^2) ) )
	if(p==1){	# similarity -1.. 0 (!) ..1
		return(ps)
	}
	if(p==2){	# dissimilarity 0..1
		pd <- (1-ps)/2		# -1->1 , 0->0.5 (!) , 1->0
		return(pd)
	}
}
#pear.cor(c(1,0,0), c(1,1,1), p=1)	# NaN
#pear.cor(c(1,2,1), c(1,0,1), p=1)	# -1




## sample function for cor or other simplifying proximity functions
sample.proxy.condM <- function(v1, v2, p=1, jitt="", fnIn="", samplN=10, main="", res=F){
	# wrapper sampling - condensing proximity measures
	# jitt: str: jitter on "v1" or "v2" or "v1v2", jitter series calls always different
	# res==T: return results
	if(fnIn==""){ stop("Give a valid function name as a string") }
	if(jitt=="v1"|jitt=="v2"|jitt=="v1v2"){
		cat("\n jitt ",jitt,"\n")
	}else{
		stop("Give valid --jitt--")
	}
	fnOp <- get(fnIn)
	v1l <- length(v1)
	v2l <- length(v2)
	samplN <- samplN+1
	erg <- vector("numeric",samplN)
	v1m <- NULL
	for(i in 1:v1l){ v1m <- cbind(v1m,rep(v1[i],samplN)) }
	v2m <- NULL
	for(i in 1:v2l){ v2m <- cbind(v2m,rep(v2[i],samplN)) }
	erg[1] <- fnOp(v1=v1, v2=v2, p=p)
	if(jitt=="v1"){ v1m[(2:samplN),] <- apply(v1m[(2:samplN),],2,jitter) }
	if(jitt=="v2"){ v2m[(2:samplN),] <- apply(v2m[(2:samplN),],2,jitter) }
	if(jitt=="v1v2"){ v1m[(2:samplN),] <- apply(v1m[(2:samplN),],2,jitter); v2m[(2:samplN),] <- apply(v2m[(2:samplN),],2,jitter) }
	for(i in 1:samplN){
		erg[i] <- fnOp(v1=v1m[i,], v2=v2m[i,], p=p)
	}
	#print(erg)
	par(mfrow=c(2,2))
	hist.plot.color(erg[2:samplN], xcolor=F, y=NULL, bin.num=15, x.range=NULL, norm=F, n.factor=1, lty=3, lwd=1,
			xlab="red, original result &\njitter distribution", ylab="", xlim=NULL, y.max=NULL, x.at=NULL, y.at=NULL, ylog=F, bar.width=1, offset=0, digits=2,
			col.grad=c("red","blue"), col.f="gray", col.sh=10, col.b="blue", cex=1, h.title=main, add=F, rC=F)
	abline(v=erg[1], lwd=2, col="red")
	if(jitt=="v1"|jitt=="v1v2"){
	for(i in 1:v1l){
	hist.plot.color(v1m[(2:samplN),i], xcolor=F, y=NULL, bin.num=15, x.range=NULL, norm=F, n.factor=1, lty=3, lwd=1,
			xlab=paste("jitter:",jitt,sep=" "), ylab="", xlim=NULL, y.max=NULL, x.at=NULL, y.at=NULL, ylog=F, bar.width=1, offset=0, digits=2,
			col.grad=c("red","blue"), col.f="gray", col.sh=10, col.b="blue", cex=1, h.title="red, original v1 value(s)", add=F, rC=F)
	abline(v=v1m[1,i], lwd=2, col="red")}}
	if(jitt=="v2"|jitt=="v1v2"){
	for(i in 1:v2l){
	hist.plot.color(v2m[(2:samplN),i], xcolor=F, y=NULL, bin.num=15, x.range=NULL, norm=F, n.factor=1, lty=3, lwd=1,
			xlab=paste("jitter:",jitt,sep=" "), ylab="", xlim=NULL, y.max=NULL, x.at=NULL, y.at=NULL, ylog=F, bar.width=1, offset=0, digits=2,
			col.grad=c("red","blue"), col.f="gray", col.sh=10, col.b="blue", cex=1, h.title="red, original v2 value(s)", add=F, rC=F)
	abline(v=v2m[1,i], lwd=2, col="red")}}
	par(mfrow=c(1,1))
	if(res){ return(erg) }else{ return() }
}
# sample.proxy.condM(v1=c(1,0,0), v2=c(1.01,0.99,1.02), p=1, jitt="v1", fnIn="pear.cor", samplN=20, main=paste("v1 c(1,0,0), v2 c(1.01,0.99,1.02)"))
# a <- sample.proxy.condM(v1=c(1,0,0), v2=c(1.01,0.99,1.02), p=1, jitt="v1", fnIn="pear.cor", samplN=20, main=paste("v1 c(1,0,0), v2 c(1.01,0.99,1.02)"), res=T)




#a <- cov(x=c(1,2,3),y=c(9,16,9))

# from asbio package
asbio.r.bw <- function(x, y=NULL){
	# asbio package - midvariances, midcovariance, midcorrelation
	# x: result: the biweight midvariance
	# x,y: result: the biweight midvariances, midcovariance, midcorrelation (robust alternative to Pearson's r)
	U.i <- (x-median(x))/(9*qnorm(.75)*mad(x))
	a.i <- ifelse(U.i<=-1|U.i>=1,0,1)
	n <- nrow(as.matrix(x))
	nx <- sqrt(n)*sqrt(sum((a.i*((x-median(x))^2))*((1-U.i^2)^4)))
	dx <- abs(sum(a.i*(1-U.i^2)*(1-5*U.i^2)))
	S.xx <- (nx/dx)^2
	if(!is.null(y)){
		V.i <- (y-median(y))/(9*qnorm(.75)*mad(y))
		b.i <- ifelse(V.i<=-1|V.i>=1,0,1)
		ny <- sqrt(n)*sqrt(sum((b.i*((y-median(y))^2))*((1-V.i^2)^4)))
		dy <- abs(sum(b.i*(1-V.i^2)*(1-5*V.i^2)))
		S.yy <- (ny/dy)^2
		S.xy <- n*sum((a.i*(x-median(x)))*((1-U.i^2)^2)*(b.i*(y-median(y)))*((1-V.i^2)^2))/((sum((a.i*(1-U.i^2))*(1-5*U.i^2)))*(sum((b.i*(1-V.i^2))*(1-5*V.i^2))))
		R.xy <- S.xy/(sqrt(S.xx*S.yy))}
	if(is.null(y))res <- data.frame(S.xx=S.xx)
	if(!is.null(y))res <- data.frame(s.xx=S.xx,s.yy=S.yy,s.xy=S.xy,r.xy=R.xy)
	return(res)
}
#a <- rnorm(20); a1 <- rnorm(20)
#asbio.r.bw(x=a); var(a)
#asbio.r.bw(x=a, y=a1)




## wrapper
w.mn.mn <- function(x, fn=NULL){
	# wrapper - apply function - samples in columns
	# matrix in: n columns -> matrix out: n col
	# x: matrix
	if(is.null(fn)){ stop("fn : missing") }else{ dfn <- get(fn) }
	x <- as.matrix(x)
	nr <- nrow(x)
	nc <- ncol(x)
	rname <- dimnames(x)[[1]]
	cname <- dimnames(x)[[2]]
	mat <- matrix(0, nrow=nr, ncol=nc)
	dimnames(mat)[[1]] <- rname
	dimnames(mat)[[2]] <- cname
	for(i in 1:nc){
		mat[,i] <- dfn(x[,i])
	}
	return(mat)
}
#w.mn.mn(x[,1], fn="norm.xi.Zscore") sum(x[,1]-mean(x[,1])) jitter
#x <- matrix(sample(jitter(as.vector(x))),10,10)
#summary(x)

w.mn.vn1 <- function(x, fn=NULL){		# e.g. cor
	# wrapper - apply function - samples in columns
	# matrix in: n columns -> vector out: n-1
	# x: matrix/ data frame
	# fn: str: a function for calculating one vector from two input vectors position based
	if(is.null(fn)){ stop("fn : missing") }else{ dfn <- get(fn) }
	
	x <- as.matrix(x)
	nr <- nrow(x)
	nc <- ncol(x)
	rname <- dimnames(x)[[1]]
	cname <- dimnames(x)[[2]]
	
	# prox
	erg.len <- (nc^2-nc)/2
	# create name combinations
	vec1 <- vector("character",erg.len)
	vecP <- vector("numeric",erg.len)
	k <- 1
	for(i in 1:(nc-1)){		#y
		for(j in (i+1):nc){		#x
			vec1[k] <- paste(cname[i], cname[j], sep=".")
			k <- k+1
		}
	}
	names(vecP) <- vec1
	# prox measure
	k <- 1
	for(i in 1:(nc-1)){		#y
		for(j in (i+1):nc){		#x
			vecP[k] <- dfn(x[,i], x[,j])
			k <- k+1
		}
	}
	return(vecP)
}

w.mn.vn <- function(x, fn=NULL){
	# wrapper - apply function - samples in columns
	# matrix in: n columns -> vector out: n
	# x: matrix/ data frame
	# fn: str: a function for calculating one vector from one input vector position based
	if(is.null(fn)){ stop("fn : missing") }else{ dfn <- get(fn) }
	
	x <- as.matrix(x)
	nc <- ncol(x)
	cname <- dimnames(x)[[2]]
	
	# reduce
	vecM <- vector("numeric", nc)
	names(vecM) <- cname
	for(i in 1:nc){
		vecM[i] <- dfn(x[,i])
	}
	return(vecM)
}

w.mn.pmn1 <- function(x, fn=NULL){
	# wrapper - apply function - samples in columns
	# matrix in: n columns -> matrix out: n-1 (proximity)
	# x: matrix/ data frame
	# fn: str: a function for calculating a proximity matrix
	if(is.null(fn)){ stop("fn : missing") }else{ dfn <- get(fn) }
	
	x <- as.matrix(x)
	nr <- nrow(x)
	nc <- ncol(x)
	rname <- dimnames(x)[[1]]
	cname <- dimnames(x)[[2]]
	
	# prox
	erg.len <- (nc^2-nc)/2
	# prox measure
	matP <- matrix(0, nrow=nc, ncol=nc)	
	dimnames(matP)[[1]] <- cname
	dimnames(matP)[[2]] <- cname
	for(i in 1:(nc-1)){		#y
		for(j in (i+1):nc){		#x
			matP[i,j] <- dfn(x[,i], x[,j])
			matP[j,i] <- matP[i,j]		# full matrix
		}
	}
	return(matP)
}

w.mn.l.pmn1.smn1 <- function(x, fn=NULL){
	# wrapper - apply function - samples in columns
	# matrix in: n columns -> list with two matrices: n-1 (proximity and sign)
	# x: matrix/ data frame
	# fn: str: a function for calculating the proximity per vector position
	if(is.null(fn)){ stop("fn : missing") }else{ dfn <- get(fn) }
	
	x <- as.matrix(x)
	nr <- nrow(x)
	nc <- ncol(x)
	rname <- dimnames(x)[[1]]
	cname <- dimnames(x)[[2]]
	
	# prox
	erg.len <- (nc^2-nc)/2
	# create name combinations
	vec1 <- vector("character",erg.len)
	k <- 1
	for(i in 1:(nc-1)){		#y
		for(j in (i+1):nc){		#x
			vec1[k] <- paste(cname[i], cname[j], sep="-")
			k <- k+1
		}
	}
	# prox measure
	matV <- matrix(0, nrow=nr, ncol=erg.len)	
	dimnames(matV)[[1]] <- rname
	dimnames(matV)[[2]] <- vec1
	matS <- matV
	k <- 1
	for(i in 1:(nc-1)){		#y
		for(j in (i+1):nc){		#x
			vs <- dfn(x[,i], x[,j])
			matV[,k] <- vs[[1]]		# value
			matS[,k] <- vs[[2]]		# sign
			k <- k+1
		}
	}
	return(list(value=matV, sign=matS))
}



## dendrograms, distances

d.upper.tri.idx.value <- function(x){
	# distance matrix, upper triangle, line by line into three columns: index X, index Y and value
	# x: symmetric dist matrix
	nc <- ncol(x)
	# number of dist elements w/o diagonal
	e.len <- (nc^2-nc)/2
	v.ix <- vector("integer",e.len)
	v.iy <- vector("integer",e.len)
	v.val <- vector("numeric",e.len)
	k <- 1
	for(i in 1:(nc-1)){		# row
		for(j in (i+1):nc){		# column
			v.ix[k] <- i
			v.iy[k] <- j
			v.val[k] <- x[i,j]
			k <- k+1
		}
	}
	return(data.frame(v.val=v.val, v.ix=v.ix, v.iy=v.iy))
}
# a <- d.upper.tri.idx.value(a.02.cor)
# a1 <- d.upper.tri.idx.value(matrix(0,10,10))

d.hier.dfn1 <- function(x, link="mean"){
	# list of hierarchical cluster nodes
	# x: symmetric dist table
	# link: linkage method: min, max, mean[avg]
	xnc <- ncol(x)
	i.n <- dimnames(x)[[2]]
	node1 <- data.frame(matrix(0,nrow=(xnc-1),ncol=7))
	names(node1) <- c("nodeIdx","value","ix1","iy1","idxN1","idxN2","idxNew")
	y <- x
	diag(y) <- NA
	for(i in 1:(xnc-1)){
		idx <- which(y==min(y, na.rm=T), arr.ind=T)
		#cat("\n",row.names(idx),idx[2,])
		node1[i,1:6] <- c(i, y[idx[2,1],idx[2,2]], idx[2,], row.names(idx))	# upper
		if(i!=(xnc-1)){	# no shrinking from 2x2 to 1x1 necessary
			ylist <- d.hier.dfn1.shrink(df=y, a1=idx[2,1], a2=idx[2,2], rnc=(xnc-i), link=link, i=i)
			y <- ylist[[1]]
			node1[i,7] <- ylist[[2]]
		}
	}
	node1 <- transform(node1, value=as.numeric(value))
	return(list(nodes=node1,names=i.n))
}
#a1 <- d.hier.dfn1(a.02.cor)

d.hier.dfn1.shrink <- function(df, a1, a2, rnc, link, i){
	# part of d.hier.dfn1()
	# df: dist() table
	# a1,a2: the positions in a dist() matrix
	# rnc: reduce the symmetric dist() matrix size
	# link: linkage method
	# error in 'out' if 2x2 matrix (out[xn,])
	fn <- get(link)		# single l. min , complete l. max , average l. mean
	# fill coordinates
	xn <- c(1:(rnc+1))[-a2]
	# copy+adjust
	out <- df[,xn]
	out <- out[xn,]
	new <- paste("n",i,"i",a1,a2,collapse="",sep="")
	#cat("\nnew",new)
	df.n <- dimnames(out)[[2]]
	df.n[a1] <- new
	dimnames(out)[[1]] <- df.n
	dimnames(out)[[2]] <- df.n
	# calc
	c1 <- apply(df[,c(a1,a2)], 1, fn)
	c1 <- c1[-a2]	# second discard
	c1[a1] <- NA	# first to NA
	# overwrite
	out[,a1] <- c1
	out[a1,] <- c1
	return(list(y=out, new=new))
}
# d.hier.dfn1.shrink(df=a.02.cor, a1=2, a2=3, rnc=9, link="mean")

d.dend.merge <- function(x){
	# create a merge table for dendrogram
	# x: d.hier.dfn1() table 
	nr <- nrow(x$nodes)
	out <- matrix(0,nr,2)
	for(i in nr:1){
		d2 <- unlist(x$nodes[i,c(5,6)])
		for(j in 1:2){
			a1 <- which(x$names %in% d2[j], arr.ind=T)
			if(!identical(a1, integer(0))){
				out[i, j] <- -a1
			}else{
				a1 <- which(x$nodes[,7] %in% d2[j], arr.ind=T)
				if(!identical(a1, integer(0))){
					out[i, j] <- a1
				}
			}
		}
	}
	mode(out) <- "integer"
	return(out)
}
# a2 <- d.dend.merge(a1)

d.dend.order <- function(x){
	# create an order vector for dendrogram
	# x: d.dend.merge() table 
	nr <- nrow(x)
	out <- vector("integer",nr)
	tmp <- out
	k <- 1
	l <- T
	m <- 1
	end <- nr+1
	for(i in nr:1){
		back <- d.dend.order.recu(i=i, x=x, k=k, end=end, out=out)
		k <- back[[1]]
		out <- back[[2]]
		if(k>=end){ break }
	}
	out <- as.integer(out)
	return(list(k=k,order=out))
}
# a3 <- d.dend.order(a2)

d.dend.order.recu <- function(i, x, k, end, out){
	# recursive function to d.dend.order()
	if(k>=end){ return(list(k,out)) }
	d2 <- unlist(x[i,c(1,2)])
	for(j in 1:2){
		if(d2[j]<0){
			out[k] <- abs(d2[j])
			k <- k+1
		}else{
			back <- d.dend.order.recu(i=d2[j], x=x, k=k, end=end, out=out)
			k <- back[[1]]
			out <- back[[2]]
		}
	}
	return(list(k,out))
}

d.dend.height <- function(n, height){
	# adjust dendrogram object height
	# n: dendrogram object,  see also agnes.cl()
	if(is.leaf(n)){
		attr(n, "height") <- height
		#print(attributes(n))
	}
	return(n)
}
# a5a <- dendrapply(a5, d.dend.height, height=0)

d.dend.plot <- function(x, link="avg"){
	# wrapper function - for ek dendrogram functions
	if(link=="avg"){ lk <- "mean" }
	a1 <- d.hier.dfn1(x, link=lk)
	a2 <- d.dend.merge(a1)
	a3 <- d.dend.order(a2)
	# create hclust like object
	a4 <- list(merge=a2,height=a1$nodes$value,order=a3$order,labels=a1$names,method=link,dist.method="input")
	attr(a4, "class") <- "hclust"
	a5 <- as.dendrogram(a4)
	a5max <- max(a1$nodes$value)
	a5min <- min(a1$nodes$value)
	h1 <- a5min-((a5max-a5min)/8)
	a5a <- dendrapply(a5, d.dend.height, height=h1)
	plot(a5a, xlim=c(0,length(a1$names)), ylim=c(h1,a5max))
	title(sub=paste("d.hier.dfn1(), ",deparse(substitute(x)),", ",link,sep=""))
	return()
}

# tests
# require(cluster)

#a1 <- d.hier.dfn1(a.02.cor)
#a2 <- d.dend.merge(a1)
#a3 <- d.dend.order(a2)
## hclust like object for plot(as.dendrogram(x))
#a4 <- list(merge=a2,height=a1$nodes$value,order=a3$order,labels=a1$names,method="average",dist.method="euclidean") #call="d.hier.dfn1()"
#attr(a4, "class") <- "hclust"
#a5 <- as.dendrogram(a4)
#a5a <- dendrapply(a5, d.dend.height, height=min(a1$nodes$value))
#plot(a5a, xlim=c(0,length(a1$names)), ylim=c(min(a1$nodes$value),max(a1$nodes$value)))
#
#plot(agnes(a.02.cor, diss=T, method="average"),  which.plot=2, main="", cex.main=0.9)
#
#a6 <- hclust(dist(t(a.00), method="euclidean"), method="average")
#a7 <- as.dendrogram(a6)
#
#
#par(mfrow=c(2,2))
#plot(as.dendrogram(hclust(dist(t(a.00), method="euclidean"), method="average")))
#title(sub="hclust(), a.00, euclidean, avg")
#plot(agnes(a.03.dis, diss=T, method="average"), main="", which.plot=2)
#mtext("agnes(), avg",side=1,line=2,cex=0.8)
#d.dend.plot(a.03.dis, link="avg")
#par(mfrow=c(1,1))


#save.image()




#### dist wrapper functions

prox.col.wise <- function(x=NULL, delta.fn=NULL, norm.fn=NULL, adj.fn=NULL){
	# wrapper - apply proximity function - return a list: value matrix and sign matrix
	# column wise - the preservation of dimensionality is assumed
	# x: matrix/ data frame, samples are columns
	# delta.fn: str: a function for calculating vector 'differences'
	# norm.fn: str: optional, an existing function normalizing x per column
	# adj.fn: str: optional, an existing function normalizing over all columns
	if(is.null(x)){ stop("x: data frame mandatory") }
	if(is.null(delta.fn)){ stop("delta.fn: missing") }else{ dfn <- get(delta.fn) }
	if(!is.null(norm.fn)){ nfn <- get(norm.fn) }
	if(!is.null(adj.fn)){ afn <- get(adj.fn) }
	
	x <- as.matrix(x)
	nr <- nrow(x)
	nc <- ncol(x)
	rname <- dimnames(x)[[1]]
	cname <- dimnames(x)[[2]]
	
	# norm
	if(!is.null(norm.fn)){
		x <- nfn(x)
	}
	# all norm
	if(!is.null(adj.fn)){
		x <- afn(x)
	}
	# diff
	erg.len <- (nc^2-nc)/2
	# create name combinations
	vec1 <- vector("character",erg.len)
	k <- 1
	for(i in 1:(nc-1)){		#y
		for(j in (i+1):nc){		#x
			vec1[k] <- paste(cname[i], cname[j], sep="-")
			k <- k+1
		}
	}
	# create proxy measure
	matV <- matrix(0, nrow=nr, ncol=erg.len)	
	dimnames(matV)[[1]] <- rname
	dimnames(matV)[[2]] <- vec1
	matS <- matV
	k <- 1
	for(i in 1:(nc-1)){		#y
		for(j in (i+1):nc){		#x
			vs <- dfn(x[,i], x[,j])
			matV[,k] <- vs[[1]]		# value
			matS[,k] <- vs[[2]]		# sign
			k <- k+1
		}
	}
	return(list(matV=matV, matS=matS))
}

#a2 <- data.frame(a=sample.int(50,5,replace=T),
#		b=sample.int(50,5,replace=T),
#		c=sample.int(50,5,replace=T),
#		d=sample.int(50,5,replace=T),
#		e=sample.int(50,5,replace=T),
#		row.names=c("a1","a2","a3","a4","a5"), stringsAsFactors=F)
#b <- prox.col.wise(x=a2, delta.fn="v.diff", norm.fn=NULL, adj.fn=NULL)

# norm
#a3 <- b$matV[,1] / sqrt(sum(b$matV[,1]^2))			# check: always between -1..1 ?
#a3 <- b$matV[,1] / sum(b$matV[,1])				# 0..1
#sum(a3)			# 1


prox.col.wise2 <- function(x=NULL, diss.fn=NULL, norm.fn=NULL, adj.fn=NULL){
	# wrapper - apply proximity function - full matrix returned for mds.stats()
	# column wise - dimensionality one reduction is assumed
	# x: matrix/ data frame, samples are columns
	# diss.fn: str: a function for calculating vector 'dissimilarities'
	# norm.fn: str: optional, an existing function normalizing x per column
	# adj.fn: str: optional, an existing function normalizing over all columns
	if(is.null(x)){ stop("x: data frame mandatory") }
	if(is.null(diss.fn)){ stop("diss.fn: missing") }else{ dfn <- get(diss.fn) }
	if(!is.null(norm.fn)){ nfn <- get(norm.fn) }
	if(!is.null(adj.fn)){ afn <- get(adj.fn) }
	
	x <- as.matrix(x)
	nr <- nrow(x)
	nc <- ncol(x)
	rname <- dimnames(x)[[1]]
	cname <- dimnames(x)[[2]]
	
	# norm
	if(!is.null(norm.fn)){
		x <- nfn(x)
	}
	# all norm
	if(!is.null(adj.fn)){
		x <- afn(x)
	}
	
	# create name combinations
	matN <- matrix(0,nc,nc)
	for(i in 1:(nc-1)){		#y
		for(j in (i+1):nc){		#x
			matN[i,j] <- paste(cname[i], cname[j], sep="-")
		}
	}
	# create proxy measure
	matV <- matrix(0,nc,nc)	
	dimnames(matV)[[1]] <- cname
	dimnames(matV)[[2]] <- cname
	for(i in 1:(nc-1)){		#y
		for(j in (i+1):nc){		#x
			tmp <- dfn(x[,i], x[,j])
			matV[i,j] <- tmp
			matV[j,i] <- tmp
		}
	}
	#vec1 <- matV[upper.tri(mat1, diag=F)]		# from top to bottom from left to right -> vector
	return(matV)
}
#b <- prox.col.wise2(x=a3, diss.fn="cor", norm.fn=NULL, adj.fn=NULL)




#### use existing graph functions

mds.stats <- function(x, group, colors=NULL, pch=1, txt=T, sicol=NULL, coladj=T, main="", kc=0, kccol=c("blue","red"), lpos="topright", roti=0){
	# mds plot - stats package
	# x: a full symmetric matrix containing the dissimilarities -> result of: prox.col.wise2() or dist()
	# group, colors: char vector of group labels as long as row/column names (group labels, colors in accordance to the x names)
	#  if colors=NULL&sicol=NULL: colors are choosen, if sicol!=NULL: one color for all
	# pch: one number  or  analogous to 'group' a number vector as long as group  ,  txt:T: text labels, F: no labels 
	# coladj:T: adjust color, F: no adjust
	# kc: number of cluster center:0 non, kccol: the colors for the centers,  rot:0: nothing, 90,180,270 clockwise
	## note: insert asp=1, in plot() fn, to ensure Euclidean distances are represented correctly
	require(scales)
	source("../0functions/cluster/color.groups.R")
	nc <- ncol(x)
	loc <- cmdscale(d=x, k=2)
	rn <- row.names(loc)
	x <- loc[, 1]
	y <- loc[, 2]
	if(length(x)<=kc){ cat("\nkc too small, set to 0 "); kc <- 0 }
	ug <- unique(group)
	ug.len <- length(ug)
	if(is.null(colors)&is.null(sicol)){
		g.col <- hue_pal()(ug.len)
		sample.col <- color.by.id(xc=g.col, group)
	}else if(is.null(colors)&!(is.null(sicol))){
		sample.col <- rep(sicol,nc)
	}else if(!(is.null(colors))&is.null(sicol)){
		sample.col <- colors
	}
	if(coladj){ sample.col <- adjustcolor(sample.col, alpha.f=0.9) }
	if(length(pch)==1){ charcode <- rep(pch,nc) }else{ charcode <- pch }
	coldata <- data.frame(condition=group, scol=sample.col, charcode=charcode)
	rownames(coldata) <- rn
	transf.k <- mds.stats.tr(x.m=x, y.m=y, rota=roti)
	if(txt){
		plot(x=transf.k$x2, y=transf.k$y2, type="n", xlab="c1", ylab="c2", xlim=range(transf.k$x2)*1.1, ylim=range(transf.k$y2)*1.1, axes=F, main=main, cex=1.5)
		axis(side=1)
		axis(side=2)
		text(x=transf.k$x2, y=transf.k$y2, labels=rn, cex=0.8, col=coldata$scol)		#pos=3,offset=0.8
	}else{
		plot(x=transf.k$x2, y=transf.k$y2, type="p", xlab="c1", ylab="c2", xlim=range(transf.k$x2)*1.1, ylim=range(transf.k$y2)*1.1, axes=F, main=main, cex=1.5, col=coldata$scol, pch=coldata$charcode)
		axis(side=1)
		axis(side=2)
	}
	# legend
	uc <- unique(colors)
	if(length(pch)==1){ px <- rep(pch,ug.len) }else{ px <- unique(pch) }
	legend(lpos, legend=ug, col=uc, pch=px)
	# kmeans center
	if(kc>0){
		cl <- kmeans(loc, centers=kc)
		transf.k <- mds.stats.tr(x.m=cl$centers[,1], y.m=cl$centers[,2], rota=roti)
		points(x=transf.k$x2, y=transf.k$y2, col=kccol, pch=8, cex=1.5)
	}
	return(loc)
}
#a <- mds.stats(ab.07.disN, group=c(rep("a",20),rep("b",25)), pch=1, txt=T, sicol=NULL, coladj=T, main="Test data")

mds.stats.tr <- function(x.m, y.m, rota){
	# rotate by 90 degree clockwise
	if(rota==0){
		x2 <- x.m
		y2 <- y.m
	}else if(rota==90){
		x2 <- y.m
		y2 <- -x.m
	}else if(rota==180){
		x2 <- -x.m
		y2 <- -y.m
	}else if(rota==270){
		x2 <- -y.m
		y2 <- x.m
	}else{	# no action between
		x2 <- x.m
		y2 <- y.m
	}
	return(list(x2=x2,y2=y2))
}