#
#   dykstra 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-24   Added monProj()
#

dHalfSpace<-function(x,a,b,w=rep(1,length(x)),test="lt",verbose=FALSE,eps=1e-8,itmax=100) {
p<-length(b); n<-length(x); xcr<-x; y<-matrix(0,p,n); itel<-1; fold<-0
repeat {
	xpr<-xcr
	for (i in 1:p) {
		xup<-projHalf(xpr-y[i,],a[i,],b[i],test,w)
		y[i,]<-xup-(xpr-y[i,])
		xpr<-xup
		}
	fnew<-sum(w*(x-xup)^2); fcyc<-sum(w*(xcr-xup)^2)
	if (verbose) 
		cat("Cycle: ",formatC(itel,width=3, format="d"),
	    	" Change: ", formatC(fcyc,digits=8,width=12,format="f"),
			" SSQ: ", formatC(c(fold,fnew),digits=8,width=12,format="f"),"\n")
	if (((fnew-fold) < eps) || (itel == itmax))
		return(list(f=fnew,x=xup))
	xcr<-xup; fold<-fnew; itel<-itel+1
	}
}

projHalf<-function(x,a,b,test="lt",w=rep(1,length(x))) {
xa<-sum(x*a); aw<-a/w; xb<-xa-b; aa<-sum(a*aw)
if ((test == "lt") && (xb < 0)) return(x)
if ((test == "gt") && (xb > 0)) return(x)
return(x-aw*(xb/aa))
}

dMonReg<-function(x,w=rep(1,length(x)),test="lt",verbose=FALSE,off=0,eps=1e-6,itmax=100) {
n<-length(x); b<-rep(off,n-1); a<--diff(diag(length(x)))
dHalfSpace(x,a,b,w,test,verbose,eps,itmax)
}

#
# monProj(y,x) projects a vector y on the intersection of the subspace spanned by x
# and the cone of increasing vectors using alternating projections
#

monProj<-function(y,x,verbose=FALSE,eps=1e-6,itmax=1000){
z<-y; sold<-Inf; itel<-1
	repeat{
		b<-colSums(z*x)
		u<-drop(x%*%b)
		smid<-sum((z-u)^2)
		z<-pava(u)
		snew<-sum((z-u)^2)
		if (verbose) 
			cat("Cycle: ",formatC(itel,digits=3,width=3),
				"After Subspace: ",formatC(smid,digits=6,width=10,format="f"),
				"After Cone: ",formatC(snew,digits=6,width=10,format="f"),
				"\n")
		if ((max(smid,snew) < eps) || (itel == itmax)) break()
		itel<-itel+1
		}
return(z)
}