function [x,histout,costdata] = gradproj(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] = gradproj(x0,f,up,low,tol,maxit)
%
% gradient projection 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(grad), 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(' gradproj requires bounds ');
end
xc=x0; ndim=length(up); kku=zeros(ndim,1); kkl=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;
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);
xt=kk_proj(xc - gc,kku,kkl);
pgc=xc - kk_proj(xc - gc,kku,kkl);
ia=0; for i=1:ndim; if(xc(i)==kku(i) | xc(i)==kkl(i)) ia=ia+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;
        xt=kk_proj(xc-lambda*gc,kku,kkl); ft=feval(f,xt);
        numf=numf+1;
	iarm=0; itc=itc+1;
        pl=xc - xt; fgoal=fc-(pl'*pl)*(alp/lambda);
%
%       simple line search
%
        q0=fc; qp0=-gc'*gc; qc=ft;
	while(ft > fgoal)
                lambda=lambda*.1;
		iarm=iarm+1;
		xt=kk_proj(xc-lambda*gc,kku,kkl);
                pl=xc-xt;
		ft=feval(f,xt); numf = numf+1;
		if(iarm > 10) 
		disp(' Armijo error in gradient projection')
                histout=ithist(1:itc,:); costdata=[numf, numg, numh];
		return; end
                fgoal=fc-(pl'*pl)*(alp/lambda);
	end
	xc=xt; [fc,gc]=feval(f,xc); numf=numf+1; numg=numg+1;
        pgc=xc-kk_proj(xc-gc,kku,kkl); 
	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);