#
#   logithom package
#   Copyright (C) 2006  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 1.0.0, 2006-05-05,		split off from scalassoc
# version 1.1.0, 2006-05-06,        ad-hoc homals can now do princals
# version 1.2.0, 2006-05-07,        special purpose matrix approximation program
# version 1.3.0, 2006-05-07,        logitHom can now do single nominal
# version 1.4.0, 2006-05-12,        now has timer and outfile like homals
#

logitHom<-function(
	dframe,					# data frame
	ndim=2,                 # dimensionality 
	outeps=1e-4,
	inseps=1e-4,
	inbeps=1e-4,
	bias=TRUE,
	rk=rep(ndim,ncol(dframe)),
	outmax=500,
	insmax=1,
	inbmax=1,
	offset=1.20,
	rowplot=FALSE,
	colplot=FALSE,
	voroplot=FALSE,
	vorolab=FALSE,
	movieplot=FALSE,
	movielab=FALSE,
	moviequick=FALSE,
	ffmpeg=moviequick,
	jointplot=FALSE,
	verbose=TRUE,
	outwrite=TRUE,
	objwrite=FALSE,
	timer=FALSE) {
if (timer) stime<-proc.time()
name<-deparse(substitute(dframe)); nvar<-dim(dframe)[2]
parts<-sapply(1:nvar,function(j) length(table(dframe[,j])))
mis<-ifelse(is.na(dframe),0,1); indi<-as.matrix(expandFrame(dframe))
up<-cumsum(parts); lw<-1+c(0,up[-nvar])
nrow<-dim(indi)[1]; ncol<-dim(indi)[2]; itel<-1
rd<-matrix(0,nrow,nvar); cd<-rep(0,ncol)
b<-rep(1,ncol); a<-mis; ab<-matrix(0,nrow,ncol)
for (j in 1:nvar) {
	idi<-lw[j]:up[j]; dmi<-indi[,idi]; cd[idi]<-colSums(dmi)
	ab[,idi]<-outer(a[,j],b[idi])
	}
dmin<-2*sum(indi)
z<-stripHomals(dframe,ndim,rk)
x<-z$x; y<-z$y
phi<-tcrossprod(x,y); ep<-exp(phi); lb<-ab*ep; pdev<-(2*(sum(indi*lb)-sum(xlogy(indi,lb))))-dmin
pterm<-sum(sapply(1:nvar,function(j) sum(mis[,j]*lb[,lw[j]:up[j]])))
pdev<-(2*(pterm-sum(xlogy(indi,lb))))-dmin
if (moviequick) system(paste("mkdir ",name,"_movie",sep=""))
if (outwrite)
	outfile<-file(paste(name,"out",sep="."),"w")
if (timer) itime<-proc.time()
if (outwrite)
	cat("Iterations:\n",file=outfile)
repeat {
	for (j in 1:nvar) {
		idi<-lw[j]:up[j]; emi<-ep[,idi]
		if (bias) b[idi]<-cd[idi]/as.vector(a[,j]%*%emi)
		a[,j]<-mis[,j]/as.vector(emi%*%b[idi])
		ab[,idi]<-outer(a[,j],b[idi])
		}	
	lb<-ab*ep
	pterm<-sum(sapply(1:nvar,function(j) sum(mis[,j]*lb[,lw[j]:up[j]])))
    devab<-(2*(pterm-sum(xlogy(indi,lb))))-dmin
    targ<-phi+2*(indi-lb)
	uu<-stripRect(targ,mis,ndim,x,y,rk,lw,up)
	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
	pterm<-sum(sapply(1:nvar,function(j) sum(mis[,j]*lb[,lw[j]:up[j]])))
	devxy<-(2*(pterm-sum(xlogy(indi,lb))))-dmin
	if (verbose)
		writeIter(itel,pdev,devab,devxy) 
	if (outwrite)
		writeIter(itel,pdev,devab,devxy,myfile=outfile) 		
	if ((itel == outmax) || ((pdev-devxy) < outeps)) break()
	itel<-itel+1; pdev<-devxy
	}
if (timer) ctime<-proc.time()
if (ffmpeg) {
	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(indi)[[1]])
if (colplot) myOnePlot(y,"colplot",name,bias,"Column Objects",labels(indi)[[2]])
if (jointplot) myTwoPlot(x,y,"jointplot",name,bias,
	"Column Objects (Red) Row Objects (Green)",labels(indi)[[1]],labels(indi)[[2]]) 
if (voroplot) {
	pdf(file=paste("voronoi_",name,".pdf",sep=""),encoding="MacRoman")
	mm1<-max(rbind(x,y)); mm2<-min(rbind(x,y)); mm<-c(mm2,mm1); lb<-labels(indi)[[2]]
	for (j in 1:nvar) {
		idi<-lw[j]:up[j]; yj<-y[idi,]; z<-rbind(x,yj); vname<-names(dframe)[j]
		xx<-c(min(z[,1]),max(z[,1])); yy<-c(min(z[,2]),max(z[,2]))
		plot(x,type="n",main=paste("Voronoi plot for",name,":",vname),
		xlab="dimension 1",ylab="dimension 2",
		xlim=offset*xx,ylim=offset*yy)
		drawEdges(yj,b[idi],far=1000)
		if (vorolab) text(x,as.character(dframe[,j]))
		text(yj,lb[idi],col="red")
		}
	dev.off()
	}

#if (outwrite) {
#	writeHead(name,vname[j],ndim,active[j],rank[j],level[j],outfile)
#	writeY(data[,j],y,"Y",outfile)
#	writeY(data[,j],c,"C",outfile)
#	writeY(data[,j],u,"Z",outfile)
#	writeA(t(a),outfile)
#	}

if (objwrite) {
	sink(paste(name,"obj",sep="."))
	rownames(x)<-labels(dframe)[[1]]
	print(x)
	sink()
	}
if (timer && outwrite) {
	otime<-proc.time()
	cat("Input time:     ",formatC(itime-stime,digits=6,width=10,format="f"),"\n",file=outfile)
	cat("Compute time:   ",formatC(ctime-itime,digits=6,width=10,format="f"),"\n",file=outfile)
	cat("Iteration time: ",formatC((ctime-itime)/itel,digits=6,width=10,format="f"),"\n",file=outfile)
	cat("Output time:    ",formatC(otime-ctime,digits=6,width=10,format="f"),"\n",file=outfile)
	}
if (outwrite) close(outfile)
return(list(x=x,y=y,a=a,b=b,dev=devxy))
}


