// Methode d'optimisation sans gradient
//
//   MŽthode MDS
//----------------------------------------
// 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
//
n=evstr(x_dialog('nombre de parametres (contours si n=2)','2'));
Nit=evstr(x_dialog('nombre d''iterations','1000'));
//
xmin=-2;xmax=2;
//
function Y=evalsimplexe(X)
  [n1,n2]=size(X);
  Y=X;
  for i=2:n1
    Y(i,n2)=f(X(i,1:n2-1));
  end
  [a,b]=gsort(-Y(:,n2));
  Y=Y(b,:);
endfunction  
//
// tracŽ 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);
//X=[4,5,0;5,6,0;5,3,0]
//
// 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
  // ordonnement du simplexe
  [a,b]=gsort(-X(:,n+1));
  X=X(b,:);
  //
  Xr=X;
  Xr(2:n+1,1:n)=2*ones(n,1)*X(1,1:n)-X(2:n+1,1:n);
  //
  Xr=evalsimplexe(Xr);
  //
  if (Xr(1,n+1)<X(1,n+1)) then 
    Xe=X;
    Xe(2:n+1,1:n)=3*ones(n,1)*X(1,1:n)-2*X(2:n+1,1:n);   
    Xe=evalsimplexe(Xe);Neval=Neval+n;
    if (Xe(1,n+1)<Xr(1,n+1)) then Xnew=Xe, else Xnew=Xr;end
  else
    Xc=X;
    Xc(2:n+1,1:n)=0.5*ones(n,1)*X(1,1:n)+0.5*X(2:n+1,1:n);   
    Xc=evalsimplexe(Xc);Neval=Neval+n;
    Xnew=Xc;
  end   
  //
  //tracŽ 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
  X=Xnew;  
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)
    
