function trabalho2()
    ponto_inicial = [0,1];
    r = 10;
    beta = 0.1;
    tol = 1e-6;


    %[minimo,iteracoes,tempo] = barreira(ponto_inicial,r,beta,tol,@univariante,0)
    %f1(minimo(1),minimo(2))
    %[minimo,iteracoes,tempo] = barreira(ponto_inicial,r,beta,tol,@powell,0)
    %f1(minimo(1),minimo(2))
    %[minimo,iteracoes] = barreira(ponto_inicial,r,beta,tol,@maxdeclive,1)
    %f1(minimo(1),minimo(2))
    %[minimo,iteracoes] = barreira(ponto_inicial,r,beta,tol,@fletcher_reeves,1)
    %f1(minimo(1),minimo(2))
    %[minimo,iteracoes] = barreira(ponto_inicial,r,beta,tol,@BFGS,1)
    %f1(minimo(1),minimo(2))
    [minimo,iteracoes,tempo] = barreira(ponto_inicial,r,beta,tol,@newton_raphson,2)
    f1(minimo(1),minimo(2))


    %{
    [X,Y]=meshgrid(0.8:0.01:1.1);
    contourf(X,Y,f(X,Y))
    hold on
    L = 0.8:0.1:1.1;
    plot(L,teste(L))
    %}
end

function y = g(L)
    global g_f;
    global g_x;
    global g_d;
    y = g_f(L*g_d(1) + g_x(1),L*g_d(2) + g_x(2));
end

function [minimo,iteracoes,tempo] = barreira(ponto_inicial,r,beta,tol,fun,ordem)
    global g_r;
    passo = 0.001;
    j = 1;
    tic;
    for k = 1:100
        g_r = r;
        if ordem == 0
            [min,iter] = fun(ponto_inicial,@phi,passo);
        elseif ordem == 1
            [min,iter] = fun(ponto_inicial,@phi,@grad_phi,passo);
        else
            [min,iter] = fun(ponto_inicial,@phi,@grad_phi,@hess_phi,passo);
        end
        min
        iter;
        converg = r*b(min(1), min(2))
        if converg < 0 || iter < 0
            passo = passo / 10;
            continue;
        else
            %passo = 0.001;
        end
        j = j + 1;
        if abs(converg) < tol
            minimo = min;
            iteracoes = j;
            break
        else
           r = r*beta;
           ponto_inicial = min;
        end
    end
    toc;
    tempo = toc;
end

function y = f1(x1,x2)
    y = (x1 - 2).^4 + (x1 - 2*x2).^2;
end

function y = c(x1,x2)
    y = x1.^2 - x2;
end

function y = b(x1,x2)
    y = -1./c(x1,x2);
end

function y = phi(x1,x2)
    global g_r;
    y = f1(x1,x2) + g_r*b(x1,x2);
end

function y = grad_phi(x1,x2)
    global g_r;
    g11 = 2*g_r*x1/(x1.^2 - x2).^2 + 2*(x1 - 2*x2) + 4*(x1 - 2).^3;
    g12 = -g_r/(x1.^2 - x2).^2 - 4*x1 + 8*x2;
    y = [g11,g12];
end

function y = hess_phi(x1,x2)
    global g_r;
    h11 = 2*g_r/(x1.^2 - x2).^2 - (8*g_r*(x1.^2))/(x1.^2 - x2).^3 + 12*(x1 - 2).^2 + 2;
    h12 = (4*g_r*x1)/(x1.^2 - x2).^3 - 4;
    h21 = (4*g_r*x1)/(x1.^2 - x2).^3 - 4;
    h22 = 8 - 2*g_r/(x1.^2 - x2).^3;
    y = [h11 h12;h21 h22];
end

