sviZajedno

function [x, fx] = goldenRatio(fun, a, b, tol)
    c = (3 - sqrt(5)) / 2;

    x1 = a + c * (b - a); f1 = feval(fun, x1);
    x2 = a + b - x1; f2 = feval(fun, x2);

    while abs(b - a) > tol
        if f1 <= f2
            b = x2;
            x2 = x1;
            f2 = f1;
            x1 = a + c * (b - a);
            f1 = feval(fun, x1);
        else
            a = x1;
            x1 = x2;
            f1 = f2;
            x2 = b - c * (b - a);
            f2 = feval(fun, x2);
        end

        if f1 < f2
            x = x1; fx = f1;
        else
            x = x2; fx = f2;
        end
    end
end


function [x, fx, it] = powell(fun, x0, epsilon)
    n = length(x0);
    u = eye(n);
    x = x0; x1 = x0 + 2 * epsilon;
    it = 0;

    while max(abs(x - x1)) > epsilon
        it = it + 1;
        ti = x;
        for i = 1:n
            g = @(s) fun(ti + s * u(:, i));
            teta = goldenRatio(g, -100, 100, epsilon);
            ti = ti + teta * u(:, i);
        end

        u = [u(:, i) ti - x];
        x1 = x;
        g = @(s) fun(ti + s * u(:, n));
        teta = goldenRatio(g, -100, 100, epsilon);
        x = ti + teta * u(:, n);
    end
    fx = feval(fun, x);
end


function [x, fx, cnt] = hookJeeves(fun, x0, d, dmin)
    n = length(x0);
    e = eye(n) * d;
    x = x0;
    fx = feval(fun, x);
    cnt = 1;

    while e(1, 1) > dmin
        t = x;
        for i = 1:n
            z = t + e(:, i);
            y = feval(fun, z);
            cnt = cnt + 1;
            if y < fx
            t = z;
            fx = y;
            else
                z = t - e(:, i);
                y = feval(fun, z);
                cnt = cnt + 1;
                if y < fx
                    t = z;
                    fx = y;
                end
            end
        end
        if all(t == x)
            e = e * 0.5;
        else
            x1 = t + (t - x);
            f1 = feval(fun, x1);
            cnt = cnt + 1;
            x = t;
            if f1 < fx
                for i = 1:n
                    z = x1 + e(:, i);
                    y = feval(fun, z);
                    cnt = cnt + 1;
                    if y < f1
                        x = z;
                        fx = y;
                        break;
                    end
                    z = x1 - e(:, i);
                    y = feval(fun, z);
                    cnt = cnt + 1;
                    if y < f1
                        x = z;
                        fx = y;
                        break;
                    end
                end
            end
        end
    end
end


function [ x, fx, cnt ] = NelderMead( fun, x, xtol, ftol, maxit )
    n = prod(size(x));
    maxit = maxit * n;

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

    for i = 1:n
        y = xin;
        if y(i) ~= 0
            y(i) = y(i) * 1.05;
        else
            y(i) = 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 = (1 + gamma) * vbar - gamma * v(:, n + 1);
                fe = feval(fun, ve); cnt = cnt + 1;
                if fe < fr
                    vk = ve; fk = fe;
                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 i = 2:n
                    v(:, i) = (v(:, i) + v(:, 1)) / 2;
                    fv(i) = feval(fun, v(:, i));
                end
            vk = (v(:, n + 1) + v(:, 1)) / 2;
            fk = feval(fun, vk); cnt = cnt + n;
            end
        end
        v(:, n + 1) = vk;
        fv(n + 1) = fk;
        [fv, i] = sort(fv);
        v = v(:, i);
    end

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