%=====================================================
%
% OS82 - 2012/2013
% 
% Optimisation par Metaheuristique
% Algorithme du Simplexe (Nelder-Mead)
%
%=====================================================


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


%===== Représentations de la fonction de Rosenbrock 2D
xmin=-2;
ymin=-2;
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


%===== Algorithme du simplexe de Nelder-Mead
fprintf('\n***** Algorithme de Nelder-Mead\n');
N=2;        %dimension du pb
alpha=2;    %coefficient pour la reflexion
beta=3;     %coefficient pour l'expansion
gamma=1/2;  %coefficient pour la contraction

xi=-1.5;    %premier point du simplexe initial
yi=-1.5;
dl=0.2;     %longueur pour 

Simplexe=zeros(N+1,N);  %simplexe : N+1 points dans un espace de dimension N
Simplexe(1,1)=xi;       %1er point du simplexe initial
Simplexe(1,2)=yi;
Simplexe(2,1)=xi+dl;    %2ieme point du simplexe initial
Simplexe(2,2)=yi;
Simplexe(3,1)=xi;       %3ieme point du simplexe initial
Simplexe(3,2)=yi+dl;

Valeurs=zeros(3,1);     %valeurs de la fonction-objectif aux 3 points

figure(4)
hold on
contourf(X,Y,Z,100,'linestyle','none');
xlabel('x1');
ylabel('x2');
title('Algorithme du simplexe de Nelder-Mead');
axis equal;

nIterMax=200;
epsilonConvergence=1E-6;

test_fin=0;
compteur=0;

while(test_fin==0)

    Valeurs(1,1)=Rosenbrock2D(Simplexe(1,1),Simplexe(1,2),1);
    Valeurs(2,1)=Rosenbrock2D(Simplexe(2,1),Simplexe(2,2),1);
    Valeurs(3,1)=Rosenbrock2D(Simplexe(3,1),Simplexe(3,2),1);

    [vmax,jmax]=max(Valeurs);
    [vmin,jmin]=max(Valeurs);
    
    fprintf('* iteration %d\n',compteur);
    fprintf('- point 1 : x1=%f\tx2=%f\tf=%f\n',Simplexe(1,1),Simplexe(1,2),Valeurs(1,1));
    fprintf('- point 2 : x1=%f\tx2=%f\tf=%f\n',Simplexe(2,1),Simplexe(2,2),Valeurs(2,1));
    fprintf('- point 3 : x1=%f\tx2=%f\tf=%f\n',Simplexe(3,1),Simplexe(3,2),Valeurs(3,1));
    
    m=zeros(1,2);%barycentre de l'ensemble des points du simplexe prive de celui qui possede la valeur la plus elevee
    for i=1:N+1
        if i~=jmax
            m=m+Simplexe(i,:);
        end
    end
    m=m/N;

    xR=Simplexe(jmax,:)+alpha*(m-Simplexe(jmax,:)); %candidat par reflexion
    xE=Simplexe(jmax,:)+beta*(m-Simplexe(jmax,:));  %candidat par expansion
    xC=Simplexe(jmax,:)+gamma*(m-Simplexe(jmax,:)); %candidat par contraction

    ValeurR=Rosenbrock2D(xR(1,1),xR(1,2),1);%valeur par reflexion
    ValeurE=Rosenbrock2D(xE(1,1),xE(1,2),1);%valeur par expansion
    ValeurC=Rosenbrock2D(xC(1,1),xC(1,2),1);%valeur par contraction
    
    plot([Simplexe(:,1);Simplexe(1,1)],[Simplexe(:,2);Simplexe(1,2)],'-ko','markerfacecolor',[1,1,1]);
    fill(Simplexe(:,1),Simplexe(:,2),[1,1,1]);
    plot(m(1,1),m(1,2),'-ko','markerfacecolor',[0,0,0]);
    plot(xR(1,1),xR(1,2),'-ko','markerfacecolor',[1,0,0]);
    plot(xE(1,1),xE(1,2),'-ko','markerfacecolor',[0,1,0]);
    plot(xC(1,1),xC(1,2),'-ko','markerfacecolor',[0,0,1]);
	
    input('Presser ENTER pour continuer !');%pour executer pas a pas
    
    if(ValeurR<vmin)
        if(ValeurE<ValeurR)
            Simplexe(jmax,:)=xE;%expansion
        else
            Simplexe(jmax,:)=xR;%reflexion
        end
    else
        if(ValeurC<vmax)
            Simplexe(jmax,:)=xC;%contraction
        else
            SimplexeOld=Simplexe;
            for i=1:N+1
                Simplexe(i,:)=0.5*(SimplexeOld(jmin,:)+SimplexeOld(i,:));%contractions multiples
            end
        end
    end
    
    compteur=compteur+1;
    if(   (compteur==nIterMax)||(  ( max(Simplexe(:,1))-min(Simplexe(:,1))<epsilonConvergence )&&( max(Simplexe(:,2))-min(Simplexe(:,2))<epsilonConvergence )  )   )
        test_fin=1;
    end

end
hold off

fprintf('* Simplexe a l issue de la derniere iteration :\n',compteur);
fprintf('- point 1 : x1=%f\tx2=%f\tf=%f\n',Simplexe(1,1),Simplexe(1,2),Valeurs(1,1));
fprintf('- point 2 : x1=%f\tx2=%f\tf=%f\n',Simplexe(2,1),Simplexe(2,2),Valeurs(2,1));
fprintf('- point 3 : x1=%f\tx2=%f\tf=%f\n',Simplexe(3,1),Simplexe(3,2),Valeurs(3,1));