function [u,iter] = univariante(x0,f,passo)
    global g_f;
    g_f = f;
    global g_x;
    global g_d;
    tol = 1e-6;
    tol2 = 1e-7;
    min_d1 = x0;
    min_d2 = x0;

    %para plotar a rota de minimização
    listax(1) = x0(1);
    listay(1) = x0(2);
    j = 2;

    for k = 1:50000
        if (norm(min_d1 - min_d2) <= tol) && (k > 1)
            u = min_d2;
            iter = k;
            return;

        else
            x0 = min_d2;
            d = [1,0];
            g_x = x0;
            g_d = d;
            [a1,b1] = passoconstante(0, @g, passo, 1);
            c1 = secaoaurea(a1, b1, tol2, @g);
            min_d1 = x0 + d * c1;


            d = [0,1];
            x0 = min_d1;
            g_x = x0;
            g_d = d;
            [a2,b2] = passoconstante(0, @g, passo, 1);
            c2 = secaoaurea(a2, b2, tol2, @g);
            min_d2 = x0 + d * c2;

            %para plotar a rota de minimização
            % listax(j) = min_d1(1);
            % listax(j+1) = min_d2(1);
            % listay(j) = min_d1(2);
            % listay(j+1) = min_d2(2);
            % j = j + 2;

        end
    end

    u = min_d2;
    iter = -1;

    %{
    [X,Y]=meshgrid(-5:0.01:5);
    contourf(X,Y,f1(X,Y))
    hold on
    plot(2,2,'r*')
    plot(x0(1),x0(2),'r*')
    plot(listax, listay, 'red')
    %}
end

function [min,iter] = powell(x0,f,passo)
    global g_f;
    g_f = f;
    global g_x;
    global g_d;

    tol = 1e-6;
    tol2 = 1e-7;

    %para plotar a rota de minimização
    %listax(1) = x0(1);
    %listay(1) = x0(2);
    %j = 2;

    u = {[1,0], [0,1]};
    n = 2;

    for k = 1:1000

        p0 = x0;
        for j=1:n
            g_x = x0;
            g_d = u{j};
            [a1,b1] = passoconstante(0, @g, passo, 1);
            c1 = secaoaurea(a1, b1, tol2, @g);
            %sprintf("%d %d %f %f %f %f %f", k, j, x0(1), x0(2), u{j}(1), u{j}(2), c1)
            x0 = g_x + g_d * c1;
        end

        for j=1:n-1
            u{j} = u{j+1};
        end
        u{n} = x0 - p0;

        g_x = x0;
        g_d = u{n};
        [a1,b1] = passoconstante(0, @g, passo, 1);
        c1 = secaoaurea(a1, b1, tol2, @g);
        %sprintf("%d %d %f %f %f %f %f", k, 3, x0(1), x0(2), u{n}(1), u{n}(2), c1)
        x0 = g_x + g_d * c1;
        %sprintf("%f %f %f %f", x0(1), x0(2), p0(1), p0(2))

        if abs(x0 - p0) < tol
            min = x0;
            iter = k;
            return
        end

        if mod(k,n+1) == 0
            %u = {[1,0], [0,1]};
        end

        % L = -2:0.0001:2;
        % plot(L,g(L))
        % ylim([-50,50])
        % pause


        %para plotar a rota de minimização
        %listax(j) = min_d1(1);
        %listax(j+1) = min_d2(1);
        %listay(j) = min_d1(2);
        %listay(j+1) = min_d2(2);
        %j = j + 2;


    end

    iter = -1;

    %[X,Y]=meshgrid(-5:0.01:5);
    %contourf(X,Y,f1(X,Y))
    %hold on
    %plot(2,2,'r*')
    %plot(x3(1),x3(2),'r*')
    %plot(listax, listay, 'red')
end

function [min,iter] = maxdeclive(x0,f,grad_f,passo)
    global g_f;
    g_f = f;
    global g_x;
    global g_d;

    tol = 1e-6;
    tol2 = 1e-7;

    for k = 1:5000
        x = x0;
        d = -grad_f(x0(1),x0(2));
        if norm(d) <= tol
            min = x0;
            iter = k;
            return
        end
        g_x = x;
        g_d = d;
        [a,b] = passoconstante(0, @g, passo, 1);
        c = secaoaurea(a, b, tol2, @g);
        min = [x(1) + d(1)*c, x(2) + d(2)*c];
        x0 = min;
    end
    iter = -1;
    min = x0;
