function [alpha, x, f, grad, fail, nsteps] = linesearch_wolfe(fgname, x0, d, f0, g0, c1, c2);
%  Nocedal-Wright line search, see Chapter 3 of their book
%  includes cubic interpolation to obtain new point
% call : function [alpha, x, f, grad, fail, nsteps] = linesearch_wolfe(fgname, x0, d, f0, g0, c1, c2);
%  written by M. Overton (overton@cs.nyu.edu)
%  input
%   x0:      intial point
%   f0:      function value at x0
%   g0:      directional derivative of function at x0 in direction d
%   d:       search direction (called p in Nocedal and Wright's book) 
%   fgname:  name of function that returns function and gradient
%            it expects as input only x
%            it is invoked by: [f,g] = feval(fgname, x)
%   c1: Wolfe parameter for the condition 
%          f(x0 + alpha d) < = f0 + c1*alpha*g0   (default 1e-4)
%   c2: Wolfe parameter for the (strong) condition 
%          abs((grad f)(x0 + alpha d)'*d) < = -c2*g0
%        these should satisfy 0 < c1 < c2 < 1     (default 0.9)
%        usually, c1 should be small but c2 should not, but set c2 small
%         if an accurate line search is desired
%      
%  output:
%   alpha: steplength satisfying Wolfe conditions
%   x:     x0 + alpha*d
%   f:     f(x0 + alpha d)
%   grad:  (grad f)(x0 + alpha d)
%   fail:  0 for normal termination (Wolfe conditions satisfied)
%          1 if too many steps were taken in lszoom (Wolfe conds not satisfied)
%          -1 if minimizer was not bracketed (too many expansion steps, 
%                                              function may be unbounded below)
%   nsteps: number of steps taken in lszoom
if nargin == 5  % check if the optional Wolfe parameters were passed
    c1 = 1e-4;
    c2 = 0.9;
elseif nargin == 6
    error('must provide neither or both Wolfe parameters')
end
if c1 <= 0 | c1 >= c2 | c2 >= 1
   warning('Wolfe parameters are supposed to satisfy 0 < c1 < c2 < 1')
end
if g0 >= 0
   error('g0 must be negative, not a descent direction')
end
old = 0; fold = f0; gold = g0;
new = 1;
nexpand = 20;
for k = 1:nexpand
   xnew = x0 + new*d;
   [fnew,gradnew] = feval(fgname, xnew);
   gnew = gradnew'*d;
   if fnew > f0 + c1*new*g0 | ((fnew >= fold) & k > 1)   % gone too far
      [alpha, x, f, grad, fail, nsteps] = lszoom(...
         old, new, fold, fnew, gold, gnew, f0, g0, x0, d, fgname, c1, c2);
      return
   end
   if abs(gnew) <= -c2*g0   % quit, conditions are satisfied first time
      alpha = new; x = xnew; f = fnew; grad = gradnew; fail = 0; nsteps = k;
      return
   end
   if gnew >= 0     % gone past minimizer, but new point is better than old
      [alpha, x, f, grad, fail, nsteps] = lszoom(...
         new, old, fnew, fold, gnew, gold, f0, g0, x0, d, fgname, c1, c2);
      return
   end
   old = new;       % minimizer not bracketed yet
   fold = fnew;
   gold = gnew;
   new = 2*new;     
end % loop
fprintf('Minimizer was not bracketed in %d expansion steps, function may be unbounded below', nexpand)
alpha = new; 
x = xnew; 
f = fnew; 
grad = gradnew; 
fail = -1;
nsteps = 0;  % zoom was never called
