code: def line_search_BFGS(
      f, xk, pk, gfk, old_fval, args=(), c1=1e-4, alpha0=1
  ):
      """Minimize over alpha, the function f(xk+alpha pk)
  
      Uses the interpolation algorithm (Armiijo backtracking) as suggested by
      Wright and Nocedal in 'Numerical Optimization', 1999, pg. 56-57
  
      :Returns: (alpha, fc, gc)
      """
  
      xk = atleast_1d(xk)
      fc = 0
      phi0 = old_fval  # compute f(xk) -- done in past loop
      phi_a0 = f(*((xk + alpha0 * pk,) + args))
      fc = fc + 1
      derphi0 = numpy.dot(gfk, pk)
  
      if phi_a0 <= phi0 + c1 * alpha0 * derphi0:
          return alpha0, fc, 0, phi_a0
  
      # Otherwise compute the minimizer of a quadratic interpolant:
  
      alpha1 = (
          -(derphi0)
          * alpha0**2
          / 2.0
          / (phi_a0 - phi0 - derphi0 * alpha0)
      )
      phi_a1 = f(*((xk + alpha1 * pk,) + args))
      fc = fc + 1
  
      if phi_a1 <= phi0 + c1 * alpha1 * derphi0:
          return alpha1, fc, 0, phi_a1
  
      # Otherwise loop with cubic interpolation until we find an alpha which
      # satifies the first Wolfe condition (since we are backtracking, we will
      # assume that the value of alpha is not too small and satisfies the second
      # condition.
  
      while 1:  # we are assuming pk is a descent direction
          factor = alpha0**2 * alpha1**2 * (alpha1 - alpha0)
          a = alpha0**2 * (
              phi_a1 - phi0 - derphi0 * alpha1
          ) - alpha1**2 * (phi_a0 - phi0 - derphi0 * alpha0)
          a = a / factor
          b = -(alpha0**3) * (
              phi_a1 - phi0 - derphi0 * alpha1
          ) + alpha1**3 * (phi_a0 - phi0 - derphi0 * alpha0)
          b = b / factor
  
          alpha2 = (-b + numpy.sqrt(abs(b**2 - 3 * a * derphi0))) / (
              3.0 * a
          )
          phi_a2 = f(*((xk + alpha2 * pk,) + args))
          fc = fc + 1
  
          if phi_a2 <= phi0 + c1 * alpha2 * derphi0:
              return alpha2, fc, 0, phi_a2
  
          if (alpha1 - alpha2) > alpha1 / 2.0 or (
              1 - alpha2 / alpha1
          ) < 0.96:
              alpha2 = alpha1 / 2.0
  
          alpha0 = alpha1
          alpha1 = alpha2
          phi_a0 = phi_a1
          phi_a1 = phi_a2
  
imports: import numpy
  from numpy import atleast_1d, eye, mgrid, argmin, zeros, shape, empty, squeeze, isscalar, vectorize, asarray, absolute, sqrt, Inf, asfarray, isinf
  from mystic import linesearch
  from collections.abc import Callable
  import builtins
label_desc: Use image node_id 14 for calling a global function with example usage: line_search_BFGS(f, xk, pk, gfk, old_fval, args, c1, alpha0) and returns: alpha0,  fc,  int,  phi_a0, alpha1,  fc,  int,  phi_a1, alpha2,  fc,  int,  phi_a2