#
#   distAssoc package
#   Copyright (C) 2005  Jan de Leeuw <deleeuw@stat.ucla.edu>
#   UCLA Department of Statistics, Box 951554, Los Angeles, CA 90095-1554
#   
#   This program is free software; you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation; either version 2 of the License, or
#   (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#   
#   You should have received a copy of the GNU General Public License
#   along with this program; if not, write to the Free Software
#   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
###################################################################
#
# version 0.0.01, 2005-11-13
# version 0.0.02, 2005-11-14 (this one actually works)
# version 0.0.03, 2005-11-15 (inner products added)
# version 0.0.04, 2005-11-15 (single mode analysis added)
# version 0.0.05, 2005-11-15 (changed deviance to chi-square)
# version 0.0.06, 2005-11-15 (separate bias estimation with IPF)
# version 0.0.07, 2005-11-15 (separate versions for one and two modes)
# version 0.0.08, 2005-11-16 (plot labels, plot names, normalize bias)
# version 0.0.09, 2005-11-18 (alternative majorization for Poisson Likelihood)
# version 0.0.10, 2005-11-20 (missing data in Poisson likelihood)
# version 0.0.11, 2005-11-20 (smacofSym fixed to handle negative targets)
# version 0.0.12, 2006-02-01 (fixed distAltOne)
# version 0.0.13, 2006-02-06 (small bugfixes)
# version 0.1.00, 2006-02-06 (switched completely to Poisson likelihoods)
# version 0.1.01, 2006-02-06 (use elegant instead of alscal for squared distances)
# version 0.2.00, 2006-02-09 (added quasi-symmetry/quasi-independence for hollow matrices)
# version 0.2.01, 2006-02-09 (use elegantRect for off-diagonal case)
# version 0.3.00, 2006-02-09 (added linear restrictions in the symmetric case)
# version 0.4.00, 2006-02-18 (added distAssColPart to deal with column partitioning)
# version 0.4.01, 2006-02-18 (small speedup by moving calculations out of loop)
# version 0.4.02, 2006-02-18 (added voronoi plots to distAssColPart)
# version 0.4.03, 2006-02-18 (nasty bug with indexing fixed in distAssColPart)
# version 0.4.04, 2006-02-20 (various speedsups, separate options for row and column effects)
# version 0.4.05, 2006-02-20 (share some common IO routines)
# version 0.4.06, 2006=03-07 (added license information)
#

require("smacof")
require("elegant")

distAssOne<-function(
	freq,					# data (single square matrix of frequencies)
	ndim=2,                 # dimensionality (default 2)
	outeps=1e-4,
	ineps=1e-4,
	reg=1e-6,
	rule="dist2",
	z="none",
	bias=c(TRUE,TRUE),
	outmax=500,
	inmax=1,
	theplot=TRUE) {
name<-deparse(substitute(freq))
nrow<-dim(freq)[1]; ncol<-dim(freq)[2]; itel<-1; freq<-as.matrix(freq)
mis<-ifelse(is.na(freq),0,1); dmis<-ifelse(is.na(freq),0,freq)
b<-rep(1,ncol); a<-rep(1,nrow)
rd<-rowSums(dmis); cd<-colSums(dmis)
e<-log(ifelse(dmis==0,1,dmis)); dmin<-2*(sum(dmis)-sum(dmis*e))
e<--.5*(e-outer(rowSums(e)/ncol,colSums(e)/nrow,"+")+(sum(e)/(ncol*nrow)))
z<-svd(e,nu=ndim,nv=0); x<-z$u;  x<-crossprod(e,x)
if (rule=="dist") phi<--distSym(x,reg) else phi<--dist2Sym(x)
ep<-exp(phi); lb<-ep; pdev<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
repeat {
	if (bias[1]) a<-rd/as.vector(((mis*ep)%*%b))
	if (bias[2]) b<-cd/as.vector((a%*%(mis*ep)))
	ab<-outer(a,b); lb<-ab*ep
	devab<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin	
	targ<--((phi)+((dmis-lb)/ab)); ww<-mis*ab; wt<-ww*targ; www<-ww+t(ww)
	targ<-(wt+t(wt))/ifelse(www==0,1,www)
	if (rule=="dist") 
		uu<-smacofSym(targ,w=www,init=x,p=ndim,itmax=inmax,reg=reg,verbose=FALSE)
    else uu<-elegantSym(targ,w=www,init=x,p=ndim,itmax=inmax,verbose=FALSE)
	x<-uu$x
	if (rule=="dist") phi<--distSym(x,reg) else phi<--dist2Sym(x)
	ep<-exp(phi); lb<-ab*ep; devxy<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
	writeIter(itel,pdev,devab,devxy) 
	if ((itel == outmax) || ((pdev-devxy)<outeps)) break()
	itel<-itel+1; pdev<-devxy
	}
if (theplot) myOnePlot(x,"symplot",rule,name,bias,"objects",labels(freq)[[2]])
return(list(x=x,a=a,b=b,dev=devxy))
}

