#
#   lowRank 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-14,		first release
# version 0.2.00, 2006-04-15,		added Huber, Talwar, Barya methods
# version 0.3.00, 2006-04-16,		added symmetric case
#


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

ipRectLAV<-function(target,w=matrix(1,dim(target)[1],dim(target)[2]),ndim=2,init=FALSE,reg=1e-6,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*abs(target-xy)); itel<-1
repeat {
	res<-target-xy; ww<-w/sqrt((res^2)+reg)
	z<-ipRectLS(target,ww,init=list(x,y),ndim=ndim,itmax=1,verbose=FALSE)
	x<-z$x; y<-z$y; xy<-tcrossprod(x,y); nloss<-sum(w*abs(target-xy))
    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))
}

ipRectRobust<-function(
	target,
	w=matrix(1,dim(target)[1],dim(target)[2]),
	ndim=2,
	init=FALSE,
	fnc=list(f=fHuber,g=gHuber),
	cut=Inf,
	par=1,
	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*(fnc$f)(target-xy,cut,par)); itel<-1
repeat {
	res<-target-xy; ntarg<-xy+(fnc$g)(res,cut,par)
	z<-ipRectLS(ntarg,w,init=list(x,y),ndim=ndim,itmax=1,verbose=FALSE)
	x<-z$x; y<-z$y; xy<-tcrossprod(x,y); nloss<-sum(w*(fnc$f)(target-xy,cut,par))
    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))
}

fHuber<-function(x,cut,par) {
ifelse(abs(x)<cut,(x*x)/2,cut*abs(x)-((cut*cut)/2))
}

gHuber<-function(x,cut,par) {
ifelse(abs(x)<cut,x,cut*sign(x))
}

fTalwar<-function(x,cut,par) {
ifelse(abs(x)<cut,(x*x)/2,(cut*cut)/2)
}

gTalwar<-function(x,cut,par) {
ifelse(abs(x)<cut,x,0)
}

fBarya<-function(x,cut,par) {
ifelse(abs(x)<cut,(x*x)/2,((par*x*x)+(1-par)*cut*cut)/2)
}

gBarya<-function(x,cut,par) {
ifelse(abs(x)<cut,x,par*x)
}

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