#
#   nlPCA 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, 2007-02-24,		first alpha
#

require("pava")
require("qpvarious")

nlpcaRank<-function(x,p) {
n<-dim(x)[1]; m<-dim(x)[2]; y<-norm(x); su<-1
repeat {
	r<-crossprod(y)/m; e<-eigen(r); v<-e$values; k<-e$vectors
	g<-k[,1:p]%*%t(k[,1:p])
	z<-y%*%g; zup<-z
	for (j in 1:m) {
	    zj<-z-outer(y[,j],g[j,])
		zz<-pava(zj[x[,j],j])[order(x[,j])]
		if (sum(zz^2)<1e-6) zup[,j]<-x[,j] else zup[,j]<-zz
		}	
	y<-norm(zup)
	print(cumsum(v[1:p])); sv<-sum(v[1:p])
	if (abs(sv-su)<1e-6) break
	su<-sv
	}
y
}

ordPCA<-function(data,ndim=2,eps=1e-6,deg=0,knots=c(),cca=FALSE,itmax=1000,movie=FALSE) {
n<-nrow(data); m<-ncol(data); sold<-Inf; itel<-1
data<-apply(data,2,function(y) y-mean(y))
data<-apply(data,2,function(y) y/sqrt(sum(y^2)))
opt<-data
repeat {
	optsvd<-svd(opt)
	d<-optsvd$d[1:ndim]; a<-matrix(optsvd$u[,1:ndim],n,ndim)
	b<-matrix(optsvd$v[,1:ndim]*outer(rep(1,m),d),m,ndim)
	if (movie) {
		pdf(paste("b",as.character(itel),".pdf",sep=""))
		plot(b,type="n",xlim=c(-3,3),ylim=c(-3,+3),
			xlab="dimension 1",ylab="dimension 2",main="")
		text(b,labels(data)[[1]],col="RED")
		dev.off()
		}
	ab<-tcrossprod(a,b)
	ssvd<-sum((opt-ab)^2)
	for (j in 1:m) {
		ind<-rank(data[,j]); ord<-order(ind); z<-ab[,j]
		if (deg == 0) zhat<-pava(z[ord])[ind]
		if (deg > 0) {
			if (length(knots) == 0) x<-outer(data[,j],0:deg,"^")
			if (length(knots) > 0) x<-bs(data[,j],knots=knots,degree=deg)
			if (cca)
				zhat<-(monRegCCA(z[ord],x[ord,],itout=FALSE)$yhat)[ind]
			if (!cca)
				zhat<-monTwoProj(z[ord],xord[ord,])[ind]
			}
		zhat<-zhat-mean(zhat); opt[,j]<-zhat/sqrt(sum(zhat^2))
		}
	snew<-sum((opt-ab)^2)
	cat("Iteration: ",formatC(itel,digits=3,width=3),
		"After SVD: ",formatC(ssvd,digits=6,width=10,format="f"),
		"After Regression: ",formatC(snew,digits=6,width=10,format="f"),
		"\n")
	if (((sold-snew) < eps) || (itel == itmax)) break()
		sold<-snew; itel<-itel+1
	}
return(list(a=a,b=b,opt=opt,s=snew,k=itel))
}


norm<-function(y) {
m<-dim(y)[2]; n<-dim(y)[1]
for (j in 1:m) {
	y[,j]<-y[,j]-(sum(y[,j])/n)
	y[,j]<-y[,j]/sqrt(sum(y[,j]^2))
	}
y
}