code: def interpn(
      points,
      values,
      xi,
      method="linear",
      bounds_error=True,
      fill_value=np.nan,
  ):
      """
      Multidimensional interpolation on regular or rectilinear grids.
  
      Strictly speaking, not all regular grids are supported - this function
      works on *rectilinear* grids, that is, a rectangular grid with even or
      uneven spacing.
  
      Parameters
      ----------
      points : tuple of ndarray of float, with shapes (m1, ), ..., (mn, )
          The points defining the regular grid in n dimensions. The points in
          each dimension (i.e. every elements of the points tuple) must be
          strictly ascending or descending.
  
      values : array_like, shape (m1, ..., mn, ...)
          The data on the regular grid in n dimensions. Complex data is
          accepted.
  
          .. deprecated:: 1.13.0
              Complex data is deprecated with ``method="pchip"`` and will raise an
              error in SciPy 1.15.0. This is because ``PchipInterpolator`` only
              works with real values. If you are trying to use the real components of
              the passed array, use ``np.real`` on ``values``.
  
      xi : ndarray of shape (..., ndim)
          The coordinates to sample the gridded data at
  
      method : str, optional
          The method of interpolation to perform. Supported are "linear",
          "nearest", "slinear", "cubic", "quintic", "pchip", and "splinef2d".
          "splinef2d" is only supported for 2-dimensional data.
  
      bounds_error : bool, optional
          If True, when interpolated values are requested outside of the
          domain of the input data, a ValueError is raised.
          If False, then `fill_value` is used.
  
      fill_value : number, optional
          If provided, the value to use for points outside of the
          interpolation domain. If None, values outside
          the domain are extrapolated.  Extrapolation is not supported by method
          "splinef2d".
  
      Returns
      -------
      values_x : ndarray, shape xi.shape[:-1] + values.shape[ndim:]
          Interpolated values at `xi`. See notes for behaviour when
          ``xi.ndim == 1``.
  
      See Also
      --------
      NearestNDInterpolator : Nearest neighbor interpolation on unstructured
                              data in N dimensions
      LinearNDInterpolator : Piecewise linear interpolant on unstructured data
                             in N dimensions
      RegularGridInterpolator : interpolation on a regular or rectilinear grid
                                in arbitrary dimensions (`interpn` wraps this
                                class).
      RectBivariateSpline : Bivariate spline approximation over a rectangular mesh
      scipy.ndimage.map_coordinates : interpolation on grids with equal spacing
                                      (suitable for e.g., N-D image resampling)
  
      Notes
      -----
  
      .. versionadded:: 0.14
  
      In the case that ``xi.ndim == 1`` a new axis is inserted into
      the 0 position of the returned array, values_x, so its shape is
      instead ``(1,) + values.shape[ndim:]``.
  
      If the input data is such that input dimensions have incommensurate
      units and differ by many orders of magnitude, the interpolant may have
      numerical artifacts. Consider rescaling the data before interpolation.
  
      Examples
      --------
      Evaluate a simple example function on the points of a regular 3-D grid:
  
      >>> import numpy as np
      >>> from scipy.interpolate import interpn
      >>> def value_func_3d(x, y, z):
      ...     return 2 * x + 3 * y - z
      >>> x = np.linspace(0, 4, 5)
      >>> y = np.linspace(0, 5, 6)
      >>> z = np.linspace(0, 6, 7)
      >>> points = (x, y, z)
      >>> values = value_func_3d(*np.meshgrid(*points, indexing='ij'))
  
      Evaluate the interpolating function at a point
  
      >>> point = np.array([2.21, 3.12, 1.15])
      >>> print(interpn(points, values, point))
      [12.63]
  
      """
      # sanity check 'method' kwarg
      if method not in [
          "linear",
          "nearest",
          "cubic",
          "quintic",
          "pchip",
          "splinef2d",
          "slinear",
      ]:
          raise ValueError(
              "interpn only understands the methods 'linear', "
              "'nearest', 'slinear', 'cubic', 'quintic', 'pchip', "
              f"and 'splinef2d'. You provided {method}."
          )
  
      if not hasattr(values, "ndim"):
          values = np.asarray(values)
  
      ndim = values.ndim
      if ndim > 2 and method == "splinef2d":
          raise ValueError(
              "The method splinef2d can only be used for "
              "2-dimensional input data"
          )
      if (
          not bounds_error
          and fill_value is None
          and method == "splinef2d"
      ):
          raise ValueError(
              "The method splinef2d does not support extrapolation."
          )
  
      # sanity check consistency of input dimensions
      if len(points) > ndim:
          raise ValueError(
              "There are %d point arrays, but values has %d "
              "dimensions" % (len(points), ndim)
          )
      if len(points) != ndim and method == "splinef2d":
          raise ValueError(
              "The method splinef2d can only be used for "
              "scalar data with one point per coordinate"
          )
  
      grid, descending_dimensions = _check_points(points)
      _check_dimensionality(grid, values)
  
      # sanity check requested xi
      xi = _ndim_coords_from_arrays(xi, ndim=len(grid))
      if xi.shape[-1] != len(grid):
          raise ValueError(
              "The requested sample points xi have dimension "
              "%d, but this RegularGridInterpolator has "
              "dimension %d" % (xi.shape[-1], len(grid))
          )
  
      if bounds_error:
          for i, p in enumerate(xi.T):
              if not np.logical_and(
                  np.all(grid[i][0] <= p), np.all(p <= grid[i][-1])
              ):
                  raise ValueError(
                      "One of the requested xi is out of bounds "
                      "in dimension %d" % i
                  )
  
      # perform interpolation
      if method in [
          "linear",
          "nearest",
          "slinear",
          "cubic",
          "quintic",
          "pchip",
      ]:
          interp = RegularGridInterpolator(
              points,
              values,
              method=method,
              bounds_error=bounds_error,
              fill_value=fill_value,
          )
          return interp(xi)
      elif method == "splinef2d":
          xi_shape = xi.shape
          xi = xi.reshape(-1, xi.shape[-1])
  
          # RectBivariateSpline doesn't support fill_value; we need to wrap here
          idx_valid = np.all(
              (
                  grid[0][0] <= xi[:, 0],
                  xi[:, 0] <= grid[0][-1],
                  grid[1][0] <= xi[:, 1],
                  xi[:, 1] <= grid[1][-1],
              ),
              axis=0,
          )
          result = np.empty_like(xi[:, 0])
  
          # make a copy of values for RectBivariateSpline
          interp = RectBivariateSpline(points[0], points[1], values[:])
          result[idx_valid] = interp.ev(
              xi[idx_valid, 0], xi[idx_valid, 1]
          )
          result[np.logical_not(idx_valid)] = fill_value
  
          return result.reshape(xi_shape[:-1])
  
imports: import itertools
  import warnings
  import numpy
  from .interpnd import _ndim_coords_from_arrays
  from ._cubic import PchipInterpolator
  from ._rgi_cython import evaluate_linear_2d, find_indices
  from ._bsplines import make_interp_spline
  from ._fitpack2 import RectBivariateSpline
label_desc: Use image node_id 3 for calling a global function with example usage: interpn(points, values, xi, method, bounds_error, fill_value) and returns: interp, result