function [x,fval,g,nfe,nge,nhe,xs] = newton_tr(func,dfunc,ddfunc,params,x,...
    Delta,Delta_max,eta,tol,trace)
% [x,fval,gval,nfe,nge,nhe,xs] = ...
%	newton_tr(FUNC,DFUNC,DDFUNC,X,DELTA,DELTA_MAX,ETA,TOL,TRACE)
% Performs one step of Newton's method from X.
%   FUNC   -- function name
%   DFUNC  -- gradient function name
%   DDFUNC -- Hessian function name
%   PARAMS -- paramters for objective function
% Returns X and FUNC(PARAMS,X) as fval, and XS, the list of X values
% for each iteration.
% Also returns NFE (# function evaluations), NGE (# gradient evaluations),
% and NHE (# Hessian evaluations)
%
%   DELTA  -- original trust region size
%   DELTA_MAX -- maximum trust region size
%   ETA    -- quality control parameter for taking steps
%   TOL    -- stopping tolerance: stop if ||grad f(x)||_2 < tol
%   TRACE  -- optional; if present and non-zero
%	then print out information about process

nfe = 0;
nge = 0;
nhe = 0;
% Check input data
if Delta > Delta_max | Delta <= 0 | eta < 0 | eta > 0.25
  fprintf('newton_tr: Error in data (eta, Delta, Delta_max)\n')
  fval = feval(func,params,x);
  nfe = nfe + 1;
  return
end

% Initialization
xs = [];
if nargin <= 8
  trace = 0;
end

while 1	% forever do
  if trace ~= 0
    xs = [xs x];
  end
  % Set up model function
  fval = feval(func,params,x);
  if trace ~= 0
    fprintf('newton_tr: function value = %g\n', fval);
  end
  nfe = nfe + 1;
  g = feval(dfunc,params,x);
  nge = nge + 1;
  [m,n] = size(g);
  if m < n
    g = g';
  end
  if trace ~= 0
    fprintf('newton_tr: ||g|| = %g\n', norm(g,2));
  end
  if norm(g,2) < tol
    return
  end
  H = feval(ddfunc,params,x);
  nhe = nhe + 1;
  
  if trace ~= 0
    fprintf('newton_tr: Delta = %g\n', Delta);
  end

  [p,at_bdry] = solve_tr(g,H,Delta,1e-2,max(trace-1,0));

  if trace ~= 0
    if at_bdry
      fprintf('newton_tr: p at trust-region boundary\n');
    else
      fprintf('newton_tr: p interior to trust region\n');
    end
  end

  % Compute "quality" parameter
  rho = (feval(func,params,x+p)-fval)/(g'*p+0.5*p'*H*p);
  if trace ~= 0
    fprintf('Quality parameter rho = %g\n', rho)
  end
  if rho < 0.25
    Delta = 0.25*norm(p,2);
  else
    if rho > 0.75 & at_bdry
      Delta = min(2*Delta,Delta_max);
    end
    if rho > eta
      if trace ~= 0
	fprintf('newton_tr: take step\n');
      end
      x = x+p;
    elseif trace ~= 0
      fprintf('newton_tr: null step\n');
    end % if
  end % if...else
end % while