#
#   multiway package
#   Copyright (C) 2008  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.0, 2008-03-27,		first release
# version 0.2.0, 2008-03-28,		added tucker and constraints
# version 0.3.0, 2008-04-01,        added non-negativity constraints to candecomp
#

require("quadprog")
require("nnls")

candecomp<-function(a,x,ortho=rep(FALSE,length(x)),ispos=rep(FALSE,length(x)),itmax=1000,eps=1e-6,verbose=FALSE) {
	ndam<-dim(a)
	nard<-length(ndam)
	ndim<-ncol(x[[1]])
	for (k in 1:nard) if (ortho[k]) x[[k]]<-procrustus(x[[k]])
	c<-lapply(x,crossprod)
	oloss<-candeValue(a,x)$loss
	itel<-1
	repeat {
		for (k in 1:nard) {
			cc<-arrHadamard(c[-k])
			for (p in 1:ndim) {
				y<-lapply(x[-k],function(z) z[,p])
				b<-arrOuter(y)
				x[[k]][,p]<-apply(a,k,function(z) sum(z*b))
				}
			xx<-x[[k]]
			if (ortho[k]) x[[k]]<-procrustus(xx)
				else {
					if (!ispos) x[[k]]<-t(solve(cc,t(xx)))
						else x[[k]]<-posProj(xx,cc)
					}
			c[[k]]<-crossprod(x[[k]])
			}
			cval<-candeValue(a,x)
			nloss<-cval$loss; ahat<-cval$ahat
		if (verbose)
			cat("Iteration: ",formatC(itel,digits=3,width=3),
				"Old Loss: ",formatC(oloss,digits=6,width=10,format="f"),
				"New Loss: ",formatC(nloss,digits=6,width=10,format="f"),
				"\n")
		if ((itel == itmax) | ((oloss - nloss) < eps)) break()
		itel<-itel+1; oloss<-nloss
		}
	return(list(x=x,ahat=ahat,loss=nloss))
	}

candeValue<-function(a,x) {
	ndim<-ncol(x[[1]])
	ahat<-array(0,dim(a))
	for (p in 1:ndim) {
		y<-sapply(x,function(z) z[,p])
		b<-arrOuter(y)
		ahat<-ahat+b
		}
	loss<-sum((a-ahat)^2)
	return(list(ahat=ahat,loss=loss))
	}
	
tucker<-function(a,x,ident=rep(FALSE,length(x)),ortho=rep(FALSE,length(x)),ispos=rep(FALSE,length(x)),isposb=FALSE,itmax=1000,eps=1e-6,verbose=FALSE) {
	ndam<-dim(a)
	for (k in 1:nard) if (ident[k]) x[[k]]<-diag(ndam[k])
	ndbm<-sapply(x,function(z) ncol(z))
	nard<-length(ndam)
	rard<-rev(1:nard)
	x<-lapply(x,procrustus)
	xx<-arrKronecker(x)
	aa<-as.vector(aperm(a,rard))
	if (!isposb) bb<-lsfit(xx,aa,intercept=FALSE)
		else bb<-nnls(xx,aa)
	b<-aperm(array(bb,rev(ndbm)),rard)
	ahat<-aperm(array(drop(xx%*%bb),rev(ndam)),rard)
	oloss<-sum((a-ahat)^2)
	itel<-1
	repeat {
		for (k in 1:nard) {	
			if (!ident[k]) {
				z<-crossprod(flatten(a,k),arrKronecker(x[-k]))%*%flatten(b,k)
				x[[k]]<-procrustus(z)
				}
			}
		xx<-arrKronecker(x)
		aa<-as.vector(aperm(a,rard))
		if (!isposb) bb<-lsfit(xx,aa,intercept=FALSE)
			else bb<-nnls(xx,aa)
		b<-aperm(array(bb,rev(ndbm)),rard)
		ahat<-aperm(array(drop(xx%*%bb),rev(ndam)),rev(1:nard))
		nloss<-sum((a-ahat)^2)
		if (verbose)
			cat("Iteration: ",formatC(itel,digits=3,width=3),
				"Old Loss: ",formatC(oloss,digits=6,width=10,format="f"),
				"New Loss: ",formatC(nloss,digits=6,width=10,format="f"),
				"\n")
		if ((itel == itmax) | ((oloss - nloss) < eps)) break()
		itel<-itel+1; oloss<-nloss
		}
	return(list(x=x,b=b,ahat=ahat,loss=nloss))
	}
		
arrHadamard<-function(c,fun=function(x,y) x*y) {
	nmat<-length(c)
	if (nmat == 0) stop("empty argument in arrHadamard")
	res<-c[[1]]
	if (nmat == 1) return(res)
	for (i in 2:nmat) res<-fun(res,c[[i]])
	return(res)
	}

arrOuter<-function(x,fun="*") {
	nmat<-length(x)
	if (nmat == 0) stop("empty argument in arrOuter")
	res<-x[[1]]
	if (length(x) == 1) return(res)
	for (i in 2:nmat) res<-outer(res,x[[i]],fun)
	return(res)
	}
	
arrKronecker<-function(x,fun="*") {
	nmat<-length(x)
	if (nmat == 0) stop("empty argument in arrKronecker")
	res<-x[[1]]
	if (length(x) == 1) return(res)
	for (i in 2:nmat) res<-kronecker(res,x[[i]],fun)
	return(res)
	}
	
flatten<-function(a,k) {
	nard<-rev(1:(length(dim(a))-1))
	apply(a,k,function(z) as.vector(aperm(z,nard)))
}

procrustus<-function(x) {
	res<-svd(x)
	return(tcrossprod(res$u,res$v))
}

posProj<-function(x,c) {
n<-nrow(x); p<-ncol(x); ip<-diag(p)
u<-matrix(0,n,p)
for (i in 1:n) {
	y<-x[i,]
	u[i,]<-solve.QP(c,y,ip)$solution
	}
return(u)
}