function mckinnon_nedler
% Derivative-free Nedler-Mead algorithm
lambda_1 = (1+sqrt(33))/8;
lambda_2 = (1-sqrt(33))/8;
% initial simplex
X = [[0;0], [lambda_1^1; lambda_2^1], [lambda_1^2; lambda_2^2]];
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(0,-1/2,'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)
%tau = 1.0; theta = 15.0; phi = 10.0;
%tau = 2.0; theta = 6.0; phi = 60.0;
tau = 3.0; theta = 6.0; phi = 400.0;
if x(1)<0,
  r = theta*phi*abs(x(1))^tau + x(2) + x(2)^2;
else
r = theta*x(1)^tau + x(2) + x(2)^2;
end