code: def _evaluate(self, indices, norm_distances, out=None):
          
          # slice for broadcasting over trailing dimensions in self.values
          vslice = (slice(None),) + (None,) * (self.values.ndim - len(indices))
  
          if out is None:
              out_shape = out_shape_from_meshgrid(norm_distances)
              out_dtype = self.values.dtype
              out = np.zeros(out_shape, dtype=out_dtype)
          else:
              out[:] = 0.0
  
          # Weights and indices (per axis)
          low_weights, high_weights, edge_indices = _create_weight_edge_lists(
              indices, norm_distances, self.schemes, self.nn_variants)
  
          # Iterate over all possible combinations of [i, i+1] for each
          # axis, resulting in a loop of length 2**ndim
          for lo_hi, edge in zip(product(*([['l', 'h']] * len(indices))),
                                 product(*edge_indices)):
              weight = 1.0
              # TODO: determine best summation order from array strides
              for lh, w_lo, w_hi in zip(lo_hi, low_weights, high_weights):
  
                  # We don't multiply in-place to exploit the cheap operations
                  # in the beginning: sizes grow gradually as following:
                  # (n, 1, 1, ...) -> (n, m, 1, ...) -> ...
                  # Hence, it is faster to build up the weight array instead
                  # of doing full-size operations from the beginning.
                  if lh == 'l':
                      weight = weight * w_lo
                  else:
                      weight = weight * w_hi
              out += np.asarray(self.values[edge]) * weight[vslice]
          return np.array(out, copy=False, ndmin=1)
docs: Evaluate linear interpolation.
  
          Modified for in-place evaluation and treatment of out-of-bounds
          points by implicitly assuming 0 at the next node.