function rosenbrock_nedler
% Derivative-free Nedler-Mead algorithm
%X = [[0; 2], [-0.5; 2], [-0.5; 1.5]];
X = [[0; 1], [1; 0], [0.499999999; 0.5]];
%X = randn(2,3);
N = size(X,2)-1;
% evaluate the function
F = zeros(N+1,1);
for i=1:N+1,
    F(i) = f(X(:,i));
end
% sort the points according to the function values
[X,F]=sort_points(X,F);

vizualize0;
HL=[];

Nmax = 3000;
for k=1:Nmax,
    [X,F]=nedler_step(X,F);
    fbar = mean(F);
    if(fbar < 1.0E-12), break; end
    printf('It: %3d, f: %15.6e\n', k, fbar);
    HL=vizualize1(X,F,HL); drawnow;
    pause(0.1);
end

function [X,F]=sort_points(X,F)
[F,I]=sort(F);
X=X(:,I);

function [X,F]=nedler_step(X,F)
N=length(F)-1;
xbar = mean(X(:,1:N),2);
xref = @(t) xbar + t*(X(:,N+1)-xbar);
x_m1 = xref(-1); f_m1 = f(x_m1);
if (F(1) <= f_m1) && (f_m1 < F(N))
    X(:,N+1) = x_m1; F(N+1) = f_m1;
    [X,F] = sort_points(X,F);
    return;
elseif f_m1 < F(1)
    x_m2 = xref(-2); f_m2 = f(x_m2);
    if f_m2 < f_m1
        X(:,N+1) = x_m2;  F(N+1) = f_m2;
        [X,F] = sort_points(X,F);
        return;
    else
        X(:,N+1) = x_m1; F(N+1) = f_m1;
        [X,F] = sort_points(X,F);
        return;
    end
elseif f_m1 >= F(N)
    if (F(N) <= f_m1) && (f_m1 < F(N+1))
        x_mh = xref(-0.5); f_mh = f(x_mh);
        if f_mh <= f_m1
            X(:,N+1) = x_mh; F(N+1) = f_mh;
            [X,F] = sort_points(X,F);
            return;
        end
    else
        x_ph = xref(0.5); f_ph = f(x_ph);
        if f_ph < F(N+1)
            X(:,N+1) = x_ph; F(N+1) = f_ph;
            [X,F] = sort_points(X,F);
            return;
        end
    end
    for i=2:N+1,
        X(:,i) = (X(:,1)+X(:,i))/2;
        F(i)   = f(X(:,i));
    end
    [X,F] = sort_points(X,F);
    return;
end

function vizualize0()
clf; hold on;
x1=linspace(-1.5,1.5,100);
y1=linspace(-1,2,100);
[x2,y2]=meshgrid(x1,y1);
f2=zeros(size(x2));
for i=1:length(x2(:)),
    f2(i) = f([x2(i),y2(i)]);
end
plot(1,1,'Color','r','Marker','*','MarkerSize',10);
fmax = max(f2(:)); fmin = 0;
p = 2;
z=linspace(fmin^(1/p),fmax^(1/p),50);
contour(x2,y2,f2,z.^p);
hold off;
HA=gca;

function [HL]=vizualize1(X,F,HL)
delete(HL);
hold on;
N = length(F)-1;
HL=line(X(1,[1:N+1,1]),X(2,[1:N+1,1]),'Color','k','Marker','x');
hold off;

function r=f(x)
r = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;

function r=df(x)
r=[-400*(x(2)-x(1)^2)*x(1) - 2*(1-x(1)); 200*(x(2)-x(1)^2)];

function r=df2(x)
r=[[400*(3*x(1)^2-x(2))+2, -400*x(1)];  [-400*x(1), 200]];