code: def gradient(x, fx, method=None, approx=True, extrap=False, **kwds):
      """find gradient of ``fx`` at ``x``, with ``fx`` a function ``z = fx(*x)`` or an array ``z``
  
      Args:
        x (ndarray): an array of shape ``(npts, dim)`` or ``(npts,)``
        fx (ndarray): array of shape ``(npts,)`` **or** function ``z = f(*x)``
        method (str, default=None): string name of the kind of gradient method
        approx (bool, default=True): if True, use local approximation method
        extrap (bool, default=False): if True, extrapolate a bounding box
        new (bool, default=False): include the extrapolated points in the output
  
      Returns:
        array of dimensions ``x.shape``, gradient of the points at ``(x,fx)``
  
      Notes:
        - if ``approx`` is True, use ``approx_fprime`` from ``mystic``, which
          uses a local gradient approximation; otherwise use ``gradient`` from
          ``numpy``, which performs a memory-intensive calculation on a grid.
        - if ``scipy`` is not installed, will use ``numpy.interp`` for 1D
          (non-rbf), or ``rbf`` from ``mystic`` otherwise. Default method is
          ``'nearest'`` for 1D, and is ``'linear'`` otherwise. ``method`` can
          be one of ``('rbf', 'linear', 'cubic', 'nearest', 'inverse',
          'gaussian', 'multiquadric', 'quintic', 'thin_plate')``.
        - if ``extrap`` is True, extrapolate using ``interpf`` with
          ``method='thin_plate'`` (or ``'rbf'`` if ``scipy`` is not installed).
          Alternately, any one of ``('rbf', 'linear', 'cubic', 'nearest',
          'inverse', 'gaussian', 'multiquadric', 'quintic', 'thin_plate')`` can
          be used. If ``extrap`` is a cost function ``z = f(x)``, then directly
          use it in the extrapolation. Using extrapolation can reduce the number
          of ``nan``.
      """  # NOTE: uses 'unique' in all cases
      # XXX: nice to have a test for smoothness, worth exploring?
      slc = (
          slice(None, None)
          if kwds.get("new", False)
          else slice(None, len(x))
      )
      import numpy as np
  
      if not hasattr(fx, "__call__"):
          x, fx = extrapolate(x, fx, method=extrap)
          fx = _interpf(x, fx, method=method)
      elif extrap:
          x = extrapolate(x)
      x = np.asarray(x)
      if approx is True:
          fx = _to_objective(fx)  # conform to gradient interface
          if x.ndim == 1:
              return _fprime(x.reshape(x.shape + (1,)), fx).reshape(
                  x.shape
              )[slc]
          return _fprime(x, fx)[slc]
      q = True  # XXX: is q=True better for memory, worse for accuracy?
      if q is True:
          x, i = _unique(x, index=True)
      else:
          i = slice(None)
      gfx = _gradient(
          x, fx(*grid(*_axes(x)))
      )  # XXX: diagonal w/o full grid?
      if type(gfx) is not list:
          gfx.shape = x.shape
          return gfx[i][slc]
      idx = np.diag_indices(*x.shape)
      return np.array([j[idx] for j in gfx]).T[i][slc]
  
imports: from mystic.math import _rbf
label_desc: Use image node_id 24 for calling a global function with example usage: gradient(x, fx, method, approx, extrap) and returns: np, _fprime, gfx, _fprime