close all; clear; clc; format short;

%%
f = @(x,y) log(1 + (x.*sign(real(x))).^(2 + sin(x+y))) + y.^2/20;
ff = @(p) f(real(p),imag(p)); % helper f - called with complex coordinates

tol = 1e-4; reflection = 1; contraction = 1/2;

p = 8.3+1i*(-8.9); % Starting point p0
P = p + exp(1i*(0:2)*2*pi/3); F = ff(P); % Initial simplex
for j = 1:1000 % Limit number of iterations to 1e3
    [~,l] = min(F); [~,h] = max(F); % Find highest and lowest points
    Pc = P; Pc(h) = []; Pc = mean(Pc); % Find centroid of h-less P
    Pn = (1+reflection)*Pc - reflection*P(h); Fn = ff(Pn); % Reflect
    if any(Fn < F)
        P(h) = Pn; F(h) = Fn; % Accept if any better
    else
        Pn = contraction*P(h) + (1-contraction)*Pc; Fn = ff(Pn); % Contract
        if Fn < F(h) 
            P(h) = Pn; F(h) = Fn; % Accept if better than F(h)
        else
            P = (P+P(l))/2; F = ff(P); % Expand about P(l)
        end
    end    
    if max( abs(diff([P P(1)])) ) < tol; break; end % Check for convergence
end

pf = mean(P); pf = [real(pf) imag(pf)];
[j pf(1) pf(2) f(pf(1),pf(2))]

