% NewtonBFGS.m Implemenation of Quasi-Newton's Method with BFGS
%
% Usage: modify the starting point, and function pointer below.
%
% Author: Sven Leyffer, Argonne National Laboratory, 2016.
% Modification on NewtonExact.m by Dev, Meenarli, and Prashant on Sep 13,
% 2016.
%-----------------------------------------------------------------------

% ... initialization: set initial point, and assing function pointer
n        = 2;
x        = [-1.5;0.5];
%x        = [0;0];
%fun      = @rosenbrock;
x        = [1.5;0.5];
fun      = @powell;
PrntAlg  = 1;      % ... indicate whether or not to print the progress

% ... set up print of progress ... CHOOSE FUNCTIO BELOW
% if (PrntAlg == 1)
%     figure(1); clf;
%     % ... print rosenbrock
%     %% xx = -2:0.04:2;
%     %% [xx,yy] = meshgrid(xx,xx);
%     %% ff = 100*(xx.^2 - yy).^2 + (xx - 1).^2;
%     %% contour(xx,yy,ff,40); hold on;
%     %% title('Progress of Newton's method on Rosenbrock Function');
%     % ... print powell
    xx = -2:0.04:2;
    [xx,yy] = meshgrid(xx,xx);
    ff = xx.^4   + xx.*yy + (1+yy).^2;
%     contour(xx,yy,ff,15); hold on;
%     title('Progress of Newton method on Powell Function');
%    
% end

% ... define the constants of the search
sigma  = 0.1;      % ... sufficient reduction parameter
OptTol = 1e-4;     % ... loose optimality tolerance
kmax   = 1000;     % ... maximum number of iterations
fail   = 0;        % ... failure code for Armijo search
MinEig = 1e-2;     % ... minimum perturbation to eigenvalues

% ... set up output format
format1 = '%5i %12.6g %12.6g %12.6g %12.6g %5i \n';
fprintf(' iter    objective      || g ||   min(eig(H))       alpha LiFail\n');

% ... initialization
k        = 0;
H        = eye(n);
[f,g]    = feval( fun, x );
fstart   = f;
xstart   = x;
s        = zeros(n,1);

while ( (norm(g) > OptTol) & (k < kmax) )
    % ... (modified) Newton direction
    MinEigHess = min( eig(H) );
    if (MinEigHess < MinEig)
        H = H + (-MinEigHess + MinEig)*eye(n);
    end
    s =-H*g;
    % ... Armijo line search
    [alpha, j, fail] = ArmijoLineSearch( x, s, g, fun );
    x = x + alpha*s;
    diff_x = alpha*s;
    % Evaluate function and its gradient at new x
    [f,g_new] = feval(fun,x);
    y = g_new - g;
    rho = 1/(y'*diff_x);
    H = (eye(2)-rho*diff_x*y')*H*(eye(2)-rho*y*diff_x') + rho*(diff_x)*diff_x'; %BFGS update formula to compute Hessian inverse
    
    % ... exist, if line search failed (info <> 0)
    if (fail ~= 0)
        break;
    end; % if
    
    % ... update iterate and function value, gradient
    k = k+1;
    g = g_new;
    
    % ... print progress of method
    if (PrntAlg == 1)
        plot([x(1),x(1)+alpha*s(1)], [x(2),x(2)+alpha*s(2)]);
        figure(1); clf; qq = f + g(1).*(xx-x(1)) + g(2).*(yy-x(2)) ...
            + H(1,1).*(xx-x(1)).^2./2 + H(2,2).*(yy-x(2)).^2./2 + H(2,1).*(xx-x(1)).*(yy-x(2));
        contour(xx,yy,ff,'LineWidth',2); hold on;
        contour(xx,yy,qq,':','LineWidth',2);
        plot(x(1),x(2),'r*');
        title('Progress of Quasi-Newton method on Powell Function');
         %%       title('Progress of Quasi-Newton method on Rosenbrock Function');
        pause;
    end % if
    
    % ... print progress
    fprintf(format1, k, f, norm(g), MinEigHess, alpha, fail);
end; % while

if (fail ~= 0)
    fprintf(format1, k+1, f, norm(g), MinEigHess, alpha, fail);
end; % if

% ... display the result
fail   = fail
x      = x
fmin   = f
fstart = fstart