distAssTwo<-function(
	freq,					# freq (single matrix of frequencies)
	ndim=2,                 # dimensionality (default 2)
	outeps=1e-4,
	ineps=1e-4,
	reg=1e-6,
	rule="dist2",
	bias=c(TRUE,TRUE),
	outmax=500,
	inmax=1,
	rowplot=FALSE,
	colplot=FALSE,
	jointplot=TRUE) {
name<-deparse(substitute(freq))
nrow<-dim(freq)[1]; ncol<-dim(freq)[2]; itel<-1; freq<-as.matrix(freq)
mis<-ifelse(is.na(freq),0,1); dmis<-ifelse(is.na(freq),0,freq)
b<-rep(1,ncol); a<-rep(1,nrow)
rd<-rowSums(dmis); cd<-colSums(dmis)
e<-log(ifelse(dmis==0,1,dmis)); dmin<-2*(sum(dmis)-sum(dmis*e))
e<--.5*(e-outer(rowSums(e)/ncol,colSums(e)/nrow,"+")+(sum(e)/(ncol*nrow)))
z<-svd(e,nu=ndim,nv=0); x<-z$u;  y<-crossprod(e,x)
if (rule=="dist") phi<--distRect(x,y,reg) else phi<--dist2Rect(x,y)
ep<-exp(phi); lb<-ep; pdev<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
repeat {
	if (bias[1]) a<-rd/as.vector(((mis*ep)%*%b))
	if (bias[2]) b<-cd/as.vector((a%*%(mis*ep)))
	ab<-outer(a,b); lb<-ab*ep
	devab<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin	
	targ<-(phi)+((dmis-lb)/ab)
	if (rule=="dist") 
		uu<-smacofRect(-targ,w=mis*ab,init=list(x,y),itmax=inmax,reg=reg,verbose=FALSE)
	else uu<-elegantRect(-targ,w=mis*ab,init=list(x,y),itmax=inmax,verbose=FALSE)
	x<-uu$x; y<-uu$y
	if (rule=="dist") phi<--distRect(x,y,reg) else phi<--dist2Rect(x,y)
	ep<-exp(phi); lb<-ab*ep
	devxy<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
	writeIter(itel,pdev,devab,devxy) 
	if ((itel == outmax) || ((pdev-devxy)<outeps)) break()
	itel<-itel+1; pdev<-devxy
	}
if (rowplot) myOnePlot(x,"rowplot",name,bias,"Row Objects",labels(freq)[[1]])
if (colplot) myOnePlot(y,"colplot",name,bias,"Column Objects",labels(freq)[[2]])
if (jointplot) myTwoPlot(x,y,"jointplot",rule,name,bias,
	"Column Objects (Red) Row Objects (Green)",labels(freq)[[1]],labels(freq)[[2]]) 
return(list(x=x,y=y,a=a,b=b,dev=devxy))
}

