%=====================================================
%
% OS82 - 2012/2013
% 
% Optimisation par descente avec gradient
%
%=====================================================


%===== Nettoyage de Matlab
clc;
clear all;
close all;


%===== Représentations de la fonction de Rosenbrock 2D
xmin=-1;
ymin=-1;
xmax=1.5;
ymax=1.5;
nx=20;
ny=20;

X=zeros(nx,1);
Y=zeros(ny,1);
for i=1:nx
    X(i,1)=xmin+(xmax-xmin)*(i-1)/(nx-1);
end
for i=1:ny
    Y(i,1)=ymin+(ymax-ymin)*(i-1)/(ny-1);
end

Z=zeros(ny,nx);%attention a l'ordre des indices - convention Matlab !
for i=1:ny
    y=Y(i,1);
    for j=1:nx
        x=X(j,1);
        Z(i,j)=Rosenbrock2D(x,y,1);
    end
end

figure(1)
hold on
surf(X,Y,Z,'linestyle','none');
xlabel('x1');
ylabel('x2');
title('Fonction de Rosenbrock 2D, k=1');
hold off

figure(2)
hold on
contourf(X,Y,Z,100,'linestyle','none');%,'linewidth',2);
xlabel('x1');
ylabel('x2');
title('Fonction de Rosenbrock 2D lissée, k=1');
axis equal;
hold off

%===== Representation du gradient
gradX=zeros(ny,nx);
gradY=zeros(ny,nx);
for i=1:ny
    y=Y(i,1);
    for j=1:nx
        x=X(j,1);
        temp=GradRosenbrock2D(x,y,1);
        gradX(i,j)=temp(1,1);
        gradY(i,j)=temp(2,1);
    end
end

figure(3)
hold on
contourf(X,Y,Z,100,'linestyle','none');%,'linewidth',2);
quiver(X,Y,gradX,gradY,2,'color',[1,0,0])
xlabel('x1');
ylabel('x2');
title('Gradient de la Fonction de Rosenbrock 2D lissée, k=1');
axis equal;
hold off


%===== Algorithme de descente à pas constant
fprintf('\n***** Algorithme de descente à pas constant\n');
xi=-0.5;
yi=1;
alpha=0.05;

gradient=zeros(2,1);

figure(4)
hold on
contourf(X,Y,Z,100,'linestyle','none');%,'linewidth',2);
quiver(X,Y,gradX,gradY,2,'color',[1,0,0])
xlabel('x1');
ylabel('x2');
title('Algorithme de descente à pas constant');
axis equal;

nIterMax=100;
epsilonConvergence=1E-5;

test_fin=0;
compteur=0;
fi=Rosenbrock2D(xi,yi,1);
while(test_fin==0)

    plot(xi,yi,'go','markerfacecolor',[0,1,0]);
    fprintf('* iteration %d\n(xi,yi)=\t\t(%f,%f)\t fi=\t%f\n',compteur,xi,yi,fi);
    
    gradient=GradRosenbrock2D(xi,yi,1);
    gradient=gradient/norm(gradient);
    xi=xi-alpha*gradient(1,1);
    yi=yi-alpha*gradient(2,1);
    fip1=Rosenbrock2D(xi,yi,1);
    
    fprintf('(xip1,yip1)=\t(%f,%f)\t fip1=\t%f\n',xi,yi,fip1);
    
    compteur=compteur+1;
    if((compteur==nIterMax)||(abs(fip1-fi)<epsilonConvergence))
        test_fin=1;
    end
    fi=fip1;
end
hold off

input('presser ENTER pour continuer');


%===== Algorithme de descente à pas variable 1
fprintf('\n***** Algorithme de descente à pas variable 1\n');
xi=-0.5;
yi=1;
alpha0=0.2;

gradient=zeros(2,1);

figure(5)
hold on
contourf(X,Y,Z,100,'linestyle','none');%,'linewidth',2);
quiver(X,Y,gradX,gradY,2,'color',[1,0,0])
xlabel('x1');
ylabel('x2');
title('Algorithme de descente à pas variable 1');
axis equal;

nIterMax=1000;
epsilonConvergence=1E-8;

test_fin_1=0;
compteur=0;
fi=Rosenbrock2D(xi,yi,1);
while(test_fin_1==0)
    plot(xi,yi,'go','markerfacecolor',[0,1,0]);
    fprintf('* iteration %d\n(xi,yi)=\t\t(%f,%f)\t fi=\t%f\n',compteur,xi,yi,fi);
    gradient=GradRosenbrock2D(xi,yi,1);
    gradient=gradient/norm(gradient);
    
    alpha=alpha0;
    test_fin_2=0;
    while (test_fin_2==0)
        xi_test=xi-alpha*gradient(1,1);
        yi_test=yi-alpha*gradient(2,1);
        fip1=Rosenbrock2D(xi_test,yi_test,1);
        if(fip1<fi)
            xi=xi_test;
            yi=yi_test;
            test_fin_2=1;
        else
            alpha=alpha/2;
        end 
    end
    
    fprintf('(xip1,yip1)=\t(%f,%f)\t fip1=\t%f\n',xi,yi,fip1);
    
    compteur=compteur+1;
    if((compteur==nIterMax)||(abs(fip1-fi)<epsilonConvergence))
        test_fin_1=1;
    end
    
    fi=fip1;
end
hold off