end

function [min,iter] = fletcher_reeves(x0,f,grad_f,passo)
    global g_f;
    g_f = f;
    global g_x;
    global g_d;

    tol = 1e-6;
    tol2 = 1e-7;
    x = x0;
    grad_x0 = grad_f(x0(1),x0(2));
    d = -grad_x0;

    %para plotar a rota de minimização
    %{
    listax(1) = x0(1);
    listay(1) = x0(2);
    j = 2;
    %}

    %tic;
    for k = 1:100
        if norm(d) <= tol
            min = x;
            iter = k;
            break
        end
        g_x = x;
        g_d = d;
        [a,b] = passoconstante(0, @g, passo, 1);
        c = secaoaurea(a, b, tol2, @g);
        %L = c-5:0.1:c+5;
        %plot(L,g(L))
        x1 = [x(1) + d(1)*c, x(2) + d(2)*c];
        grad_x1 = grad_f(x1(1),x1(2));
        if norm(grad_x1) <= tol
            min = x1;
            iter = k;
            break
        end
        beta = (grad_x1*(grad_x1.'))/(grad_x0*(grad_x0.'));
        d = -grad_x1 + beta.*d;
        x = x1;
        grad_x0 = grad_x1;

        %para plotar a rota de minimização
        %{
        listax(j) = x1(1);
        listax(j+1) = x(1);
        listay(j) = x1(2);
        listay(j+1) = x(2);
        j = j + 2;
        %}
    end

    %{
    toc;
    t = toc;

    [X,Y]=meshgrid(-5:0.01:5);
    contourf(X,Y,f1(X,Y))
    hold on
    plot(2,2,'r*')
    plot(x1(1),x1(2),'r*')
    plot(listax, listay, 'red')
    %}
end

