function [x,f,iterCount,X,Z] = newtonmin(fcn,h,x0)
%NEWTONMIN  Cholesky Strategy in modified Newton method.CMSC 764 HW1.
%   [x, iterCount] = newtonmin(fcn,hessian,x0) finds the minimizer of
%         function f, which is the 1st component of fcn, provided the 
%         gradian  g, which is the 2nd component of fcn, and the 
%         hessian  h, starting with the initial guess x0.
%   newtonmin returns minimizer x, optimal function value f at x, 
%                     the number of iterations iterCount,
%
%   The algorithm [1]:
%     Until x(k) is a good enough solution,
%       Step 1: Find a search direction p(k).          
%       Step 2: Set x(k+1) = x(k) + alpha_k * p(k).
%     
%     Criteria to accept x:
%          (1)(norm(gk)/ (1+abs(fk))) <= EPS,e.g.,10^(-8)
%       or (2) itCount >= ITMAX ,e.g.,1000%   
%       or (3) norm( xkp1 - xk ) / max([typX ,norm(xk)]) <= EPSX,eg,10^(-4)
% 
%     Step 2:line search[2]:
%            use cvsrch from the class website.
%
%   Examples:
%      newtonmin('fg1','Hf1',[1;1;1])
%      newtonmin('fg2','Hf2',[0;0])


% Reference: 
% [1]algorithm:http://www.cs.umd.edu/users/oleary/a607/607unc1hand.pdf 
% [2]cvsrch:http://www.cs.umd.edu/users/oleary/a607/cvsrch.m
% [3]cvsrch parameters:csmr.ca.sandia.gov/~tgkolda/pubs/SIAM-30645.pdf
% 
% Dianne P. O'Leary oleary@cs.umd.edu 
% Geping Liu  geping@math.umd.edu
% 2006/09/21
% $Revision:$  $Date: 2006/09/24  $


% stopping criteria
EPS = 10^(-8);
ITMAX = 1000;

typX = 1;
typF = 1;
EPSX = 10^(-4);


% parameters for cvsrch from [3]
STP = 1.0;
XTOL = 1e-15; 
STPMIN = 1e-15;
STPMAX = 1e15;
MAXFEV = 20;
FTOL = 1e-4;
GTOL = .9;

% initialize outputs 
xk = x0;
iterCount = 0;
X = [];
Z = [];

n = length(x0);
[fk,gk] = feval(fcn,n,xk);
X = [X, xk];
Z = [Z, fk];

relChangeX = 1;
relChangeGrad = norm(gk) / (typF+abs(fk));
while ( relChangeGrad >= EPS ) && (iterCount <= ITMAX )...
        && (relChangeX >= EPSX )
    
    %Step 1:find search direction pk
    Hk = feval(h,xk);
    delta = findDeltaLM(Hk);
    
    pk = findDirction(gk,n,Hk,delta);
            
    %Step 2:line search      
    [xkp1,fkp1,gkp1,stp,info,nfev] = cvsrch(fcn,n,xk,fk,gk,pk,...
                               STP,FTOL,GTOL,XTOL,STPMIN,STPMAX,MAXFEV);                     
    
    %output 
    %debug disp('search direction is ');disp(pk);
    %debug disp('step length is ');disp(stp);
    %debug disp('new estimate is ');disp(xkp1);
    
    %update
    relChangeX = norm( xkp1 - xk ) / max([typX ,norm(xk)]);
    relChangeGrad = norm(gkp1) / (typF+abs(fkp1));
    
    xk = xkp1;
    fk = fkp1;
    gk = gkp1;
    
    iterCount = iterCount + 1;    
    X = [X,xk];    
    Z = [Z,fk];
    
end

x = xk;
f = fk;