%Federica Tinti (ftinti@dm.unipd.it)
%Department of Mathematics of Ferrara (ITALY)
%
%Hyperplane projection method (Solodov-Svaiter)
%
clear all
feasible=0;         %feasible starting point         
%feasible=1;        %unfeasible starting point
stringa=input('test problem=', 's');
eval(stringa);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%parameters
P=eye(n);
eps=1e-4;   
eta=1.;
gamma=0.5;
sigma=0.3;
theta=4;
mu=min([theta*eta,1]);
np=0;
%
%unfiseable starting point
if feasible==1
    options=optimset('LargeScale','off','display','off');
    f=zeros(size(x));
    [x]=quadprog(P,f,A,b,Aeq,beq,lb,ub,x,options);
    x=abs(x);
    np=np+1;
end;
x0=x;
%
nf=1;           
fx=feval(fs,x);
f=-(x-mu*fx);
options=optimset('LargeScale','off','display','off');
[xx]=quadprog(P,f,A,b,Aeq,beq,lb,ub,x,options);
xx=abs(xx);
np=np+1;
rx=x-xx;
iter=0;
tempo=cputime;
while norm(rx)>eps | iter==0
   iter=iter+1;
   k=-1;
   while k==-1 |sini< destra
      k=k+1;
      zz=x-gamma^k*mu*rx;
      fzz=feval(fs,zz);
      nf=nf+1;   
      sini=fzz'*rx;
      destra=sigma/mu*(rx'*rx);
   end;
   eta=gamma^k*mu;
   z=x-eta*rx;
   fz=feval(fs,z);
   nf=nf+1;
   f=-x;
   Aa=[A;fz'];
   ba=[b; fz'*z];
   options=optimset('LargeScale','off','display','off');
   [x]=quadprog(P,f,Aa,ba,Aeq,beq,lb,ub,x,options);
   x=abs(x);
   np=np+1;
   fx=feval(fs,x);
   nf=nf+1;
   mu=min([theta*eta,1]);
   f=-(x-mu*fx);
   options=optimset('LargeScale','off','display','off');
   [xx]=quadprog(P,f,A,b,Aeq,beq,lb,ub,x,options);
   xx=abs(xx);
   np=np+1;
   rx=x-xx;
   if iter==1 
      x0=x;
   end;     
end;
t=cputime-tempo
fprintf('solution=%g \n',x);          
fprintf('number of projections=%g \n',np);          
fprintf('number of function evalutions=%g \n',nf);         
fprintf('number of iterations=%g \n',iter);