#
#   esopo 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-04-06   Initial Alpha Release
# version 0.2.00, 2007-02-17   Added monRegCCA()
#

ccaGen<-function(y,a,b=rep(0,dim(a)[1]),w=diag(length(y)),eps=1e-6,verbose=FALSE) {
x<-y; l<-rep(0,length(b)); m<-dim(a)[1]; r<-((a%*%x)-b); g<-solve(w,t(a))
c<-a%*%g; d<-diag(c)
repeat{
    ops<-0
	for (i in 1:m) {
		t<-max(-l[i],r[i]/d[i])
		l[i]<-l[i]+t; x<-x-t*g[,i]; r<-r-t*c[,i]
	    ops<-max(ops,t)
		}
    if (ops < eps) break()
	}
list(x=x,l=l,r=r,ssq=sum(x^2))
}

ccaPairs<-function(y,pairs,b=rep(0,dim(a)[1]),w=diag(length(y)),eps=1e-6,verbose=FALSE) {
x<-y; n<-dim(pairs)[1]
repeat{
    ops<-0
	for (i in 1:n) {
	    k<-pairs[i,1]; l<-pairs[i,2]
		t<-max(0,(x[k]-x[l])/2)
		x[k]<-x[k]-t; x[l]<-x[l]+t
	    ops<-max(ops,t)
		}
    if (ops < eps) break()
	}
list(x=x,ssq=sum(x^2))
}

ccaMono<-function(y,eps=1e-6) {
x<-y; n<-length(y)-1
repeat{
    ops<-0
	for (i in 1:n) {
		t<-max(0,(x[i]-x[i+1])/2)
		x[i]<-x[i]-t; x[i+1]<-x[i+1]+t
	    ops<-max(ops,t)
		}
    if (ops < eps) break()
	}
list(x=x,ssq=sum(x^2))
}

ccaTriangular<-function(d,eps=1e-6) {
n<-dim(d)[1]; dhat<-d
repeat{
	ops<-0
	for (i in 1:(n-2)) {
		for (j in i:(n-1)) {
			for (k in j:n) {
			t<-max(0,dhat[i,j]-(dhat[i,k]+dhat[j,k]))/3
			dhat[i,j]<-dhat[i,j]-t; dhat[i,k]<-dhat[i,k]+t; dhat[j,k]<-dhat[j,k]+t
			ops<-max(ops,t)
			}
		}
	}
    if (ops < eps) break()
    print(sum((d-dhat)^2))
	}
list(dhat=dhat,ssq=sum((d-dhat)^2))
}

mkDiffMat<-function(n) {
a<-diag(n)[1:(n-1),]
return(a+ifelse(outer(1:(n-1),1:n,"-")==-1,-1,0))
}

monRegCCA<-function(y,x=diag(length(y)),a=diff(diag(length(y))),eps=1e-6,itmax=100,verbose=FALSE,inter=TRUE){
qr<-qr(x); q<-qr.Q(qr); r<-qr.R(qr); m<-nrow(a); n<-ncol(x)
b<-colSums(y*q); aq<-a%*%q; lbd<-rep(0,m); gam<-b; yhat<-drop(q%*%b)
aa<-rowSums(aq^2); itel<-1; yy<-sum(y^2); sold<-yy-sum(b^2)
repeat {
    snew<-sold
	for (i in 1:m) {
		aqi<-aq[i,]; aai<-aa[i]
		lopt<-max(0,-sum(gam*aqi))/aai
		lbd[i]<-lbd[i]+lopt
		gam<-gam+lopt*aqi
		snew<-snew+aai*(lopt^2)
		if (verbose) 
			cat("Cycle: ",formatC(itel,digits=3,width=3),
				"Coordinate: ",formatC(i,digits=3,width=3),
				"Loss: ",formatC(snew,digits=6,width=10,format="f"),
				"\n")
		}
	if (inter)
		cat("Cycle: ",formatC(itel,digits=3,width=3),
			"****************",
			"Previous Loss: ",formatC(sold,digits=6,width=10,format="f"),
			"Current Loss: ",formatC(snew,digits=6,width=10,format="f"),
			"\n")
	if (((snew-sold) < eps) || (itel == itmax)) break()
	sold<-snew; itel<-itel+1
	}
return(list(yhat=q%*%gam,s=snew,k=itel))
}