function [x,fval,g,nfe,nge,nhe,xs] = ...
    newton_wc(func,dfunc,ddfunc,params,x,c1,c2,tol,itmax,trace)
% [X,FVAL,GVAL,NFE,NGE,NHE,XS] = ...
%	newton_wc(FUNC,DFUNC,DDFUNC,PARAMS,X,C1,C2,TOL,ITMAX,TRACE)
% Performs Newton's method from X.
% Line search termination criteria are the (strong) Wolfe conditions:
%
%   f(x+alpha*p) <= f(x) + c1*alpha*g'*p
%   |grad f(x+alpha*p)'*g| <= c2|g'*p|
%
% where p is the Newton search direction.
% Recommended values: c1 = 1e-4 and c2 = 0.9
%   FUNC   -- function name
%   DFUNC  -- gradient function name
%   DDFUNC -- Hessian function name
%   TOL    -- tolerance for the gradient
%   ITMAX  -- maximum number of iterations
% Returns X and FUNC(PARAMS,X), # function eval'ns (NFE), # gradient
% eval'ns (NGE), and # Hessian eval'ns (NHE).
% If TRACE is non-zero, then print information about the process

if nargin <= 9
  trace = 0;
end
nge = 0;
nhe = 0;

fval = feval(func,params,x);
nfe = 1;
xs = x;

for k = 1:itmax
  % Find gradient.
  g = feval(dfunc,params,x);
  nge = nge + 1;
  [m,n] = size(g);
  if m < n
    g = g';
  end
  if norm(g) < tol
    return
  end
  H = feval(ddfunc,params,x);
  nhe = nhe + 1;
  if trace ~= 0
    fprintf('newton_wc: Min eigenvalue of H is: %g\n', min(eig(H)));
  end
  p = - H\g;
  slope = g'*p;

  if slope >= 0	% Slope is **up** hill !!! --- H not +ve def. !!!
    p = -g;	% Use steepest descent instead
    slope = g'*p;
    if trace ~= 0
      fprintf('newton_wc: Using steepest descent direction\n');
    end
  end

  alpha = 1;
  if trace > 1
    p
  end

  [alpha,fval,gval,nfe_ws,nge_ws] = ...
      wolfe_search2(func,dfunc,params,x,p,alpha,c1,c2,fval,g, ...
      max(trace-1,0));
  nfe = nfe + nfe_ws;
  nge = nge + nge_ws;
  if trace ~= 0
    fprintf('newton_wc: function value = %g\n', fval);
  end
  x = x + alpha*p;
  if trace ~= 0
    fprintf('newton_wc: Newton step with alpha = %g\n',alpha);
    xs = [xs, x];
  end
end