code: def numerical_grad(
      f,
      inputs,
      grad_outputs,
      eps=1e-3,
      detect_nondifferentiable=False,
      diff_atol=0,
      diff_rtol=1e-2,
      center_outputs=None,
  ):
      """Computes numerical gradient by finite differences.
  
      This function is used to implement gradient check. For usage example, see
      unit tests of :mod:`chainer.functions`.
  
      By default, ``numerical_grad`` computes the gradient to the first order of
      ``eps``.
  
      Args:
          f (callable): Python function with no arguments that runs forward
              computation and returns the result.
          inputs (tuple of arrays): Tuple of arrays that should be treated as
              inputs. Each element of them is slightly modified to realize
              numerical gradient by finite differences.
          grad_outputs (tuple of arrays or scalars): Tuple of arrays or scalars
              that are treated as output gradients.
          eps (float): Epsilon value of finite differences.
          detect_nondifferentiable (bool):
              ``False`` by default.
              If ``True``, ``numerical_grad`` checks whether ``f`` is
              differentiable at ``inputs``.
              It requires evaluation of ``f`` at 5 points instead of 2.
              As a side effect, the accuracy of numerical gradient will be
              increased to the third order of ``eps``.
              If it turns out that ``f`` is non-differentiable at ``input``,
              ``numerical_grad`` raises
              :class:`~chainer.gradient_check.NondifferentiableError`.
          diff_atol (float):
              Absolute tolerance of fitting error of non-differentiable point
              detection.
          diff_rtol (float):
              Tolerance of fitting error of non-differentiable point detection
              relative to the output values of ``f``.
          center_outputs (tuple of arrays or None):
              Only used if ``detect_nondifferentiable`` is ``True``.
              If specified, these arrays are used as the outputs of ``f`` at
              ``inputs``.
              Otherwise, it is calculated.
              It can be used to reduce the computation if these arrays are
              already calculated before calling ``numerical_grad``.
  
      Returns:
          tuple: Numerical gradient arrays corresponding to ``inputs``.
  
      """
      # TODO(niboshi): Deprecate `center_outputs` argument.
      # If dtype of this argument is not float64, often the resolution is
      # insufficient for numerical gradient calculation. We might use it only
      # when its dtype is float64, but it would be better to simply remove it.
      center_outputs = None
  
      assert eps > 0
      assert isinstance(inputs, (tuple, list))
      for x in inputs:
          if x.dtype.kind != "f":
              raise RuntimeError(
                  "The dtype of input arrays must be kind of float"
              )
  
      inputs = tuple(inputs)
      # Cast grad_outputs to float64
      grad_outputs = tuple(
          [
              None
              if g is None
              else numpy.float64(g)
              if numpy.isscalar(g)
              else g.astype(numpy.float64)
              for g in grad_outputs
          ]
      )
  
      if not chainer.is_arrays_compatible(
          [a for a in inputs + grad_outputs if not numpy.isscalar(a)]
      ):
          raise RuntimeError(
              "Do not mix GPU and CPU arrays in `numerical_grad`"
          )
  
      device = backend.get_device_from_array(*(inputs + grad_outputs))
      xp = device.xp
  
      if xp is cuda.cupy:
          numerical_grad_kernel_1 = cuda.reduce(
              "T y1, T y2, U gy, T eps",
              "V gxi",
              "(y1 - y2) * gy",
              "a + b",
              "gxi += a / (eps * 2)",
              "0",
              "numerical_grad_kernel_1",
          )
          numerical_grad_kernel_3 = cuda.reduce(
              "T y1, T y2, T y3, T y4, U gy, T eps",
              "V gxi",
              "(-y1 + 8 * y2 - 8 * y3 + y4) * gy",
              "a + b",
              "gxi += a / (eps * 6)",
              "0",
              "numerical_grad_kernel_3",
          )
  
      if xp is chainerx:
          grads = [
              xp.zeros(x.shape, numpy.float64, device=x.device)
              for x in inputs
          ]
      else:
          grads = [xp.zeros(x.shape, numpy.float64) for x in inputs]
  
      if detect_nondifferentiable:
          if center_outputs is None:
              ys0 = _copy_arrays(f())
          else:
              ys0 = center_outputs
          nout = len(ys0)
          shapes = [y0.shape for y0 in ys0]
          sizes = numpy.array([y0.size for y0 in ys0])
          cumsizes = numpy.cumsum(sizes)
  
      # Evaluate func at a single input
      def eval_func(x, x_ind, delta, orig):
          x[x_ind] = orig + delta
          ys = _copy_arrays(f())
          assert len(ys) == len(grad_outputs)
          assert all(
              [
                  gy is None
                  for y, gy in zip(ys, grad_outputs)
                  if y is None
              ]
          )
          assert all(
              [
                  gy is None
                  or numpy.isscalar(gy)
                  or y.shape == gy.shape
                  for y, gy in zip(ys, grad_outputs)
              ]
          )
          x[x_ind] = orig
          return ys
  
      # An iteration on a single input displacement
      def iterate_single_input(i_in, x, orig_x, x_ind):
          orig = orig_x[x_ind]
          # `yss` holds a list of output arrays for each of 2 or 5 sampling
          # points.
          if detect_nondifferentiable:
              yss = [
                  eval_func(x, x_ind, -eps * 1.0, orig),
                  eval_func(x, x_ind, -eps * 0.5, orig),
                  ys0,
                  eval_func(x, x_ind, +eps * 0.5, orig),
                  eval_func(x, x_ind, +eps * 1.0, orig),
              ]
          else:
              yss = [
                  eval_func(x, x_ind, -eps * 1, orig),
                  eval_func(x, x_ind, +eps * 1, orig),
              ]
  
          assert all(
              [
                  y is None
                  or (
                      y.shape == yss[0][i].shape
                      and y.dtype == yss[0][i].dtype
                  )
                  for ys in yss
                  for i, y in enumerate(ys)
              ]
          )
  
          # If all the outputs are 0-size, skip non-differentiable check.
          if all([y is None or y.size == 0 for y in yss[0]]):
              detect_nondifferentiable_ = False
          else:
              detect_nondifferentiable_ = detect_nondifferentiable
  
          if detect_nondifferentiable_:
              # Detect non-differentiable point by quadratic fitting
  
              # Check for non-finite output.
              # If any single element in the output arrays has different
              # finiteness among sampled points, that means this is a
              # non-differentiable point.
              # If the function consistently generates non-finite values
              # around the point, we do not treat the point as
              # non-differentiable.
              # (Example: x<0 region for the logarithm function)
              any_nonfinite = False
              for i_out in range(nout):
                  isfinites = [xp.isfinite(ys[i_out]) for ys in yss]
                  if any(
                      (isfinites[0] != isfinites[i]).any()
                      for i in range(1, len(yss))
                  ):
                      s = six.StringIO()
                      s.write(
                          "Tried to compute the numeric gradient on a "
                          "non-differentiable point.\n\n"
                      )
                      s.write("i_in: {}\n".format(i_in))
                      s.write("i_out: {}\n".format(i_out))
                      s.write("x: {}\n".format(inputs[i_in]))
                      s.write("index on x: {}\n".format(x_ind))
                      s.write("eps: {}\n".format(eps))
                      s.write("y[x-eps  ]: {}\n".format(yss[0][i_out]))
                      s.write("y[x-eps/2]: {}\n".format(yss[1][i_out]))
                      s.write("y[x      ]: {}\n".format(yss[2][i_out]))
                      s.write("y[x+eps/2]: {}\n".format(yss[3][i_out]))
                      s.write("y[x+eps  ]: {}\n".format(yss[4][i_out]))
                      raise NondifferentiableError(s.getvalue())
  
                  any_nonfinite |= not all((_).all() for _ in isfinites)
  
              if not any_nonfinite:
                  # Stack flattened outputs to make (5, *)-shaped 2D array
                  ystack = xp.vstack(
                      [xp.hstack([y.ravel() for y in ys]) for ys in yss]
                  )
                  assert ystack.ndim == 2 and ystack.shape[0] == len(
                      yss
                  )
                  # Fit to quadratic
                  if xp is not numpy:
                      ystack = _cpu._to_cpu(ystack)
                  polyfit = numpy.polynomial.polynomial.polyfit
                  _, (residuals, _, _, _) = polyfit(
                      range(len(yss)), ystack, deg=2, full=True
                  )
                  if xp is not numpy:
                      residuals = device.send(residuals)
                  residuals = xp.sqrt(residuals / len(yss))
  
                  # Check for error for each output array
                  for i_out in range(nout):
                      size = sizes[i_out]
                      cumsize = cumsizes[i_out]
                      shape = shapes[i_out]
                      # TODO(niboshi): The following two lines could be
                      # rewritten using xp.stack, which is supported in
                      # NumPy>=1.10
                      ymax = xp.concatenate(
                          [ys[i_out][None] for ys in yss]
                      ).max(axis=0)
                      ymin = xp.concatenate(
                          [ys[i_out][None] for ys in yss]
                      ).min(axis=0)
                      # Restore the shape of flattened residual
                      res = residuals[cumsize - size : cumsize]
                      res = res.reshape(shape)
                      det = utils.force_array(
                          diff_atol + diff_rtol * (ymax - ymin) < res
                      )
                      # Constant output = not nondifferentiable
                      det[ymax == ymin] = False
                      if det.any():
                          s = six.StringIO()
                          s.write(
                              "Tried to compute the numeric gradient on a "
                              "non-differentiable point.\n\n"
                          )
                          s.write("i_in: {}\n".format(i_in))
                          s.write("i_out: {}\n".format(i_out))
                          s.write("x: {}\n".format(inputs[i_in]))
                          s.write("index on x: {}\n".format(x_ind))
                          s.write("eps: {}\n".format(eps))
                          s.write("diff_rtol: {}\n".format(diff_rtol))
                          s.write("diff_atol: {}\n".format(diff_atol))
                          s.write("ymax: {}\n".format(ymax))
                          s.write("ymin: {}\n".format(ymin))
                          s.write(
                              "diff_atol + diff_rtol * (ymax-ymin): {}\n".format(
                                  diff_atol + diff_rtol * (ymax - ymin)
                              )
                          )
                          s.write("fitting errors: {}\n".format(res))
                          s.write(
                              "y[x-eps  ]: {}\n".format(yss[0][i_out])
                          )
                          s.write(
                              "y[x-eps/2]: {}\n".format(yss[1][i_out])
                          )
                          s.write(
                              "y[x      ]: {}\n".format(yss[2][i_out])
                          )
                          s.write(
                              "y[x+eps/2]: {}\n".format(yss[3][i_out])
                          )
                          s.write(
                              "y[x+eps  ]: {}\n".format(yss[4][i_out])
                          )
                          raise NondifferentiableError(s.getvalue())
  
          # Calculate numerical gradient
          for i_out, gy in enumerate(grad_outputs):
              if gy is None:
                  continue
              if not numpy.isscalar(gy):
                  gy = gy.astype(numpy.float64, copy=False)
              gpu_ = xp is cuda.cupy and all(
                  isinstance(ys[i_out], cuda.ndarray) for ys in yss
              )
              # If any output sample is None, all others must be.
              assert all(
                  [
                      (yss[0][i_out] is None) == (yss[j][i_out] is None)
                      for j in range(len(yss))
                  ]
              )
              # If outputs samples are None, the part of numeric gradient for
              # this output is considered as zero: skip the accumulation.
              if yss[0][i_out] is None:
                  continue
  
              if len(yss) == 2:  # 1st order
                  y0 = yss[0][i_out]
                  y1 = yss[1][i_out]
                  if gpu_:
                      numerical_grad_kernel_1(
                          y1, y0, xp.asarray(gy), eps, gx[x_ind]
                      )
                  else:
                      dot = ((y1 - y0) * gy).sum()
                      gx[x_ind] = gx[x_ind] + dot / (2 * eps)
              elif len(yss) == 5:  # 3rd order
                  y0 = yss[0][i_out]
                  y1 = yss[1][i_out]
                  y2 = yss[3][i_out]
                  y3 = yss[4][i_out]
                  if gpu_:
                      numerical_grad_kernel_3(
                          y3, y2, y1, y0, gy, eps, gx[x_ind]
                      )
                  else:
                      num = -y3 + 8 * y2 - 8 * y1 + y0
                      dot = (num * gy).sum()
                      gx[x_ind] = gx[x_ind] + dot / (6 * eps)
              else:
                  assert False
  
      # Calculate numeric gradient
      with configuration.using_config("type_check", False):
          for i_in, (x, gx) in enumerate(six.moves.zip(inputs, grads)):
              orig_x = x.copy()  # hold original value
              for x_ind in numpy.ndindex(x.shape):
                  iterate_single_input(i_in, x, orig_x, x_ind)
  
      return [
          g.astype(x.dtype, copy=False)
          for g, x in six.moves.zip(grads, inputs)
      ]
  
imports: from __future__ import absolute_import
  import warnings
  import numpy
  import six
  import chainer
  from chainer import backend
  from chainer.backends import _cpu
  from chainer.backends import cuda
  from chainer import configuration
  from chainer import testing
  from chainer import utils
  from chainer import variable
  import chainerx
label_desc: Use image node_id 5 for calling a global function with example usage: numerical_grad(f, inputs, grad_outputs, eps, detect_nondifferentiable, diff_atol, diff_rtol, center_outputs) and returns: unknown, ys