function [min,iter] = BFGS(x0,f,grad_f,passo)
    global g_f;
    g_f = f;
    global g_x;
    global g_d;
    tol = 1e-6;
    tol2 = 1e-7;
    x = x0;
    S0 = [1 0;0 1];
    g_x0 = grad_f(x0(1),x0(2));
    d = -(S0*g_x0.').';

    %para plotar a rota de minimização
    %{
    listax(1) = x(1);
    listay(1) = x(2);
    j = 2;

    tic;
    %}
    for k = 1:100
        if norm(d) <= tol
            min = x;
            iter = k;
            break
        end
        g_x = x;
        g_d = d;
        [a,b] = passoconstante(0, @g, passo, 1);
        c = secaoaurea(a, b, tol2, @g);
        x1 = [x(1) + d(1)*c, x(2) + d(2)*c];
        g_x1 = grad_f(x1(1),x1(2));
        if norm(g_x1) <= tol
            min = x1;
            iter = k;
            break
        end
        d_xk = (x1 - x0).';
        d_gk = (g_x1 - g_x0).';
        primeiro_termo = (((d_xk.')*d_gk + (d_gk.')*S0*d_gk)*d_xk*(d_xk.'))/((d_xk.')*d_gk)^2;
        segundo_termo = (S0*d_gk*(d_xk.') + d_xk*((S0*d_gk).'))/((d_xk.')*d_gk);
        S1 = S0 + primeiro_termo - segundo_termo;

        d = -(S1*g_x1.').';
        x = x1;
        S0 = S1;
        g_x0 = g_x1;

        %para plotar a rota de minimização
        %{
        listax(j) = x(1);
        listax(j+1) = x1(1);
        listay(j) = x(2);
        listay(j+1) = x1(2);
        j = j + 2;
        %}
    end

    %{
    toc;
    t = toc;

    [X,Y]=meshgrid(-5:0.01:5);
    contourf(X,Y,f1(X,Y))
    hold on
    plot(2,2,'r*')
    plot(x1(1),x1(2),'r*')
    plot(listax, listay, 'red')
    %}
end

function [min,iter] = newton_raphson(x0,f,grad_f,hess_f,passo)
    global g_f;
    g_f = f;
    global g_x;
    global g_d;
    tol = 1e-5;
    tol2 = 1e-5;
    x = x0;
    grad_x0 = grad_f(x0(1),x0(2));
    h = hess_f(x0(1),x0(2));
    inv_h = inv(h);
    d = -(inv_h*grad_x0.').';

    %para plotar a rota de minimização
    %{
    listax(1) = x(1);
    listay(1) = x(2);
    j = 2;
    %}
    %tic;

    for k = 1:1000
        if norm(d) <= tol
            min = x;
            iter = k;
            return
        end
        g_x = x;
        g_d = d;
        [a,b] = passoconstante(0, @g, passo, 1);
        c = secaoaurea(a, b, tol2, @g);
        x1 = [x(1) + d(1)*c, x(2) + d(2)*c];
        grad_x1 = grad_f(x1(1),x1(2));
        if norm(grad_x1) <= tol
            min = x1;
            iter = k;
            %listax(j) = x1(1);
            %listay(j) = x1(2);
            return
        end

        h = hess_f(x1(1),x1(2));
        inv_h = inv(h);
        d = -(inv_h*grad_x1.').';
        x = x1;

        %para plotar a rota de minimização
        %{
        listax(j) = x(1);
        listax(j+1) = x1(1);
        listay(j) = x(2);
        listay(j+1) = x1(2);
        j = j + 2;
        %}
    end

    iter = -1;
    min = x;
    %toc;
    %t = toc;

    %{
    [X,Y]=meshgrid(-5:0.01:5);
    contourf(X,Y,f1(X,Y))
    hold on
    plot(2,2,'r*')
    plot(x1(1),x1(2),'r*')
    plot(listax, listay, 'red')
    %}
end

function [y,w] = passoconstante(x0,f,da,d)
    % y0 = f(x0);

    % for k = 1:10000000
    %     x1 = x0 + da*d;
    %     y1 = f(x1);
    %     if k == 1
    %         if y1 > y0
    %             d = -d;
    %             continue
    %         elseif y1 == y0
    %             y = min([x0 - da*d,x0 + da*d]);
    %             w = max([x0 - da*d,x0 + da*d]);
    %             break
    %         end
    %     else
    %         if y1 >= y0
    %             y = min([x0 - da*d,x0 + da*d]);
    %             w = max([x0 - da*d,x0 + da*d]);
    %             break
    %         end
    %     end
    %     x0 = x1;
    %     y0 = y1;
    % end

    x1 = x0;
    h = da * d;
    c = 1.618033989;
    f1 = f(x1);
    x2 = x1 + h; f2 = f(x2);
    % Determine downhill direction & change sign of h if needed.
    if f2 > f1
        h = -h;
        x2 = x1 + h; f2 = f(x2);
        % Check if minimum is between x1 - h and x1 + h
        if f2 > f1
            y = x2; w = x1 - h; return
        end
    end
    % Search loop
    for i = 1:100
        h = c*h;
        x3 = x2 + h; f3 = f(x3);
        if f3 > f2
            y = x1; w = x3; return
        end
        x1 = x2; f1 = f2; x2 = x3; f2 = f3;
    end
end

function y = secaoaurea(aL,aU,tol,f)
    Ra = (sqrt(5)-1)/2;

    da = aU - aL;

    aE = (1-Ra)*da + aL;
    aD = Ra*da + aL;
    faE = f(aE);
    faD = f(aD);

    for k = 1:5000
        if abs(da) <= tol
            y = (aL + aU)/2;
            break
        end

        if faE >= faD
            aL = aE;
            da = aU - aL;

            aE = aD;
            faE = faD;
            aD = Ra*da + aL;
            faD = f(aD);
        else
            aU = aD;
            da = aU - aL;

            aD = aE;
            faD = faE;
            aE = (1-Ra)*da + aL;
            faE = f(aE);
        end
    end
end