%Federica Tinti (ftinti@dm.unipd.it)
%Department of Mathematics of Ferrara (ITALY)
%
%Extragradient method with adaptive steplegth 
%       (second modified version)
%
clear all
%
stringa=input('test problem=', 's');
eval(stringa);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%parameters
P=eye(n);
eps=1e-4;
rx=0;
icont=0;
icont=icont+1;
alpha(icont)=1.;
alphacap=1e-3;
%
beta=0.7;     %beta=0.5; %beta=0.8;
csi=0.8;      %csi=0.5;
flag=0;
%
iter=0;      
nf=0;               
np=0;				
icontvisual=0;
tic;
while norm(rx)>eps | iter==0 
    if iter==0
        fx=feval(fs,x);
        nf=nf+1;
        f=-(x-alpha(icont)*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;        
        fxx=feval(fs,xx);
        nf=nf+1;
        diff=x-xx;
        fdiff=fx-fxx;
        num=fdiff'*fdiff;
        den=diff'*diff;
        while (alpha(icont)>beta*sqrt(den/num))   
             icont=icont+1;
             alpha(icont)=max(alphacap,min(alpha(icont-1)*csi,(beta*sqrt(den/num))));                        
             f=-(x-alpha(icont)*fx);   
             options=optimset('LargeScale','off','display','off');
             [xx]=quadprog(P,f,A,b,Aeq,beq,lb,ub,x,options);  
             xx=abs(xx);
             fxx=feval(fs,xx);
             np=np+1;
             nf=nf+1;
             diff=x-xx;
             fdiff=fx-fxx;
             num=fdiff'*fdiff;
             den=diff'*diff;
        end;    
    else
          icont=icont+1;
          alpha(icont)=beta*norm(xnew-xx)/norm(fxnew-fxx);
          x=xnew;
          fx=fxnew;
          f=-(x-alpha(icont)*fx);
          options=optimset('LargeScale','off','display','off');
          [xx]=quadprog(P,f,A,b,Aeq,beq,lb,ub,x,options);
          xx=abs(xx);
          fxx=feval(fs,xx);
          nf=nf+1;
          if (alpha(icont) > 0.9*norm(x-xx)/norm(fx-fxx))
              fprintf('non verifica');
              pause;
          end;    
      end;
        %second projection
        iter=iter+1;
        alphab(iter)=alpha(icont);
        f=-(x-alpha(icont)*fxx);
        options=optimset('LargeScale','off','display','off');
        [xnew]=quadprog(P,f,A,b,Aeq,beq,lb,ub,xx,options);  
        xnew=abs(xnew);
        np=np+1;
        fxnew=feval(fs,xnew);
        nf=nf+1;
        rx=xnew-x;
  end;
  toc;
  t=toc;
  fprintf('number of projections=%g \n',np);          
  fprintf('number of function evaluations=%g \n',nf);         
  fprintf('number of iterations=%g \n',iter);           
  fprintf('input parameter beta=%g \n',beta); 
  fprintf('input parameter csi=%g \n',csi); 