//
// Methode d'optimisation sans gradient
//
//   Methode Nelder Mead
//----------------------------------------
// definition de la fonction ˆ minimiser
//
//deff('y=f(x)','y=sum(x.*x)');                  // fonction quadratique
//deff('y=f(x)','y=sum(x.*x-cos(2*%pi*x))+n');   // fonction de Rastrigin
function y=f(x)                                     // fonction Rosenbrock
    n=length(x)
y=100*(x(n)-x(n-1)^2)^2+norm(x-ones(x))^2-(x(n)-1)^2;
endfunction
//
//
// un cas qui ne marche pas pour Nelder Mead (Mac Kinnon)
//deff('y=f(x)','if x(1)<0 then y=360*x(1)^2+x(2)+x(2)^2;else y=6*x(1)^2+x(2)+x(2)^2;end');
//X=[1,1;(1+sqrt(33))/8,(1-sqrt(33))/8;0,0]; // simplexe initial
n=evstr(x_dialog('nombre de parametres (contours si n=2)','2'));
Nit=evstr(x_dialog('nombre d''iterations','2000'));
//
xmin=-2;xmax=2;

//
// trace des lignes de niveau
//
if (n==2)
N=300;
xplot=xmin:((xmax-xmin)/(N-1)):xmax;
yplot=xplot;
zplot=zeros(N,N);
for i=1:N
for j=1:N
zplot(i,j)=f([xplot(i),yplot(j)]);
end
end
xset('window',0)
clf()
contour2d(xplot,yplot,zplot,20);
end
// 
// creation du premier simplexe
//
X=xmin+(xmax-xmin)*rand(n+1,n+1);
// un cas qui ne marche pas pour Nelder Mead (Mac Kinnon)
//X=[1,1;(1+sqrt(33))/8,(1-sqrt(33))/8;0,0]; // simplexe initial
//
//
// evaluation du premier simplexe
for i=1:n+1
    X(i,n+1)=f(X(i,1:n))
end
Neval=n+1;  
//
// debut boucle principale
for k=1:Nit
  sh=0;
  // ordonnement du simplexe
  [a,b]=gsort(-X(:,n+1));
  X=X(b,:);
  //
  xc=1/n*sum(X(1:n,1:n),'r');
  d=xc-X(n+1,1:n);
  xr=X(n+1,1:n)+2*d;
  fr=f(xr);fnew=fr;Neval=Neval+1;
  // cas 1
  if (X(1,n+1)<fr)&(fr<X(n,n+1)) then
    xnew=xr
  end  
  // cas 2
  if (fr<X(1,n+1)) then
    xe=X(n+1,1:n)+3*d;fe=f(xe);
    if (fe<fr) then xnew=xe;fnew=fe;Neval=Neval+1;else xnew=xr;end;
  end
  // cas 3
  if (X(n,n+1)<fr)&(fr<X(n+1,n+1)) then
    xtilde=X(n+1,1:n)+3/2*d;ftilde=f(xtilde);
    if (ftilde<fr) then xnew=xtilde;fnew=ftilde;Neval=Neval+1;else xnew=xr;end;
  end
  // cas 4
  if (X(n+1,n+1)<fr) then
    xhat=X(n+1,1:n)+1/2*d;fhat=f(xhat);
    if (fhat<X(n+1,n+1)) then xnew=xhat;fnew=fhat;Neval=Neval+1;else 
      //    retrecissement
      sh=1;
      for j=2:n+1
        X(j,1:n)=X(1,1:n)+1/2*(X(j,1:n)-X(1,1:n));
        X(j,n+1)=f(X(j,1:n));
      end   
      Neval=Neval+n; 
    end;
  end  
  // ajout du nouveau point
  //
  if (sh==0) then
    X(n+1,1:n)=xnew;X(n+1,n+1)=fnew;
  end  
  //
  //trace si n=2
  if (n==2) then
    plot2d([X(:,1);X(1,1)],[X(:,2);X(1,2)],k,rect=[xmin,xmin,xmax,xmax]);
  end  
end
disp('meilleur point:'); disp(X(1,1:n));
disp('valeur de f en ce point:'); disp(X(1,n+1));
disp('nombre d''evaluations:');disp(Neval)
//