expandFrame<-function(tab,clean=TRUE,zero=TRUE,returnFrame=TRUE) {
n<-dim(tab)[1]; m<-dim(tab)[2]; g<-matrix(0,n,0); l<-rep("",0)
lab1<-labels(tab)[[1]]; lab2<-labels(tab)[[2]]
for (j in 1:m) {
	y<-as.factor(tab[,j]); h<-levels(y)
	g<-cbind(g,ifelse(outer(y,h,"=="),1,0))
	l<-c(l,paste(lab2[j],"_",h,sep=""))
	}
if (zero) g<-ifelse(is.na(g),0,as.matrix(g))
if (clean) {
	g<-g[which(rowSums(g)>0),which(colSums(g)>0)]
	g<-g[,which(colSums(g)<n)]
	}
if (!returnFrame) return(g)
g<-as.data.frame(g,row.names=lab1)
names(g)<-l
return(g)
}

stripRect<-function(target,w,ndim=2,x,y,rk,lw,up,itmax=100,eps=1e-6,verbose=FALSE) {
nobj<-nrow(target); nvar<-length(rk); xy<-tcrossprod(x,y); oloss<-0; itel<-1
for (j in 1:nvar) {
	indi<-lw[j]:up[j]
	oloss<-oloss+sum(w[,j]*(target[,indi]-xy[,indi])^2)
	}
repeat {
    sj<-matrix(0,nobj,ndim)
	for (j in 1:nvar) {
		indi<-lw[j]:up[j]; yj<-y[indi,]
		tj<-target[,indi]; sj<-sj+w[,j]*tj%*%yj
		}
	for (i in 1:nobj) {
		cj<-matrix(0,ndim,ndim)
	    for (j in 1:nvar) {
			indi<-lw[j]:up[j]; yj<-y[indi,]
			cj<-cj+w[i,j]*crossprod(yj)
			}
		x[i,]<-solve(cj,drop(sj[i,]))
		}
    xy<-tcrossprod(x,y); xloss<-0
	for (j in 1:nvar) {
		indi<-lw[j]:up[j]
		xloss<-xloss+sum(w[,j]*(target[,indi]-xy[,indi])^2)
		}
	for (j in 1:nvar) {
		indi<-lw[j]:up[j]; tj<-target[,indi];
		sj<-crossprod(x,w[,j]*tj); cj<-crossprod(x,w[,j]*x)
		yj<-t(solve(cj,sj))
		if (rk[j] < ndim) {
			rj<-chol(cj); zj<-lowRank(tcrossprod(yj,rj),rk[j])
			yj<-t(solve(rj,t(zj)))
			}
		y[indi,]<-yj
		}
    xy<-tcrossprod(x,y); yloss<-0
	for (j in 1:nvar) {
		indi<-lw[j]:up[j]
		yloss<-yloss+sum(w[,j]*(target[,indi]-xy[,indi])^2)
		}
    if (verbose)
	    cat("Iteration: ",formatC(itel,digits=6,width=6),
			"Previous Loss: ",formatC(oloss,digits=6,width=12,format="f"),
			"Loss after X: ",formatC(xloss,digits=6,width=12,format="f"),
			"Loss after Y: ",formatC(yloss,digits=6,width=12,format="f"),
			"\n")
    if (((oloss-yloss) < eps) || (itel == itmax)) break()
    oloss<-yloss; itel<-itel+1
    }
q<-qr(x); x<-qr.Q(q); y<-tcrossprod(y,qr.R(q))
return(list(x=x,y=y))
}

