function XN=optimG(J,gradJ,bet,alphain,tau,X0,N)// dernier itéré (gradient)
 XN=X0;
 for i=1:N 
   alpha=alphain;d=-gradJ(XN);
   while(J(XN+alpha*d)>J(XN)+bet*alpha*d'*gradJ(XN))
      alpha=tau*alpha
   end
   XN=XN+alpha*d     
 end         
endfunction
//
function XN2=optimG2(J,gradJ,bet,alphain,tau,X0,N) // conserve l'historique (gradient)
 XN=X0;XN2=XN;
 for i=1:N 
   alpha=alphain;d=-gradJ(XN);
   while(J(XN+alpha*d)>J(XN)+bet*alpha*d'*gradJ(XN))
      alpha=tau*alpha
   end
   XN=XN+alpha*d   
   XN2=[XN2,XN]     
 end         
endfunction
//
function XN=optimH(J,gradJ,HJ,X0,N)// dernier itéré (Hessien)
 XN=X0;
 for i=1:N 
   alpha=1;d=-inv(HJ(XN))*gradJ(XN);
   XN=XN+alpha*d     
 end         
endfunction
//
function XN2=optimH2(J,gradJ,HJ,X0,N)// tout l'historique (Hessien)
 XN=X0;XN2=XN
 for i=1:N 
   alpha=1;d=-inv(HJ(XN))*gradJ(XN);
   XN=XN+alpha*d  
   XN2=[XN2,XN]     
 end         
endfunction
//

function XN=optimBFGS(J,gradJ,X0,N)// dernier itéré (BFGS)
 XN=X0;n=length(X0);B=eye(n,n)
 for i=1:N 
   Xold=XN;gxold=gradJ(XN);
   alpha=1;d=-inv(B)*gradJ(XN);
   XN=XN+alpha*d  
   y=(gradJ(XN)-gxold);s=XN-Xold;
   if (s<>0)&(y<>0) then
     B=B+(1/(s'*y))*y*y'-(1/(s'*(B*s)))*B*s*s'*B;  
   else
     B=B
   end      
 end         
endfunction
//
function XN2=optimBFGS2(J,gradJ,X0,N)// tout l'historique (BFGS)
 XN=X0;XN2=XN;n=length(X0);B=eye(n,n)
 for i=1:N 
   Xold=XN;gxold=gradJ(XN);
   alpha=1;d=-inv(B)*gradJ(XN);
   XN=XN+alpha*d  
   y=(gradJ(XN)-gxold);s=XN-Xold;
   if (s<>0)&(y<>0) then
     B=B+(1/(s'*y))*y*y'-(1/(s'*(B*s)))*B*s*s'*B;  
   else
     B=B
   end 
   XN2=[XN2,XN]   
 end     
endfunction
//
function y=rosen(x)
    y=100*(x(2)-x(1)^2)^2+(x(1)-1)^2
endfunction
//
function y=rosengrad(x)
y=[-400*x(1)*(x(2)-x(1)^2)+2*(x(1)-1);200*(x(2)-x(1)^2)]
endfunction
//
function y=rosenhess(x)
     y=[-400*x(2)+1200*x(1)^2+2, -400*x(1);-400*x(1),200]; 
endfunction
//
bet=0.1
alphain=1;
tau=0.3;
X0=[0;1];
//
N1=1000;
N2=5;
N3=36;
XNG=optimG(rosen,rosengrad,bet,alphain,tau,X0,N1)
XNG2=optimG2(rosen,rosengrad,bet,alphain,tau,X0,N1)
XNH=optimH(rosen,rosengrad,rosenhess,X0,N2)
XNH2=optimH2(rosen,rosengrad,rosenhess,X0,N2)
XNBFGS=optimBFGS(rosen,rosengrad,X0,N3)
XNBFGS2=optimBFGS2(rosen,rosengrad,X0,N3)

disp('nombre d''iterations gradient');disp(N1)
disp('valeur apres N iterations (rosen)');disp(XNG)
disp('nombre d''iterations hessien');disp(N2)
disp('valeur apres N iterations (rosen)');disp(XNH)
disp('nombre d''iterations BFHS');disp(N3)
disp('valeur apres N iterations BFGS (rosen)');disp(XNBFGS)
//
x=linspace(-1,2,100)
y=linspace(-1,2,100)
z=zeros(100,100)
for i=1:100
    for j=1:100
    z(i,j)=rosen([x(i),y(j)])
    end        
end
figure(1)
clf()
contour(x,y,z,0:20)
plot2d(XNG2(1,:),XNG2(2,:),2,rect=[-1,-1,2,2])
plot2d(XNH2(1,:),XNH2(2,:),5,rect=[-1,-1,2,2])
plot2d(XNBFGS2(1,:),XNBFGS2(2,:),3,rect=[-1,-1,2,2])
xtitle('méthode de gradient (bleu) vs mathode de newton (rouge) et BFGS (vert)')
//
valjG=zeros(N1+1,1)
valjH=zeros(N2+1,1)
valjBFGS=zeros(N3+1,1)
for i=1:N1+1
    valjG(i)=rosen(XNG2(:,i))
end
for i=1:N2+1
    valjH(i)=rosen(XNH2(:,i))
end
for i=1:N3+1
    valjBFGS(i)=rosen(XNBFGS2(:,i))
end
//
figure(2)
clf()

plot(log(1:(N1+1)),log(valjG))
plot(log(1:(N2+1)),log(valjH),'r-')
plot(log(1:(N3+1)),log(valjBFGS),'g-')
xtitle('méthode de gradient (bleu) vs mathode de newton (rouge), BFGS (vert) historique de convergence')
//
//
