function nmm

% Nelder-Mead algorithm

%  Inputs:
%        x = starting vector
%        xa, xb = x-interval used in contour plot
%        ya, yb = y-interval used in contour plot
%        tol = tolerance for stopping iteration

%  Required (end of file)
%        F    is function F(x,y) to be minimized

tol=0.0001;

%  Specify which example to run
ifun=1

clf
hold on
% get(gcf)
set(gcf,'Position', [4 431 560 420])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

if ifun==1
    ns=50;
    x1(1)=-1;   x1(2)=-1;
    x2(1)=-1;   x2(2)=-0.72;
    x3(1)=-0.7; x3(2)=-1;
    xa=-1.1; xb=1.1; ya=-1.1; yb=1.1;
    xd=linspace(xa,xb,50);
    yd=linspace(ya,yb,50);
    [X,Y]=meshgrid(xd,yd);
    Z=100*(X.^2-Y).^2+(X-1).^2;
    exact=[1, 1];
    hold on
    set(gca,'ytick',[-1 -0.5 0 0.5 1]);
    set(gca,'xtick',[-1 -0.5 0 0.5 1]);
elseif ifun==2
    ns=5;
    x1(1)=-1;   x1(2)=-1;
    x2(1)=-1;   x2(2)=-0.7;
    x3(1)=-0.7; x3(2)=-1;
    xa=-1.1; xb=1.1; ya=-1.1; yb=1.1;
    xd=linspace(xa,xb,50);
    yd=linspace(ya,yb,50);
    [X,Y]=meshgrid(xd,yd);
    Z=(X-Y).^4+8*X.*Y-X+Y+3;
    exact=[.55357993584438379685, -.55357993584438379685];  
    hold on
    set(gca,'ytick',[-1 -0.5 0 0.5 1]);
    set(gca,'xtick',[-1 -0.5 0 0.5 1]);
elseif ifun==3
    ns=4;
    x1(1)=-1;   x1(2)=-1;
    x2(1)=-1;   x2(2)=-0.4;
    x3(1)=-0.5; x3(2)=-1;
    xa=-1.1; xb=1.1; ya=-1.1; yb=1.1;
    xd=linspace(xa,xb,200);
    yd=linspace(ya,yb,200);
    [X,Y]=meshgrid(xd,yd);
    Z=10*X.^2+Y.^2;
    exact=[0, 0];
    hold on
    set(gca,'ytick',[-1 -0.5 0 0.5 1]);
    set(gca,'xtick',[-1 -0.5 0 0.5 1]);
end

% DD = [F value, x value, y value] sorted by values of F (small to big)
DD = sortrows( [ [F(x1(1),x1(2),ifun), x1(1), x1(2) ];
    [F(x2(1), x2(2),ifun), x2(1), x2(2) ];
    [F(x3(1), x3(2),ifun) , x3(1), x3(2)]]);

contour(X,Y,Z,30,'b','LineWidth',1)
xpoints1=[DD(1,2), DD(2,2)]; ypoints1=[DD(1,3), DD(2,3)];
xpoints2=[DD(2,2), DD(3,2)]; ypoints2=[DD(2,3), DD(3,3)];
xpoints3=[DD(3,2), DD(1,2)]; ypoints3=[DD(3,3), DD(1,3)];
box on
if ifun==1
    plot(exact(1),exact(2),'r*','MarkerSize',14,'LineWidth',2)
end
plot(xpoints1,ypoints1,'ro-','LineWidth',2)
plot(xpoints2,ypoints2,'ro-','LineWidth',2)
plot(xpoints3,ypoints3,'ro-','LineWidth',2)

xlabel('x-axis','FontSize',16,'FontWeight','bold');
ylabel('y-axis','FontSize',16,'FontWeight','bold');

say=['1'];
xx=(DD(1,2)+DD(2,2)+DD(3,2))/3;
yy=(DD(1,3)+DD(2,3)+DD(3,3))/3;
text(xx,yy,say,'FontSize',24,'FontWeight','bold')

set(gca,'FontSize',16,'FontWeight','bold')

pause

% surfc(X,Y,Z)
% xlabel('x-axis','FontSize',14,'FontWeight','bold');
% ylabel('y-axis','FontSize',14,'FontWeight','bold');
% zlabel('z-axis','FontSize',14,'FontWeight','bold');
% set(gca,'FontSize',14,'FontWeight','bold');
% pause