distAssLinOne<-function(
	freq,					# data (single square matrix of frequencies)
	ndim=2,                 # dimensionality (default 2)
	outeps=1e-4,
	ineps=1e-4,
	z="none",
	reg=1e-6,
	rule="dist2",
	bias=c(TRUE,TRUE),
	outmax=500,
	inmax=1,
	theplot=TRUE) {
name<-deparse(substitute(freq))
nrow<-dim(freq)[1]; ncol<-dim(freq)[2]; itel<-1; freq<-as.matrix(freq)
if (is.character(z)) z<-diag(n); zz<-qr(z); z<-qr.Q(zz)
mis<-ifelse(is.na(freq),0,1); dmis<-ifelse(is.na(freq),0,freq)
b<-rep(1,ncol); a<-rep(1,nrow)
rd<-rowSums(dmis); cd<-colSums(dmis)
e<-log(ifelse(dmis==0,1,dmis)); dmin<-2*(sum(dmis)-sum(dmis*e))
e<--.5*(e-outer(rowSums(e)/ncol,colSums(e)/nrow,"+")+(sum(e)/(ncol*nrow)))
if ((rule=="dist") || (rule == "dist2")) e<--e/2
ze<-svd(e,nu=ndim,nv=0); x<-ze$u;  x<-crossprod(e,x)
a<-crossprod(z,x); x<-z%*%a; d<-dist2Sym(x); r<-delta-d; rr<-rowSums(w*r)
if (rule=="dist") phi<--distSym(x,reg) else phi<--dist2Sym(x)
ep<-exp(phi); lb<-ep; pdev<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
repeat {
	if (bias[1]) a<-rd/as.vector(((mis*ep)%*%b))
	if (bias[2]) b<-cd/as.vector((a%*%(mis*ep)))
	ab<-outer(a,b); lb<-ab*ep
	devab<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin	
	targ<--((phi)+((dmis-lb)/ab)); ww<-mis*ab; wt<-ww*targ; www<-ww+t(ww)
	targ<-(wt+t(wt))/ifelse(www==0,1,www)
	if (rule=="dist") 
		uu<-smacofLinSym(targ,w=www,init=x,z=z,p=ndim,itmax=inmax,reg=reg,verbose=FALSE)
	else uu<-elegantLinSym(targ,w=www,init=x,z=z,p=ndim,itmax=inmax,verbose=FALSE)
	x<-uu$x
	if (rule=="dist") phi<--distSym(x,reg) else phi<--dist2Sym(x)
	ep<-exp(phi); lb<-ab*ep; devxy<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
	writeIter(itel,pdev,devab,devxy) 
	if ((itel == outmax) || ((pdev-devxy)<outeps)) break()
	itel<-itel+1; pdev<-devxy
	}
if (theplot) myOnePlot<-function(x,"symplot",name,bias,"objects",labels(freq)[[2]])
return(list(x=x,a=a,b=b,dev=devxy))
}

