function [x,histout,costdata] = projbfgs(x0,f,up,low,tol,maxit)
%
% C. T. Kelley, June 11, 1998
%
% This code comes with no guarantee or warranty of any kind.
%
% function [x,histout,costdata] = projbfgs(x0,f,up,low,tol,maxit)
%
% projected BFGS with Armijo rule, simple linesearch 
% 
%
% Input: x0 = initial iterate
%        f = objective function,
%            the calling sequence for f should be
%            [fout,gout]=f(x) where fout=f(x) is a scalar
%              and gout = grad f(x) is a COLUMN vector
%        up = vector of upper bounds
%        low = vector of lower bounds
%        tol = termination criterion norm(grad) < tol
%              optional, default = 1.d-6
%        maxit = maximum iterations (optional) default = 1000
%
% Output: x = solution
%         histout = iteration history   
%             Each row of histout is
%   [norm(projgrad), f, number of step length reductions, iteration count,
%            relative size of active set]
%         costdata = [num f, num grad, num hess] (for steep, num hess=0)
%
if nargin < 4
error(' projbfgs requires bounds ');
end
xc=x0; ndim=length(up); kku=zeros(ndim,1); kkl=zeros(ndim,1);
%
% list of active indices 
%
alist=zeros(ndim,1);
%
for i=1:ndim
        kku(i)=up(i); kkl(i)=low(i);
	if kkl(i) > kku(i)
        error('lower bound exceeds upper bound')
        end
end
%
% put initial iterate in feasible set
%
if norm(xc - kk_proj(xc,kku,kkl)) > 0
      disp(' initial iterate not feasibile ');
      xc=kk_proj(xc,kku,kkl);
end
alp=1.d-4;
nsmax=5; ystore=zeros(ndim,nsmax); sstore=ystore; ns=0;
if nargin < 6
maxit=1000; 
end
if nargin < 5
tol=1.d-6;
end
itc=1; 
[fc,gc]=feval(f,xc); 
numf=1; numg=1; numh=0;
ithist=zeros(maxit,5);
pgc=xc - kk_proj(xc - gc,kku,kkl);
ia=0; alist=zeros(ndim,1);
tst=kku-kkl; lim1=.5*min(tst);
epsilon=min(lim1,norm(pgc));
for i=1:ndim; if(xc(i)==kku(i) | xc(i)==kkl(i)) ia=ia+1; end; end;
for i=1:ndim; 
         if( min(kku(i)-xc(i),xc(i)-kkl(i)) < epsilon) 
                 alist(i)=1; 
         end;
end;
ithist(1,5)=ia/ndim;
ithist(1,1)=norm(pgc); ithist(1,2) = fc; ithist(1,4)=itc-1; ithist(1,3)=0; 
while(norm(pgc) > tol & itc <= maxit)
        lambda=1;
        dsd=-gc; dsd=bfgsrp(ystore,sstore,ns,dsd,alist);
        dsd=dsd+proja(-gc,alist);
        xt=kk_proj(xc+lambda*dsd,kku,kkl); ft=feval(f,xt);
        numf=numf+1;
	iarm=0; itc=itc+1;
        pl=xc - xt; fgoal=fc-(gc'*pl)*alp;
%
%       simple line search
%
	while(ft > fgoal)
                lambda=lambda*.1;
		iarm=iarm+1;
                xta=xt;
		xt=kk_proj(xc+lambda*dsd,kku,kkl);
                pl=xc-xt;
		ft=feval(f,xt); numf = numf+1;
		if(iarm > 10) 
		disp(' Armijo error in steepest descent ')
                histout=ithist(1:itc,:); costdata=[numf, numg, numh];
                x=xc;
		return; end
                fgoal=fc-(gc'*pl)*alp;
	end
	[fc,gp]=feval(f,xt); numf=numf+1; numg=numg+1;
        y=gp-gc; s=xt-xc;
        gc=gp; xc=xt; pgc=xc-kk_proj(xc-gc,kku,kkl); 
        epsilon=min(lim1,norm(pgc));
        alist=zeros(ndim,1); ial=0;
                for i=1:ndim; 
                     if( min(kku(i)-xc(i),xc(i)-kkl(i)) < epsilon) 
                         alist(i)=1; ial=ial+1; 
                     end;
                end;
        y=proji(y,alist); s=proji(s,alist); 
%
%   restart if y'*s is not positive or we're out of room
%
        yts=y'*s;
        if yts <=0 ns=0;
          end
        if ns==nsmax 
            ns=0;
        elseif yts > 0
                ns=ns+1;
                alpha=sqrt(yts);
                sstore(:,ns)=s/alpha; ystore(:,ns)=y/alpha;
        end
	ithist(itc,1)=norm(pgc); ithist(itc,2) = fc; 
	ithist(itc,4)=itc-1; ithist(itc,3)=iarm;
ia=0; for i=1:ndim; if(xc(i)==kku(i) | xc(i)==kkl(i)) ia=ia+1; end; end;
ithist(itc,5)=ia/ndim;
end
x=xc; 
histout=ithist(1:itc,:); costdata=[numf, numg, numh];
%
% projection onto active set
%
function px = kk_proj(x,kku,kkl)
ndim=length(x);
px=zeros(ndim,1);
px=min(kku,x); 
px=max(kkl,px);
%
% bfgsrp
%
% C. T. Kelley, Dec 20, 1996
%
% This code comes with no guarantee or warranty of any kind.
%
% This code is used in projbfgs.m
%
% There is no reason to ever call this directly.
%
% form the product of the generalized inverse of the
% bfgs approximate Hessian
% with a vector using the recursive approach
%
function dnewt=bfgsrp(ystore,sstore,ns,dsd,alist)
dnewt=proji(dsd,alist);
if (ns==0)
return;
end;
sstore(:,ns)=proji(sstore(:,ns),alist);
ystore(:,ns)=proji(ystore(:,ns),alist);
beta=sstore(:,ns)'*dsd; dnewt=dsd-beta*ystore(:,ns);
ndim=length(dsd); xlist=zeros(ndim,1);
dnewt=bfgsrp(ystore,sstore,ns-1,dnewt,xlist);
dnewt=dnewt+(beta-ystore(:,ns)'*dnewt)*sstore(:,ns); 
dnewt=proji(dnewt,alist);
%
% projection onto epsilon-inactive set
%
function px=proji(x,alist)
ndim=length(x); px=x;
for k=1:ndim if alist(k) == 1 px(k)=0; end; end
%
%  projection onto epsilon-active set
%
function px=proja(x,alist)
ndim=length(x); px=x;
for k=1:ndim if alist(k) == 0 px(k)=0; end; end