// méthode gradient a pas constant et  methode de Newton 
// appliquée à la fonction de Rosenbrock
//
// TP2 optimisation 2018
//
function y=R(x)  // Rosenbrock
    y=10*(x(2)-x(1)^2)^2+(x(1)-1)^2
endfunction
//
function y=nablaR(x)  // gradient de Rosenbrock
    y=[-40*x(1)*(x(2)-x(1)^2)+2*(x(1)-1);20*(x(2)-x(1)^2)];
endfunction
//
function y=HR(x) // Hessien de Rosenbrock
    y=[-40*(x(2)-x(1)^2)+80*x(1)^2+2,-40*x(1);-40*x(1),20];
endfunction
//
function [utot, nb_iter] = GPF2 (J,nablaJ,u0,epsilon,max_iter,rho)
nb_iter=0;
utot=u0;resid=nablaJ(u0);u=u0;   
while(nb_iter<max_iter) & (norm(resid)>epsilon)
  u=u-rho*resid;
  utot=[utot,u]
  resid=nablaJ(u); 
  nb_iter=nb_iter+1 
end  
endfunction
//
function [utot, nb_iter] = newton(J,nablaJ,HJ,u0,epsilon,max_iter)
nb_iter=0;
utot=u0;resid=nablaJ(u0);u=u0;   
while(nb_iter<max_iter) & (norm(resid)>epsilon)
  u=u-inv(HJ(u))*nablaJ(u);
  utot=[utot,u]
  resid=nablaJ(u); 
  nb_iter=nb_iter+1 
end  
endfunction
//
//////////////// affichage des résultats///////////////////
u0=[1.1;0.1];
epsilon=1E-4
max_iter=20000;
rho1=0.01;rho2=0.02;rho3=0.03;
[utot1, nb_iter1] = GPF2 (R,nablaR,u0,epsilon,max_iter,rho1)
disp('valeur finale obtenue rho=0.01');disp(utot1(:,$))
disp('nombre d''iterations rho=0.01');disp(nb_iter1)
[utot2, nb_iter2] = GPF2 (R,nablaR,u0,epsilon,max_iter,rho2)
disp('valeur finale obtenue rho=0.01');disp(utot2(:,$))
disp('nombre d''iterations rho=0.02');disp(nb_iter2)
[utot3, nb_iter3] = GPF2 (R,nablaR,u0,epsilon,max_iter,rho3)
disp('valeur finale obtenue rho=0.01');disp(utot3(:,$))
disp('nombre d''iterations rho=0.03');disp(nb_iter3)
[utot4, nb_iter4] = newton(R,nablaR,HR,u0,epsilon,max_iter)
//
////////////sortie graphique /////////////////////////////
//
x=-1:0.03:2;y=-1:0.03:2;
z=zeros(length(x),length(y))
for i=1:length(x)
    for j=1:length(y)
      z(i,j)=R([x(i);y(j)])  
    end  
end
//
figure(1)
clf()
subplot(221)
contour(x,y,z,[0:0.2:2,3:10])
xtitle('lignes de niveau Rosenbrock')
subplot(222)
contour(x,y,z,10)
plot2d(utot1(1,:),utot1(2,:),1)
xtitle('trajection, rho=0.01')
subplot(223)
contour(x,y,z,10)
plot2d(utot2(1,:),utot2(2,:),2)
xtitle('trajection, rho=0.02')
subplot(224)
contour(x,y,z,10)
plot2d(utot3(1,:),utot3(2,:),3)
xtitle('trajection, rho=0.03')
//
figure(2)
clf()
err1=utot1-ones(utot1);
normerr1=sqrt(err1(1,1:100).^2+err1(2,1:100).^2)
plot2d(1:100,log(normerr1),1)
//
err2=utot2-ones(utot2);
normerr2=sqrt(err2(1,1:100).^2+err2(2,1:100).^2)
plot2d(1:100,log(normerr2),2)
//
err3=utot3-ones(utot3);
normerr3=sqrt(err3(1,1:100).^2+err3(2,1:100).^2)
plot2d(1:100,log(normerr3),3)
//
err4=utot4-ones(utot4);
normerr4=sqrt(err4(1,:).^2+err4(2,:).^2)
plot2d(1:length(normerr4),log(normerr4),4)
legends(['gradient 1';'gradient 2';'gradient 3';'newton'],[1,2,3,4])
xtitle('vitesse de convergence (gradient vs Newton)','Niter','log(erreur)')
