#
#   scalAssoc 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.1.00, 2006-04-13,		first release
# version 0.1.01, 2006-04-14,       bug fixes
# version 0.1.02, 2006-04-14,       added movieplot
# version 0.1.03, 2006-04-14,       added biplot scaling
# version 0.1.04, 2006-04-17,       labels in movieplot
# version 0.1.05, 2006-04-17,		labels in voronoi plot
# version 0.1.06, 2006-04-17,       creates quicktime movie
# version 0.2.00, 2006-04-19,		added scalAssocFreq
# version 0.3.00, 2006-04-20,		added scalAssocQSym
# version 0.3.01, 2006-04-20,		bug fixes
# version 0.3.02, 2006-04-20,		improve initial using ca()
# version 0.4.0,  2006-05-04,		new pseudo-voronoi
# version 0.4.1,  2006-05-04,		homals initial estimate
# version 1.0.0,  2006-05-05,		removed indicators to logithom
#
# to be added asap:  
# 	inner iterations for a and b
#

scalAssocFreq<-function(
	freq,					# cross table
	ndim=2,                 # dimensionality 
	outeps=1e-4,
	inseps=1e-4,
	inbeps=1e-4,
	bias=c(TRUE,TRUE),
	outmax=500,
	insmax=1,
	inbmax=1,
	rowplot=FALSE,
	colplot=FALSE,
	movieplot=FALSE,
	movielab=FALSE,
	moviequick=FALSE,
	jointplot=FALSE) {
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); ab<-outer(a,b)
rd<-rowSums(dmis); cd<-colSums(dmis); ww<-mis*rd
e<-log(ifelse(dmis==0,1,dmis)); dmin<-2*(sum(dmis)-sum(dmis*e))
z<-ca(freq,pdim=ndim,scale="YCX"); x<-z$x; y<-z$y; phi=tcrossprod(x,y)
ep<-exp(phi); lb<-ab*ep; pdev<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
if (moviequick) system(paste("mkdir ",name,"_movie",sep=""))
repeat {
	if (bias[2]) b<-cd/as.vector((a%*%(mis*ep)))
	if (bias[1]) a<-rd/as.vector(((mis*ep)%*%b))
	ab<-outer(a,b); lb<-ab*ep
	devab<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin	
    targ<-phi+2*(dmis-lb)/rd
	uu<-ipRect(targ,ww,ndim,init=list(x,y),itmax=insmax,eps=inseps,verbose=FALSE)
	if (movieplot) myMoviePlot(x,uu$x,movielab)
	if (moviequick) {
	    setwd(paste(name,"_movie",sep=""))
		jpeg(paste("file",as.character(itel),".jpg",sep=""))
		myMoviePlot(x,uu$x,movielab)
		dev.off()
		setwd("..")
		}
	x<-uu$x; y<-uu$y; phi<-tcrossprod(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 (moviequick) {
	setwd(paste(name,"_movie",sep=""))
	system(paste("ffmpeg -i file%d.jpg ",name,".mov",sep=""))
	system("mv *.mov ..")
	system("rm *.jpg")
	setwd("..")
	system(paste("rm -rf ",name,"_movie",sep=""))
	}
ez<-eigen(crossprod(y)); ezval<-sqrt(sqrt(ez$values)); ezvec<-ez$vectors
x<-sqrt(sqrt(nrow/ncol))*x%*%ezvec%*%diag(ezval); y<-sqrt(sqrt(ncol/nrow))*y%*%ezvec%*%diag(1/ezval)
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",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))
}

scalAssocQSym<-function(
	freq,					# cross table
	ndim=2,                 # dimensionality 
	outeps=1e-4,
	inseps=1e-4,
	inbeps=1e-4,
	bias=c(TRUE,TRUE),
	outmax=500,
	insmax=1,
	inbmax=1,
	theplot=FALSE,
	movieplot=FALSE,
	movielab=FALSE,
	moviequick=FALSE) {
require("ca")
name<-deparse(substitute(freq))
nrow<-ncol<-dim(freq)[1]; 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); ab<-outer(a,b)
rd<-rowSums(dmis); cd<-colSums(dmis); ww<-mis*rd
e<-log(ifelse(dmis==0,1,dmis)); dmin<-2*(sum(dmis)-sum(dmis*e))
z<-ca(freq,pdim=ndim,scale="GOO"); x<-z$x; phi<-tcrossprod(x)
ep<-exp(phi); lb<-ab*ep; pdev<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin
if (moviequick) system(paste("mkdir ",name,"_movie",sep=""))
repeat {
	if (bias[2]) b<-cd/as.vector((a%*%(mis*ep)))
	if (bias[1]) a<-rd/as.vector(((mis*ep)%*%b))
	ab<-outer(a,b); lb<-ab*ep
	devab<-(2*(sum(mis*lb)-sum(dmis*log(lb))))-dmin	
    targ<-phi+2*(dmis-lb)/rd
	uu<-ipSymLS(targ,ww,ndim,init=x,itmax=insmax,eps=inseps,verbose=FALSE)
	if (movieplot) myMoviePlot(x,uu$x,movielab)
	if (moviequick) {
	    setwd(paste(name,"_movie",sep=""))
		jpeg(paste("file",as.character(itel),".jpg",sep=""))
		myMoviePlot(x,uu$x,movielab)
		dev.off()
		setwd("..")
		}
	x<-uu$x; phi<-tcrossprod(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 (moviequick) {
	setwd(paste(name,"_movie",sep=""))
	system(paste("ffmpeg -i file%d.jpg ",name,".mov",sep=""))
	system("mv *.mov ..")
	system("rm *.jpg")
	setwd("..")
	system(paste("rm -rf ",name,"_movie",sep=""))
	}
if (theplot) myOnePlot(x,"rowplot",name,bias,"Objects",labels(freq)[[1]])
return(list(x=x,a=a,b=b,dev=devxy))
}

ipRect<-function(target,w=matrix(1,dim(target)[1],dim(target)[2]),ndim=2,init=FALSE,itmax=100,eps=1e-6,verbose=FALSE) {
if (is.list(init)) {
	x<-init[[1]]; y<-init[[2]]
	}
else {
	z<-svd(target); x<-z$u[,1:ndim]; y<-crossprod(target,x)
	}
n<-dim(target)[1]; m<-dim(target)[2]; xy<-tcrossprod(x,y)
oloss<-sum(w*(target-xy)^2); itel<-1
repeat {
	for (i in 1:n) {
		wi<-w[i,]; wti<-wi*target[i,]
		x[i,]<-solve(crossprod(y,wi*y),as.vector(wti%*%y))
		}
	for (j in 1:m) {
		wj<-w[,j]; wtj<-wj*target[,j]
		y[j,]<-solve(crossprod(x,wj*x),as.vector(wtj%*%x))
		}
    xy<-tcrossprod(x,y); nloss<-sum(w*(target-xy)^2)
    if (verbose)
	    cat("Iteration: ",formatC(itel,digits=6,width=6),
			"Previous Loss: ",formatC(oloss,digits=6,width=12,format="f"),
			"Current Loss: ",formatC(nloss,digits=6,width=12,format="f"),
			"\n")
    if (((oloss-nloss) < eps) || (itel == itmax)) break()
    oloss<-nloss; itel<-itel+1
    }
q<-qr(x); x<-qr.Q(q); y<-tcrossprod(y,qr.R(q))
return(list(x=x,y=y))
}

ipSymLS<-function(target,w=matrix(1,dim(target)[1],dim(target)[2]),ndim=2,init=FALSE,itmax=100,eps=1e-6,verbose=FALSE) {
if (is.matrix(init)) x<-init
else {
	z<-eigen(target); x<-z$vectors[,1:ndim]; x<-x%*%diag(sqrt(z$values[1:ndim]))
	}
n<-dim(target)[1]; xx<-tcrossprod(x); oloss<-sum(w*(target-xx)^2); itel<-1
repeat {
    for (i in 1:n) {
    	ai<-crossprod(x,w[i,]*x)-w[i,i]*outer(x[i,],x[i,])
    	bi<-colSums((w[i,]*target[i,])*x)-w[i,i]*target[i,i]*x[i,]
    	if (w[i,i] == 0) x[i,]<-solve(ai,bi)
    	else {
    	    li<-sum(x[i,]^2); zi<-x[i,]/sqrt(li); wi<-w[i,i]; ci<-target[i,i]
    		a4<-wi; a3<-0; a2<-2*(sum(zi*(ai%*%zi))-(wi*ci)); a1<--4*sum(zi*bi) 
    		b1<-a1/(2*a4); b2<-a2/(3*wi); bb<-as.complex((b1^4)+2*(b1^2)*(b2^3))
    		s<-(b1^2)+(b2^3)+sqrt(bb); kf<-s^(1/3)+(b2^2)*s^(-1/3)+b2; li<-Re(-b1/kf)
    		x[i,]<-li*zi; aa<-2*(li^2)*ai; bb<-2*li*bi; ei<-eigen(aa)
    		l<-ei$values; k<-ei$vector; kb<-drop(crossprod(k,bb)); ml<-min(l)
    		ul<-ml-sqrt(sum(kb^2)); uu<-ml-sqrt(sum(kb[which(l==ml)]^2)); u<-uu
    		fu<-function(u) sum((kb/(l-u))^2)-1
    		gu<-function(u) 2*sum((kb^2)/((l-u)^3))
    		repeat {
    			if (abs(fu(u))<1e-6) break()
    			u<-u-fu(u)/gu(u)
    			}
    		x[i,]<-li*(k%*%(kb/(l-u)))
     		}
     	}
    xx<-tcrossprod(x); nloss<-sum(w*(target-xx)^2)
    if (verbose)
	    cat("Iteration: ",formatC(itel,digits=6,width=6),
			"Previous Loss: ",formatC(oloss,digits=6,width=12,format="f"),
			"Current Loss: ",formatC(nloss,digits=6,width=12,format="f"),
			"\n")
    if (((oloss-nloss) < eps) || (itel == itmax)) break()
    oloss<-nloss; itel<-itel+1
    }
return(list(x=x))
}

iterBias<-function() {}

myOnePlot<-function(x,head,name,bias,titl,labs) {
if (bias) tb<-"B" else tb<-"N"
pdf(file=paste(head,"_",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,name,bias,titl,labs1,labs2) {
if (bias) tb<-"B" else tb<-"N"
pdf(file=paste(head,"_",name,"_",tb,".pdf",sep=""),encoding="MacRoman")
plot(rbind(x,y),type="n",xlab="dimension 1",ylab="dimension 2")
title(main=paste(titl,name))
text(x,as.character(labs1),col="RED")
text(y,as.character(labs2),col="GREEN")
dev.off()
}

myMoviePlot<-function(x,y,movielab=FALSE) {
plot(x,xlab="dimension 1",ylab="dimension 2",type="n")
if (movielab) {
	text(x,as.character(1:dim(x)[1]),pch="+",col="GREEN")
	text(y,as.character(1:dim(x)[1]),pch="+",col="RED")
	}
else {
	points(x,pch="+",col="GREEN")
	points(y,pch="+",col="RED")
	}
sapply(1:dim(x)[1],function(i) lines(matrix(c(x[i,1:2],y[i,1:2]),2,2,byrow=TRUE)))
}

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")
}