distAssColPart<-function(
	freq,					# freq (single matrix of frequencies)
	parts=dim(freq)[2],     # column partition
	ndim=2,                 # dimensionality 
	outeps=1e-4,
	ineps=1e-4,
	reg=1e-6,
	rule="dist2",
	bias=TRUE,
	outmax=500,
	inmax=1,
	rowplot=FALSE,
	colplot=FALSE,
	voronoi=TRUE,
	jointplot=FALSE) {
name<-deparse(substitute(freq)); nvar<-length(parts)
up<-cumsum(parts); lw<-1+c(0,up[-nvar])
nrow<-dim(freq)[1]; ncol<-dim(freq)[2]; itel<-1; freq<-as.matrix(freq)
rd<-matrix(0,nrow,nvar); cd<-rep(0,ncol)
mis<-ifelse(is.na(freq),0,1); dmis<-ifelse(is.na(freq),0,freq)
for (i in 1:nvar) {
	idi<-lw[i]:up[i]; dmi<-dmis[,idi]
	rd[,i]<-rowSums(dmi); cd[idi]<-colSums(dmi)
	}
b<-rep(1,ncol); a<-matrix(1,nrow,nvar)
e<-log(ifelse(dmis==0,1,dmis)); dmin<-2*(sum(dmis)-sum(dmis*e))
e<--.5*(e-outer(rowSums(e)/ncol,colSums(e)/nrow,"+")+(sum(e)/(ncol*nrow)))
if ((rule=="dist") || (rule == "dist2")) e<--e/2
z<-svd(e,nu=ndim,nv=0); x<-z$u;  y<-crossprod(e,x)
if (rule=="dist") phi<--distRect(x,y,reg) else phi<--dist2Rect(x,y)
ep<-exp(phi); lb<-ep; pdev<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
repeat {
    ab<-matrix(0,nrow,ncol)
	for (i in 1:nvar) {
		idi<-lw[i]:up[i]; emi<-mis[,idi]*ep[,idi]
		a[,i]<-rd[,i]/as.vector(emi%*%b[idi])
		if (bias) b[idi]<-cd[idi]/as.vector(a[,i]%*%emi)
		ab[,idi]<-outer(a[,i],b[idi])
		}
	lb<-ab*ep; devab<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin	
	targ<-(phi)+((dmis-lb)/ab)
	if (rule=="dist") 
		uu<-smacofRect(-targ,w=mis*ab,init=list(x,y),itmax=inmax,reg=reg,verbose=FALSE)
	else uu<-elegantRect(-targ,w=mis*ab,init=list(x,y),itmax=inmax,verbose=FALSE)
	x<-uu$x; y<-uu$y
	if (rule=="dist") phi<--distRect(x,y,reg) else phi<--dist2Rect(x,y)
	ep<-exp(phi); lb<-ab*ep
	devxy<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
	writeIter(itel,pdev,devab,devxy) 
	if ((itel == outmax) || ((pdev-devxy)<outeps)) break()
	itel<-itel+1; pdev<-devxy
	}
if (rowplot) myOnePlot(x,"rowplot",name,bias,"Row Objects",labels(freq)[[1]])
if (colplot) myOnePlot(y,"colplot",name,bias,"Column Objects",labels(freq)[[2]])
if (jointplot) myTwoPlot(x,y,"jointplot",rule,name,bias,
	"Column Objects (Red) Row Objects (Green)",labels(freq)[[1]],labels(freq)[[2]]) 
if (voronoi) {
	pdf(file=paste("voronoi_",rule,"_",name,".pdf",sep=""),encoding="MacRoman")
	mm1<-max(rbind(x,y)); mm2<-min(rbind(x,y)); mm<-c(mm2,mm1); lb<-labels(freq)[[2]]
	for (j in 1:nvar) {
		idi<-lw[j]:up[j]; yj<-y[idi,]
		plot(deldir(yj[,1],yj[,2],rw=c(mm,mm)), wlines = "tess", wpoints = "none", 
		col = c(1, 2, 1), lty = c(1, 1),
		xlab="dimension 1",ylab="dimension 2")
		title(main=paste("Voronoi",name,as.character(j)))
		text(yj,lb[idi],col="red")
		}
	dev.off()
	}
return(list(x=x,y=y,a=a,b=b,dev=devxy))
}

