code: def _ndppoly_evaluate(c, xs, ks, xp, dx, extrapolate, out):
      """
      Evaluate a piecewise tensor-product polynomial.
  
      Parameters
      ----------
      c : ndarray, shape (k_1*...*k_d, m_1*...*m_d, n)
          Coefficients local polynomials of order `k-1` in
          `m_1`, ..., `m_d` intervals. There are `n` polynomials
          in each interval.
      xs : d-tuple of ndarray of shape (m_d+1,) each
          Breakpoints of polynomials
      ks : ndarray of int, shape (d,)
          Orders of polynomials in each dimension
      xp : ndarray, shape (r, d)
          Points to evaluate the piecewise polynomial at.
      dx : ndarray of int, shape (d,)
          Orders of derivative to evaluate.  The derivative is evaluated
          piecewise and may have discontinuities.
      extrapolate : int, optional
          Whether to extrapolate to out-of-bounds points based on first
          and last intervals, or to return NaNs.
      out : ndarray, shape (r, n)
          Value of each polynomial at each of the input points.
          For points outside the span ``x[0] ... x[-1]``,
          ``nan`` is returned.
          This argument is modified in-place.
      """
      num_samples = xp.shape[0]
      total_xp = xp.size
      ndims = len(xs)
      num_ks = ks.size
  
      # Compose all n-dimensional breakpoints into a single array
      xs_sizes = cupy.asarray([x.size for x in xs], dtype=cupy.int64)
      xs_offsets = cupy.cumsum(xs_sizes)
      xs_offsets = cupy.r_[0, xs_offsets[:-1]]
      xs_complete = cupy.r_[xs]
  
      xs_sizes_m1 = xs_sizes - 1
      xs_strides = cupy.cumprod(xs_sizes_m1[:0:-1])
      xs_strides = cupy.r_[xs_strides[::-1], 1]
  
      # Map each element on the input to their corresponding dimension
      intervals = cupy.empty(xp.shape, dtype=cupy.int64)
      dim_seq = cupy.arange(ndims, dtype=cupy.int64)
      xp_dims = cupy.broadcast_to(
          cupy.expand_dims(dim_seq, 0), (num_samples, ndims)
      )
      xp_dims = xp_dims.copy()
  
      # Compute n-dimensional intervals
      interval_kernel = INTERVAL_MODULE.get_function(
          "find_breakpoint_position_nd"
      )
      interval_kernel(
          ((total_xp + 128 - 1) // 128,),
          (128,),
          (
              xs_complete,
              xp,
              intervals,
              extrapolate,
              total_xp,
              xp_dims,
              xs_sizes,
              xs_offsets,
          ),
      )
  
      # Compute coefficient displacement stride (in elements)
      c_shape = cupy.asarray(c.shape, dtype=cupy.int64)
      c_strides = (
          cupy.asarray(c.strides, dtype=cupy.int64) // c.itemsize
      )
      c2 = cupy.zeros(
          (num_samples * c.shape[0], 1, 1), dtype=_get_dtype(c)
      )
  
      # Compute order strides
      ks_strides = cupy.cumprod(cupy.r_[1, ks])
      ks_strides = ks_strides[:-1]
  
      ppoly_kernel = _get_module_func(PPOLY_MODULE, "eval_ppoly_nd", c)
      ppoly_kernel(
          ((num_samples + 128 - 1) // 128,),
          (128,),
          (
              c,
              xs_complete,
              xp,
              intervals,
              dx,
              ks,
              c2,
              c_shape,
              c_strides,
              xs_strides,
              xs_offsets,
              ks_strides,
              num_samples,
              ndims,
              num_ks,
              out,
          ),
      )
  
imports: import math
  import cupy
  from cupy._core import internal
  from cupy._core._scalar import get_typename
  from cupy_backends.cuda.api import runtime
  from cupyx.scipy import special
  from cupyx.scipy.interpolate._bspline import BSpline, _get_dtype
  import numpy
label_desc: Use image node_id 3 for calling a global function with example usage: _ndppoly_evaluate(c, xs, ks, xp, dx, extrapolate, out) without return types