stripHomals<-function(dframe,ndim=2,rk=rep(ndim,ncol(dframe)),eps1=1e-6,eps2=1e-6,itmax=100,verbose=FALSE) {            
nobj<-dim(dframe)[1]; nvar<-dim(dframe)[2]; iter<-0; pops<-0; 
mis<-rowSums(ifelse(is.na(dframe),0,1))
for (j in 1:nvar) dframe[,j]<-as.factor(dframe[,j]) 
x<-cbind(orthogonalPolynomials(mis,1:nobj,ndim))
x<-normX(centerX(x,mis),mis)$q
repeat {
	iter<-iter+1; y<-matrix(0,0,ndim); tot<-matrix(0,nobj,ndim)
	for (j in 1:nvar) {
		gg<-dframe[,j]; ii<-which(!is.na(gg)); dd<-as.vector(table(gg))
		yj<-apply(x, 2, function(z) tapply(z,gg,mean))
		if (rk[j] < ndim) yj<-lowRank(sqrt(dd)*yj,rk[j])/sqrt(dd)
		y<-rbind(y,yj); tot[ii,]<-tot[ii,]+yj[gg[ii],]
		}	
	qv<-normX(centerX((1/mis)*tot,mis),mis)
	z<-qv$q;r<-qv$r;ops=sum(r); aps<-sum(La.svd(crossprod(x,mis*z),0,0)$d)/ndim	
    if (verbose)
	    cat("Iteration: ",formatC(iter,digits=6,width=6),
			"Previous Gain: ",formatC(pops,digits=6,width=12,format="f"),
			"Current  Gain: ",formatC(ops,digits=6,width=12,format="f"),
			"Current  Corr: ",formatC(aps,digits=6,width=12,format="f"),
			"\n")
	if (((ops - pops) < eps1) || ((1.0 - aps) < eps2) || (iter == itmax)) break()
	x<-z; pops<-ops; iter<-iter+1
	}
return(list(x=z,y=y))
}

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

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

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

writeIter<-function(itel,pdev,devab,devxy,myfile="") {
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",file=myfile)
}

writeHead<-function(name,vname,p,a,r,c,ofile) {
s<-NULL; sl<-35+10*p; for (i in 1:sl) s<-paste(s,"*",sep="")
cat("\n",formatC(s,format="s"),"\n",file=ofile,sep="")
cat(formatC(vname,format="s"),
if (a) "(Active)," else "(Passive),",
"Rank =", formatC(r,width=1),"and",
"Level is", 
switch(c,"NO"="nominal","OR"="ordinal","NU"="numerical","PO"="polynomial"),
"\n",file=ofile)
cat(formatC(s,format="s"),"\n",file=ofile)
}

writeY<-function(g,y,t,ofile) {
d<-as.vector(table(g));l<-levels(g)
s<-NULL; sl<-35+10*dim(y)[2]; for (i in 1:sl) s<-paste(s,"*",sep="")
u<-NULL; sl<-35+10*dim(y)[2]; for (i in 1:sl) u<-paste(u,"-",sep="")
switch(t,"Z"=cat("Lower Rank Quantifications\n",file=ofile),"C"=cat("Category Centroids\n",file=ofile),
"Y"=cat("Rank-restricted Category Quantifications\n",file=ofile))
cat(formatC(u,format="s"),"\n",file=ofile)
for (k in 1:length(d))
cat(formatC(l[k],width=10,format="s"),
formatC(d[k],width=5)," *** ",
formatC(y[k,],digits=6,width=10,format="f")," *** ",
formatC(sum(d[k]*y[k,]^2),digits=6,width=10,format="f"),"\n",sep="",file=ofile)
cat(formatC(u,format="s"),"\n",file=ofile)
cat(formatC(" ",format="s",width=10),
formatC(sum(d),width=5)," *** ",
formatC(d%*%(y^2),digits=6,width=10,format="f")," *** ",
formatC(sum(d%*%(y^2)),digits=6,width=10,format="f"),"\n",sep="",file=ofile)
cat(formatC(s,format="s"),"\n",file=ofile)
}