distAssLinTwo<-function(
	freq,					# freq (single matrix of frequencies)
	ndim=2,                 # dimensionality (default 2)
	outeps=1e-4,
	ineps=1e-4,
	z="none",
	reg=1e-6,
	rule="dist2",
	bias=TRUE,
	outmax=500,
	inmax=1,
	rowplot=FALSE,
	colplot=FALSE,
	jointplot=TRUE) {
name<-deparse(substitute(freq))
nrow<-dim(freq)[1]; ncol<-dim(freq)[2]; itel<-1; freq<-as.matrix(freq)
if (is.character(z)) z<-diag(n); zz<-qr(z); z<-qr.Q(zz)
mis<-ifelse(is.na(freq),0,1); dmis<-ifelse(is.na(freq),0,freq)
b<-rep(1,ncol); a<-rep(1,nrow)
rd<-rowSums(dmis); cd<-colSums(dmis)
e<-log(ifelse(dmis==0,1,dmis)); dmin<-2*(sum(dmis)-sum(dmis*e))
e<--.5*(e-outer(rowSums(e)/ncol,colSums(e)/nrow,"+")+(sum(e)/(ncol*nrow)))
if ((rule=="dist") || (rule == "dist2")) e<--e/2
z<-svd(e,nu=ndim,nv=0); x<-z$u;  y<-crossprod(e,x)
if (rule=="dist") phi<--distRect(x,y,reg) else phi<--dist2Rect(x,y)
ep<-exp(phi); lb<-ep; pdev<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
repeat {
	if (bias[1]) a<-rd/as.vector(((mis*ep)%*%b))
	if (bias[2]) b<-cd/as.vector((a%*%(mis*ep)))
	ab<-outer(a,b); lb<-ab*ep
	devab<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin	
	targ<-(phi)+((dmis-lb)/ab)
	if (rule=="dist") 
		uu<-smacofRect(-targ,w=mis*ab,init=list(x,y),itmax=inmax,reg=reg,verbose=FALSE)
    else uu<-elegantRect(-targ,w=mis*ab,init=list(x,y),itmax=inmax,verbose=FALSE)
	x<-uu$x; y<-uu$y
	if (rule=="dist") phi<--distRect(x,y,reg) else phi<--dist2Rect(x,y)
	ep<-exp(phi); lb<-ab*ep
	devxy<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
	writeIter(itel,pdev,devab,devxy) 
	if ((itel == outmax) || ((pdev-devxy)<outeps)) break()
	itel<-itel+1; pdev<-devxy
	}
if (rowplot) myOnePlot(x,"rowplot",name,bias,"Row Objects",labels(freq)[[1]])
if (colplot) myOnePlot(y,"colplot",name,bias,"Column Objects",labels(freq)[[2]])
if (jointplot) myTwoPlot(x,y,"jointplot",rule,name,bias,
	"Column Objects (Red) Row Objects (Green)",labels(freq)[[1]],labels(freq)[[2]]) 
return(list(x=x,y=y,a=a,b=b,dev=devxy))
}

quasiSym<-function(freq,itmax=100,eps=1e-6,sim=TRUE,verbose=TRUE) {
n<-dim(freq)[1]; diag(freq)<-0; itel<-1
lbd<-matrix(1,n,n); diag(lbd)<-0; pdev<--Inf
repeat {
	lbd<-lbd*(rowSums(freq)/rowSums(lbd))
	lbd<-t(t(lbd)*(colSums(freq)/colSums(lbd)))
	if (sim) lbd<-lbd*((freq+t(freq))/(lbd+t(lbd)+diag(n)))
	dev<-sum(freq*log(lbd+diag(n)))
	cat("Iteration: ",formatC(itel,digits=6,width=6),
	    "Previous Deviance: ",formatC(pdev,digits=6,width=12,format="f"),
		"Current Deviance: ",formatC(dev,digits=6,width=12,format="f"),
	    "\n")
	if ((itel == itmax) || ((dev-pdev)<eps)) break()
	itel<-itel+1; pdev<-dev
	}
return(list(lbd=lbd,dev=dev))
}

myOnePlot<-function(x,head,rule,name,bias,titl,labs) {
tb<-paste(ifelse(b[1],"R","r"),ifelse(b[2],"C","c"),sep="")
pdf(file=paste(head,"_",rule,"_",name,"_",tb,".pdf",sep=""),encoding="MacRoman")
plot(x,type="n",xlab="dimension 1",ylab="dimension 2")
title(main=paste(titl,name))
text(x,as.character(labs))
dev.off()
}

myTwoPlot<-function(x,y,head,rule,name,bias,titl,labs1,labs2) 
tb<-paste(ifelse(b[1],"R","r"),ifelse(b[2],"C","c"),sep="")
pdf(file=paste(head,"_",rule,"_",name,"_",tb,".pdf",sep=""),encoding="MacRoman")
plot(rbind(x,y),type="n",xlab="dimension 1",ylab="dimension 2")
title(main=paste("title",name))
text(x,as.character(labs1),col="RED")
text(y,as.character(labs2),col="GREEN")
dev.off()
}

writeIter<-function(itel,pdev,devab,devxy) {
cat("Iteration: ",formatC(itel,digits=6,width=6),
	"Previous Deviance: ",formatC(pdev,digits=6,width=12,format="f"),
	"After a and/or b: ",formatC(devab,digits=6,width=12,format="f"),
	"After X and/or Y: ",formatC(devxy,digits=6,width=12,format="f"),
	"\n")
}