hookjeev / neld

% HOOK-JEEV METODA:
function [x, f, cnt] = HookJeev(fun, x0, d, dmin)
% HookJeev metod
% x0 - Pocetno pogadjanje
% d - Pocetni korak
% dmin - Krajnja duzina koraka (uslov zaustavljanja)
% x - Minimum funkcije
% f - Vrijednost minimuma
% cnt - Broj poziva funkcije
n = length(x0);     % Broj promenljivih
e = eye(n) * d;     % Pravci kretanja
x = x0;             % Prva tacka
f = feval(fun, x);
cnt = 1;

while e(1, 1) > dmin
    t = x;
    for i = 1:n
        % Pozitivan smijer
        z = t + e(:, i);
        y = feval(fun, z);
        cnt = cnt + 1;

        if y < f
            t = z;
            f = y;
        else
            % Negativan smijer
            z = t - e(:, i);
            y = feval(fun, z);
            cnt = cnt + 1;

            if y < f
                t = z;
                f = y;
            end
        end
    end

    if all(t == x)
        e = e * 0.5;
    else
        x1 = t + (t - x);
        y1 = feval(fun, x1);
        cnt = cnt + 1;
        x = t;

        if y1 < f
            for i = 1:n
                z = x1 + e(:, i);
                y = feval(fun, z);
                cnt = cnt + 1;

                if y < y1
                    x = z;
                    f = y;
                    break; % Imamo novu tacku!
                end
                z = x1 - e(:, i);
                y = feval(fun, z);
                cnt = cnt + 1;

                if y < y1
                    x = z;
                    f = y;
                    break;
                end
            end
        end
    end
end
end

% PLOT:
[x1, x2] = meshgrid(-5:.5:15,-30:.5:10); % -5 : 15 (desna strana), -30 : 10 (lijeva strana) dijagrama
y = fun2(x1, x2);
surf(x1,x2,y)

% FUN2.M:
function [fx] = fun2(x1, x2)
% fx =x1.^2+x1.*x2+0.5*x2.^2+x1+10*x2;
fx = 2.*x1.^2 + x2.^2-3;
end

% FUN.M:
function [fx] = fun(x)
% fx =x(1).^2+x(1).*x(2)+0.5*x(2).^2+x(1)+10*x(2);
fx = 2.*x(1).^2 + x(2).^2-3;
end

% POZIV HOOK-JEEV:
d = 1; % valjda
x0 = [0; 0];
dmin = 0.5; % valjda
[x, f, cnt] = HookJeev('fun', x0, d, dmin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% NELDER-MEAD METODA:
%
function [x, fx, cnt] = NeldMead(fun, x, xtol, ftol, maxit)
% NelderMead metod
% x - Pocetna tacka
% xtol, ftol, maxit - Kriterijum zaustavljanja

n = prod(size(x));
maxit = maxit * n;

% Pocetni simplex:
xin = x(:);
v = xin;
fv = feval(fun, v);

for j = 1:n
    y = xin;

    if y(j) ~= 0
        y(j) = 1.05 * y(j);
    else
        y(j) = 0.00025;
    end

    v = [v y];
    fv = [fv feval(fun, y)];
end

[fv, j] = sort(fv);
v = v(:, j);
cnt = n + 1;
alpha = 1;
beta = 1/2;
gamma = 2;
onesn = ones(1, n);
ot = 2:n+1;
on = 1:n;

while cnt < maxit
    if max(max(abs(v(:, ot) - v(:, onesn)))) <= xtol && max(abs(fv(1) - fv(ot))) <= ftol
        break;
    end

    vbar = (sum(v(:, on)') / n)';
    vr = (1 + alpha) * vbar - alpha * v(:, n+1);
    fr = feval(fun, vr);
    cnt = cnt + 1;
    vk = vr;
    fk = fr;

    if fr < fv(n)
        if fr < fv(1)
            ve = gamma * vr + (1 - gamma) * vbar;
            fe = feval(fun, ve);
            cnt = cnt + 1;

            if fe < fv(1)
                vk = ve;
                fk = fr;
            end
        end
    else
        vt = v(:, n + 1);
        ft = fv(n + 1);

        if fr < ft
            vt = vr;
            ft = fr;
        end

        vc = beta * vt + (1 - beta) * vbar;
        fc = feval(fun, vc);
        cnt = cnt + 1;

        if fc < fv(n)
            vk = vc;
            fk = fc;
        else
            for j = 2:n
                v(:, j) = (v(:, 1) + v(:, j)) / 2;
                fv(j) = feval(fun, v(:, j));
            end

            vk = (v(:, 1) + v(:, n + 1)) / 2;
            fk = feval(fun, vk);
            cnt = cnt + n;
        end
    end

    v(:, n + 1) = vk;
    fv(n + 1) = fk;
    [fv, j] = sort(fv);
    v = v(:, j);
end

x(:) = v(:, 1);
fx = fv(1);

if cnt == maxit
    disp(['Upozorenje: Maksimalan broj iteracija (',int2str(maxit),') je premasen!']);
end
end

% PLOT:
% prilagoditi koordinate!
[x1, x2] = meshgrid(-5:.5:15,-30:.5:10); % -5 : 15 (desna strana), -30 : 10 (lijeva strana) dijagrama
y = fun2(x1, x2);
surf(x1,x2,y)

% FUN2.M:
function [fx] = fun2(x1, x2)
% fx =x1.^2+x1.*x2+0.5*x2.^2+x1+10*x2;
fx = 2.*x1.^2 + x2.^2-3;
end

% FUN.M:
function [fx] = fun(x)
% fx =x(1).^2+x(1).*x(2)+0.5*x(2).^2+x(1)+10*x(2);
fx = 2.*x(1).^2 + x(2).^2-3;
end

% POZIV NELDMEAD:
xtol = 1; % valjda
ftol = 1; % valjda
x0 = [-2; -2];
maxit = 200;
[x, fx, cnt] = NeldMead('fun', x0, xtol, ftol, maxit)