#
#   lssim 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 0.1.00, 2006-10-22,		first release
# version 0.1.01, 2006-10-22,		bug fixes
# version 0.1.02, 2006-10-22,		iteration output improvements
# version 0.2.00, 2006-10-22,		added reduced rank regression
# version 0.3.00, 2006-10-22,		lsfa can now also do pca (i.e. ignore errors)
# version 0.4.00, 2006-10-22,		lsrrr can now ignore errors
# version 0.4.01, 2006-10-25,		minor efficiency improvements
# version 0.4.02, 2006-10-25,		lsrrr can now have nonorthogonal errors
#
# todo
#
# -- put in optimal scaling
# -- alternative algorithm
# -- extend to other sem
# -- put in general block algorithms
#

lsfa<-function(x,p,itmax=100,eps=1e-6,ignore=FALSE) {
	n<-nrow(x); m<-ncol(x)
	u<-ranort(n,p); if (!ignore) e<-ranort(n,m); sold<-Inf; itel<-1
	repeat {
		a<-crossprod(u,x); if (!ignore) d<-colSums(e*x)
		y<-u%*%a; if (!ignore) y<-y+t(t(e)*d)
		sa<-sum((x-y)^2)
		w<-tcrossprod(x,a); if (!ignore) w<-cbind(w,t(t(x)*d))
		v<-procrus(w); u<-v[,1:p]; if (!ignore) e<-v[,p+1:m]
		y<-u%*%a; if (!ignore) y<-y+t(t(e)*d)
		sb<-sum((x-y)^2)
		cat("Iteration: ",formatC(itel,digits=6,width=6),
			"   Loss Start:   ",formatC(sold,digits=6,width=12,format="f"),
			"   Loss Loads:   ",formatC(sa,digits=6,width=12,format="f"),
			"   Loss Scors: ",formatC(sb,digits=6,width=12,format="f"),"\n")
		if (((sold-sb)<eps) || (itel==itmax)) break()
		itel<-itel+1; sold<-sb
	}
if (!ignore) return(list(s=sb,u=u,a=a,e=e,d=d)) 
	else return(list(s=sb,u=u,a=a))
}

lsrrr<-function(x,y,p,eps=1e-6,itmax=100,de=TRUE,ignore=FALSE){
	n<-nrow(x); m<-ncol(x); k<-ncol(y)
	if (!ignore) e<-ortcrus(ranort(n,k),x); z<-matrix(0,m,k); sold<-Inf; itel<-1
	repeat {
		if (!ignore} {
			d<-colSums(e*(y-x%*%z)); f<-t(t(e)*d)
		}
		if (!ignore) r<-redrank(x,y-f,p) else r<-redrank(x,y,p)
		a<-r$a; b<-r$b; z<-tcrossprod(a,b); g<-x%*%z; 
		if (!ignore) h<-g+f else h<-g
		sa<-sb<-sum((y-h)^2)
		if (!ignore) {
			if (de) e<-ortcrus(t(t(y-g)*d),x)
				else {
					f<-regcrus(y-g,x); d<-sqrt(colSums(f*f)); e<-t(t(f)/d)
				}
			h<-g+f
			sb<-sum((y-h)^2)
		}
		cat("Iteration: ",formatC(itel,digits=6,width=6),
			"   Loss Start :   ",formatC(sold,digits=6,width=12,format="f"),
			"   Loss Coefs :   ",formatC(sa,digits=6,width=12,format="f"),
			"   Loss Scors: ",formatC(sb,digits=6,width=12,format="f"),"\n")
		if (((sold-sb)<eps) || (itel==itmax)) break()
		itel<-itel+1; sold<-sb
	}
if (!ignore) return(list(s=sb,e=e,d=d,a=a,b=b,z=z)) 
	else return(list(s=sb,a=a,b=b,z=z))
}

lsmimic<-function(x,y,p,eps=1e-6,itmax=100,fox=TRUE,eou=TRUE,foe=TRUE,df=TRUE,de=TRUE,ignore=FALSE) {}

lskalman<-function(){}

#
# procrus returns U of the same dimension as X, 
# such that tr U'X is maximized over all U satisfying U'U=I.
#

procrus<-function(x) {
	n<-nrow(x); m<-ncol(x)
	if (m > n) stop("matrix must be tall")
	s<-svd(x); u<-s$u; v<-s$v; k<-u[,1:m]; l<-v[,1:m] 
	return(tcrossprod(k,l))
}

#
# ortcrus returns U of the same dimension as X, 
# such that tr U'X is maximized,
# over all U satisfying U'U=I and U'Y=0.
#

ortcrus<-function(x,y) {
	return(procrus(x-qr.fitted(qr(y),x)))
}
	
#
# ortregs returns U of the same dimension as X, 
# such that SSQ(U-X) is minimized,
# over all U satisfying U'Y=0.
#

ortregs<-function(x,y) {
	return((x-qr.fitted(qr(y),x))
}

#
# ranort makes a random orthonormal matrix of dimension n x p
#

ranort<-function(n,p) {
	x<-matrix(rnorm(n*p),n,p)
	return(qr.Q(qr(x)))
}

projspace<-function(x,y) {}

projcone<-function(x,y) {}

redrank<-function(x,y,p) {
	q<-qr(x); c<-crossprod(y,qr.fitted(q,y))
	e<-eigen(c); b<-as.matrix(e$vectors[,1:p])
	a<-qr.coef(q,y)%*%b
	s<-sum((y-x%*%tcrossprod(a,b))^2)
	return(list(a=a,b=b,s=s))
}