class_docstr: Piecewise tensor product polynomial
  
  The value at point ``xp = (x', y', z', ...)`` is evaluated by first
  computing the interval indices `i` such that::
  
      x[0][i[0]] <= x' < x[0][i[0]+1]
      x[1][i[1]] <= y' < x[1][i[1]+1]
      ...
  
  and then computing::
  
      S = sum(c[k0-m0-1,...,kn-mn-1,i[0],...,i[n]]
              * (xp[0] - x[0][i[0]])**m0
              * ...
              * (xp[n] - x[n][i[n]])**mn
              for m0 in range(k[0]+1)
              ...
              for mn in range(k[n]+1))
  
  where ``k[j]`` is the degree of the polynomial in dimension j. This
  representation is the piecewise multivariate power basis.
  
  Parameters
  ----------
  c : ndarray, shape (k0, ..., kn, m0, ..., mn, ...)
      Polynomial coefficients, with polynomial order `kj` and
      `mj+1` intervals for each dimension `j`.
  x : ndim-tuple of ndarrays, shapes (mj+1,)
      Polynomial breakpoints for each dimension. These must be
      sorted in increasing order.
  extrapolate : bool, optional
      Whether to extrapolate to out-of-bounds points based on first
      and last intervals, or to return NaNs. Default: True.
  
  Attributes
  ----------
  x : tuple of ndarrays
      Breakpoints.
  c : ndarray
      Coefficients of the polynomials.
  
  See also
  --------
  PPoly : piecewise polynomials in 1D
  
  Notes
  -----
  High-order polynomials in the power basis can be numerically
  unstable.
class_docstr_tok: ["Piecewise","tensor","product","polynomial","The","value","at","point","``","xp","=","(","x","'",",","y","'",",","z","'",",","...",")","``","is","evaluated","by","first","computing","the","interval","indices","`","i","`","such","that",":",":","x","[","0","]","[","i","[","0","]","]","<","=","x","'","<","x","[","0","]","[","i","[","0","]","+1","]","x","[","1","]","[","i","[","1","]","]","<","=","y","'","<","x","[","1","]","[","i","[","1","]","+1","]","...","and","then","computing",":",":","S","=","sum","(","c","[","k0-m0-1",",","...",",","kn-mn-1",",","i","[","0","]",",","...",",","i","[","n","]","]","*","(","xp","[","0","]","-","x","[","0","]","[","i","[","0","]","]",")","*","*","m0","*","...","*","(","xp","[","n","]","-","x","[","n","]","[","i","[","n","]","]",")","*","*","mn","for","m0","in","range","(","k","[","0","]","+1",")","...","for","mn","in","range","(","k","[","n","]","+1",")",")","where","``","k","[","j","]","``","is","the","degree","of","the","polynomial","in","dimension","j",".","This","representation","is","the","piecewise","multivariate","power","basis",".","Parameters","--","--","--","--","--","c",":","ndarray",",","shape","(","k0",",","...",",","kn",",","m0",",","...",",","mn",",","...",")","Polynomial","coefficients",",","with","polynomial","order","`","kj","`","and","`","mj+1","`","intervals","for","each","dimension","`","j","`",".","x",":","ndim-tuple","of","ndarrays",",","shapes","(","mj+1",",",")","Polynomial","breakpoints","for","each","dimension",".","These","must","be","sorted","in","increasing","order",".","extrapolate",":","bool",",","optional","Whether","to","extrapolate","to","out-of-bounds","points","based","on","first","and","last","intervals",",","or","to","return","NaNs",".","Default",":","True",".","Attributes","--","--","--","--","--","x",":","tuple","of","ndarrays","Breakpoints",".","c",":","ndarray","Coefficients","of","the","polynomials",".","See","also","--","--","--","--","PPoly",":","piecewise","polynomials","in","1D","Notes","--","--","-","High-order","polynomials","in","the","power","basis","can","be","numerically","unstable","."]
docstr: Evaluate the piecewise polynomial or its derivative
  
  Parameters
  ----------
  x : array-like
      Points to evaluate the interpolant at.
  nu : tuple, optional
      Orders of derivatives to evaluate. Each must be non-negative.
  extrapolate : bool, optional
      Whether to extrapolate to out-of-bounds points based on first
      and last intervals, or to return NaNs.
  
  Returns
  -------
  y : array-like
      Interpolated values. Shape is determined by replacing
      the interpolation axis in the original array with the shape of x.
  
  Notes
  -----
  Derivatives are evaluated piecewise for each polynomial
  segment, even if the polynomial is not differentiable at the
  breakpoints. The polynomial intervals are considered half-open,
  ``[a, b)``, except for the last interval which is closed
  ``[a, b]``.
code: def __call__(self, x, nu=None, extrapolate=None):
      """
      Evaluate the piecewise polynomial or its derivative
  
      Parameters
      ----------
      x : array-like
          Points to evaluate the interpolant at.
      nu : tuple, optional
          Orders of derivatives to evaluate. Each must be non-negative.
      extrapolate : bool, optional
          Whether to extrapolate to out-of-bounds points based on first
          and last intervals, or to return NaNs.
  
      Returns
      -------
      y : array-like
          Interpolated values. Shape is determined by replacing
          the interpolation axis in the original array with the shape of x.
  
      Notes
      -----
      Derivatives are evaluated piecewise for each polynomial
      segment, even if the polynomial is not differentiable at the
      breakpoints. The polynomial intervals are considered half-open,
      ``[a, b)``, except for the last interval which is closed
      ``[a, b]``.
      """
      if extrapolate is None:
          extrapolate = self.extrapolate
      else:
          extrapolate = bool(extrapolate)
  
      ndim = len(self.x)
  
      x = _ndim_coords_from_arrays(x)
      x_shape = x.shape
      x = cupy.ascontiguousarray(
          x.reshape(-1, x.shape[-1]), dtype=cupy.float64
      )
  
      if nu is None:
          nu = cupy.zeros((ndim,), dtype=cupy.int64)
      else:
          nu = cupy.asarray(nu, dtype=cupy.int64)
          if nu.ndim != 1 or nu.shape[0] != ndim:
              raise ValueError("invalid number of derivative orders nu")
  
      dim1 = int(np.prod(self.c.shape[:ndim]))
      dim2 = int(np.prod(self.c.shape[ndim : 2 * ndim]))
      dim3 = int(np.prod(self.c.shape[2 * ndim :]))
      ks = cupy.asarray(self.c.shape[:ndim], dtype=cupy.int64)
  
      out = cupy.empty((x.shape[0], dim3), dtype=self.c.dtype)
      self._ensure_c_contiguous()
  
      _ndppoly_evaluate(
          self.c.reshape(dim1, dim2, dim3),
          self.x,
          ks,
          x,
          nu,
          bool(extrapolate),
          out,
      )
  
      return out.reshape(x_shape[:-1] + self.c.shape[2 * ndim :])
  
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 5 for calling the NdPPoly obj's underlying member method code with example usage: obj.__call__(x, nu, extrapolate) and returns: out