writeA<-function(a,ofile) {
s<-NULL; sl<-35+10*dim(a)[2]; for (i in 1:sl) s<-paste(s,"*",sep="")
u<-NULL; sl<-35+10*dim(a)[2]; for (i in 1:sl) u<-paste(u,"-",sep="")
cat("Category Loadings\n",file=ofile)
cat(formatC(u,format="s"),"\n",file=ofile)
for (k in 1:dim(a)[1])
cat(formatC(" ",format="s",width=10),formatC(k,width=5)," *** ",
formatC(a[k,],digits=6,width=10,format="f")," *** ",
formatC(sum(a[k,]^2),digits=6,width=10,format="f"),"\n",sep="",file=ofile)
cat(formatC(u,format="s"),"\n",file=ofile)
cat(formatC(" ",format="s",width=15)," *** ",
formatC(apply(a^2,2,sum),digits=6,width=10,format="f")," *** ",
formatC(sum(a^2),digits=6,width=10,format="f"),"\n",sep="",file=ofile)
cat(formatC(s,format="s"),"\n",file=ofile)
} 

weightedGramSchmidt<-function(x,w) {
ss<-NULL; 
for (j in 1:dim(x)[2]) {
if (j > 1) {xx<-x[,1:(j-1)]; x[,j]<-x[,j]-xx%*%(crossprod(xx,(w*x[,j])))}
s<-sqrt(sum(w*x[,j]^2)); ss<-c(ss,s); x[,j]<-x[,j]/s;}
list(pol=x,fac=ss)
}

orthogonalPolynomials<-function(w,x,p) {
	z<-weightedGramSchmidt(outer(x,0:p,"^"),w)$pol[,2:(p+1)]
	}

centerX<-function(x,w) apply(x,2,function(z) z-weighted.mean(z,w))

normX<-function(x,w) {
	qq<-La.svd(sqrt(w)*x); list(q=(1/sqrt(w))*(qq$u),r=qq$d)
	}

makeEdges<-function(a,b) {
c<--log(b); n<-length(b); lns<-matrix(0,0,8)
for (i in 1:(n-1)) {
	for(j in (i+1):n) {
	    dd<-a[i,]-a[j,]; dc<-c[i]-c[j]; ss<-sum(dd^2)
	    if (is.nul(ss)) next() 
	    ee<-dc*dd/ss; ff<-c(-dd[2],dd[1])
		xlw<--Inf; xup<-Inf
		for (k in (1:n)[-c(i,j)]) {
		    dd<-a[i,]-a[k,]; dc<-c[i]-c[k]
		    mum<-sum(dd*ff); mom<-dc-sum(dd*ee)
		    if (is.nul(mum) & (mom > 0)) {
		    	xlw<-Inf; xup<--Inf
		    	}
		    if (mum>0) xlw<-max(xlw,mom/mum)
		    if (mum<0) xup<-min(xup,mom/mum)
			}
		if (xlw<xup) lns<-rbind(lns,c(i,j,ee,ff,xlw,xup))
		}
	}
return(lns)
}

drawEdges<-function(a,b,far=1000) {
lns<-makeEdges(a,b)
p<-dim(lns)[1]; n<-dim(a)[1]
for (i in 1:p) {
	ee<-lns[i,3:4]; ff<-lns[i,5:6]
	xlw<-lns[i,7]; if (xlw == -Inf) xlw<--far
	plw<-ee+xlw*ff
	xup<-lns[i,8]; if (xup == Inf) xup<-far
	pup<-ee+xup*ff
	lines(rbind(plw,pup),col="RED")
	}
}

lowRank<-function(x,p) {
z<-svd(x)
return((z$u[,1:p])%*%((z$d[1:p])*t(z$v[,1:p])))
}

xlogy<-function(x,y) {
return(ifelse(((x==0) & (y==0)),0,x*log(y)))
}

is.nul<-function(x) {
return(abs(x)<1e-15)
}