alpha=1;  beta=0.5;
for is=1:ns
    DD = sortrows( DD);
    
    xB=DD(1,2)-DD(3,2);  yB=DD(1,3)-DD(3,3);
    xC=DD(2,2)-DD(3,2);  yC=DD(2,3)-DD(3,3);
    area=0.5*(xB*yC-xC*yB);
    if abs(area)<tol
        break
    end
    
    cen=0.5*[ DD(1,2)+DD(2,2), DD(1,3)+DD(2,3) ];
    r=cen+alpha*[ cen(1)-DD(3,2), cen(2)-DD(3,3)];
    Fr=F(r(1),r(2),ifun);
    
    if Fr<DD(2,1) & DD(1,1) <= Fr
        DD(3,1)=Fr; DD(3,2)=r(1); DD(3,3)=r(2);
    elseif Fr<DD(1,1)
        rr=cen+2*alpha*[ cen(1)-DD(3,2), cen(2)-DD(3,3)];
        Frr=F(rr(1),rr(2),ifun);
        if Frr<Fr
            DD(3,1)=Frr; DD(3,2)=rr(1); DD(3,3)=rr(2);
        else
            DD(3,1)=Fr; DD(3,2)=r(1); DD(3,3)=r(2);
        end
    elseif Fr<DD(3,1)
        rA=cen+beta*[ cen(1)-DD(3,2), cen(2)-DD(3,3)];
        FA=F(rA(1),rA(2),ifun);
        if FA <= Fr
            DD(3,1)=FA; DD(3,2)=rA(1); DD(3,3)=rA(2);
        else
            DD(2,2)=DD(1,2)+beta*(DD(2,2)-DD(1,2));
            DD(2,3)=DD(1,3)+beta*(DD(2,3)-DD(1,3));
            DD(2,1)=F(DD(2,2),DD(2,3),ifun);
            DD(3,2)=DD(1,2)+beta*(DD(3,2)-DD(1,2));
            DD(3,3)=DD(1,3)+beta*(DD(3,3)-DD(1,3));
            DD(3,1)=F(DD(3,2),DD(3,3),ifun);
        end
        % step 6
    else
        rA=[DD(3,2), DD(3,3)]+beta*[ cen(1)-DD(3,2), cen(2)-DD(3,3)];
        FA=F(rA(1),rA(2),ifun);
        if FA<DD(3,1)
            DD(3,1)=FA; DD(3,2)=rA(1); DD(3,3)=rA(2);
        else
            DD(2,2)=DD(1,2)+beta*(DD(2,2)-DD(1,2));
            DD(2,3)=DD(1,3)+beta*(DD(2,3)-DD(1,3));
            DD(2,1)=F(DD(2,2),DD(2,3),ifun);
            DD(3,2)=DD(1,2)+beta*(DD(3,2)-DD(1,2));
            DD(3,3)=DD(1,3)+beta*(DD(3,3)-DD(1,3));
            DD(3,1)=F(DD(3,2),DD(3,3),ifun);
        end
    end
    
    xpoints1=[DD(1,2), DD(2,2)]; ypoints1=[DD(1,3), DD(2,3)];
    xpoints2=[DD(2,2), DD(3,2)]; ypoints2=[DD(2,3), DD(3,3)];
    xpoints3=[DD(3,2), DD(1,2)]; ypoints3=[DD(3,3), DD(1,3)];
    plot(xpoints1,ypoints1,'ro-','LineWidth',2)
    plot(xpoints2,ypoints2,'ro-','LineWidth',2)
    plot(xpoints3,ypoints3,'ro-','LineWidth',2)
    
    iss=is+1;
    say=[num2str(iss)];
    xx=(DD(1,2)+DD(2,2)+DD(3,2))/3;
    yy=(DD(1,3)+DD(2,3)+DD(3,3))/3;
    text(xx,yy,say,'FontSize',24,'FontWeight','bold')
    
   
    pause
end


% error=sqrt((exact(1)-xx)^2+(exact(2)-yy)^2);
% fprintf('\n\n  Number of Steps = %i',counter)
% fprintf('\n  Computed Solution =  %e  %e',xx,yy);
% fprintf('\n  Error = %e\n\n',error);


function z=F(x,y,ifun)
if ifun==1
    z=100*(x^2-y)^2+(x-1)^2;
elseif ifun==2
    z=(x-y)^4+8*x*y-x+y+3;
elseif ifun==3
    z=10*x^2+y^2;
end;