Python scipy.interpolate.LinearNDInterpolator() Examples

def _interpolator(points, coord_source, coord_target, method="idw"):

    coord_source_i, coord_source_j = coord_source
    coord_target_i, coord_target_j = coord_target

    # reshape
    trg = np.vstack((coord_source_i.ravel(), coord_source_j.ravel())).T
    src = np.vstack((coord_target_i.ravel(), coord_target_j.ravel())).T

    if method == "nearest":
        interpolator = NearestNDInterpolator(src, points.ravel(),
                                             tree_options={"balanced_tree": False})
        points_interpolated = interpolator(trg)
    elif method == "linear":
        interpolator = LinearNDInterpolator(src, points.ravel(), fill_value=0)
        points_interpolated = interpolator(trg)
    elif method == "idw":
        interpolator = ipol.Idw(src, trg)
        points_interpolated = interpolator(points.ravel())

    # reshape output
    points_interpolated = points_interpolated.reshape(points.shape)

    return points_interpolated.astype(points.dtype) 

def interpolate_look_angles(data):
    log.info('Interpolating look angles from radar coordinates...')
    log.debug('Radar coordinates extent width %d; length %d',
              data.px_width, data.px_length)
    log.debug('Radar coordinates data: length %d - %d; width %d - %d',
              data.radar_coords[1].min(), data.radar_coords[1].max(),
              data.radar_coords[0].min(), data.radar_coords[0].max())
    log.debug('Binned radar coordinate ranges: length %d - %d; width %d - %d',
              num.nanmin(data.bin_radar_i), num.nanmax(data.bin_radar_i),
              num.nanmin(data.bin_radar_j), num.nanmax(data.bin_radar_j))

    width_coords = num.linspace(0, data.px_width, 50)
    len_coords = num.linspace(0, data.px_length, 50)
    coords = num.asarray(num.meshgrid(width_coords, len_coords))\
        .reshape(2, 2500)

    radar_coords = num.vstack(
        [data.bin_radar_j.ravel() - data.radar_coords[0].min(),
         data.bin_radar_i.ravel() - data.radar_coords[1].min()])

    interp = interpolate.LinearNDInterpolator(coords.T, data.look_angles)
    data.bin_look_angles = interp(radar_coords.T).reshape(
        *data.bin_ps_mean_v.shape)
    return interp 

def default_loss(ip):
    """Loss function that combines `deviations` and `areas` of the triangles.

    Works with `~adaptive.Learner2D` only.

    Parameters
    ----------
    ip : `scipy.interpolate.LinearNDInterpolator` instance

    Returns
    -------
    losses : numpy.ndarray
        Loss per triangle in ``ip.tri``.
    """
    dev = np.sum(deviations(ip), axis=0)
    A = areas(ip)
    losses = dev * np.sqrt(A) + 0.3 * A
    return losses 

def impute_atmosphere_grid(self, grid):
        """
        This function imputes the passed atmosphere grid by linear N-D interpolation.
        As grid is passed by reference, it is not necessary to re-assign the table to
        the return value of this function; the return value is provided for convenience
        only, but the grid is changed in place.
        """

        valid_mask = ~np.isnan(grid[...,0])
        coords = np.array(np.nonzero(valid_mask)).T
        values = grid[valid_mask][:,0]
        it = interpolate.LinearNDInterpolator(coords, values, fill_value=0)
        filled = it(list(np.ndindex(grid[...,0].shape))).reshape(grid[...,0].shape)
        filled[filled==0] = np.nan
        grid[...,0] = filled
        return grid 

def areas(ip):
    """Returns the area per triangle of the triangulation inside
    a `LinearNDInterpolator` instance.

    Is useful when defining custom loss functions.

    Parameters
    ----------
    ip : `scipy.interpolate.LinearNDInterpolator` instance

    Returns
    -------
    areas : numpy.ndarray
        The area per triangle in ``ip.tri``.
    """
    p = ip.tri.points[ip.tri.vertices]
    q = p[:, :-1, :] - p[:, -1, None, :]
    areas = abs(q[:, 0, 0] * q[:, 1, 1] - q[:, 0, 1] * q[:, 1, 0]) / 2
    return areas 

def test_remember_last():
    """
    Check we get the same answer when using masked arrays
    """

    # First set up 4 grid points and fill them randomly
    interpolation_points = {
        (0, 0): np.random.rand(2, 2),
        (0, 1): np.random.rand(2, 2),
        (1, 0): np.random.rand(2, 2),
        (1, 1): np.random.rand(2, 2),
    }

    # Create UnstructuredInterpolator and LinearNDInterpolator with these points
    interpolator = UnstructuredInterpolator(interpolation_points, remember_last=True)

    # Create some random coordinates in this space
    random_nums = np.random.rand(2, 2)
    points_mask = ma.masked_array(random_nums, mask=[[True, False], [True, False]])

    # And interpolate...
    interpolated_points = interpolator(random_nums).T[0]
    interpolated_points_mask = interpolator(points_mask).T[0]

    # Check everything agrees to a reasonable precision
    assert np.all(np.abs(interpolated_points - interpolated_points_mask) < 1e-10) 

def lin_interp(shape, xyd):
    from scipy.interpolate import LinearNDInterpolator

    m, n = shape
    ij, d = xyd[:, 1::-1], xyd[:, 2]
    f = LinearNDInterpolator(ij, d, fill_value=0)
    J, I = np.meshgrid(np.arange(n), np.arange(m))
    IJ = np.vstack([I.flatten(), J.flatten()]).T
    disparity = f(IJ).reshape(shape)
    return disparity 

def _ip(self):
        """A `scipy.interpolate.LinearNDInterpolator` instance
        containing the learner's data."""
        # XXX: take our own triangulation into account when generating the _ip
        return interpolate.LinearNDInterpolator(self.points, self.values) 

def _interpolator_combined(self):
        """A `scipy.interpolate.LinearNDInterpolator` instance
        containing the learner's data *and* interpolated data of
        the `pending_points`."""
        if self._ip_combined is None:
            points, values = self._data_combined()
            points = self._scale(points)
            self._ip_combined = interpolate.LinearNDInterpolator(points, values)
        return self._ip_combined 

def interpolator(self, *, scaled=False):
        """A `scipy.interpolate.LinearNDInterpolator` instance
        containing the learner's data.

        Parameters
        ----------
        scaled : bool
            Use True if all points are inside the
            unit-square [(-0.5, 0.5), (-0.5, 0.5)] or False if
            the data points are inside the ``learner.bounds``.

        Returns
        -------
        interpolator : `scipy.interpolate.LinearNDInterpolator`

        Examples
        --------
        >>> xs, ys = [np.linspace(*b, num=100) for b in learner.bounds]
        >>> ip = learner.interpolator()
        >>> zs = ip(xs[:, None], ys[None, :])
        """
        if scaled:
            if self._ip is None:
                points, values = self._data_in_bounds()
                points = self._scale(points)
                self._ip = interpolate.LinearNDInterpolator(points, values)
            return self._ip
        else:
            points, values = self._data_in_bounds()
            return interpolate.LinearNDInterpolator(points, values) 

def uniform_loss(ip):
    """Loss function that samples the domain uniformly.

    Works with `~adaptive.Learner2D` only.

    Parameters
    ----------
    ip : `scipy.interpolate.LinearNDInterpolator` instance

    Returns
    -------
    losses : numpy.ndarray
        Loss per triangle in ``ip.tri``.

    Examples
    --------
    >>> from adaptive.learner.learner2D import uniform_loss
    >>> def f(xy):
    ...     x, y = xy
    ...     return x**2 + y**2
    >>>
    >>> learner = adaptive.Learner2D(
    ...     f,
    ...     bounds=[(-1, -1), (1, 1)],
    ...     loss_per_triangle=uniform_loss,
    ... )
    >>>
    """
    return np.sqrt(areas(ip)) 

def deviations(ip):
    """Returns the deviation of the linear estimate.

    Is useful when defining custom loss functions.

    Parameters
    ----------
    ip : `scipy.interpolate.LinearNDInterpolator` instance

    Returns
    -------
    deviations : list
        The deviation per triangle.
    """
    values = ip.values / (ip.values.ptp(axis=0).max() or 1)
    gradients = interpolate.interpnd.estimate_gradients_2d_global(
        ip.tri, values, tol=1e-6
    )

    p = ip.tri.points[ip.tri.vertices]
    vs = values[ip.tri.vertices]
    gs = gradients[ip.tri.vertices]

    def deviation(p, v, g):
        dev = 0
        for j in range(3):
            vest = v[:, j, None] + (
                (p[:, :, :] - p[:, j, None, :]) * g[:, j, None, :]
            ).sum(axis=-1)
            dev += abs(vest - v).max(axis=1)
        return dev

    n_levels = vs.shape[2]
    devs = [deviation(p, vs[:, :, i], gs[:, :, i]) for i in range(n_levels)]
    return devs 

def __call__(self, vals, fill_value=np.nan):
        """
        Evaluate interpolator for values given at the source points.

        You can interpolate multiple datasets of source values (``vals``) at
        once: the ``vals`` array should have the shape (number of source
        points, number of source datasets). If you want to interpolate only one
        set of source values, ``vals`` can have the shape (number of source
        points, 1) or just (number of source points,) - which is a flat/1-D
        array. The output will have the same number of dimensions as ``vals``,
        i.e. it will be a flat 1-D array in case ``vals`` is a 1-D array.

        Parameters
        ----------
        vals : ndarray of float, shape (numsourcepoints, ...)
            Values at the source points which to interpolate
        fill_value : float
            is needed if linear interpolation fails; defaults to np.nan

        Returns
        -------
        output : ndarray of float with shape (numtargetpoints,...)

        """
        self._check_shape(vals)
        isnan = np.isnan(vals)
        if self.remove_missing & np.count_nonzero(isnan):
            ip = sinterp.LinearNDInterpolator(
                self.src[~isnan, ...], vals[~isnan], fill_value=fill_value
            )
        else:
            ip = sinterp.LinearNDInterpolator(self.src, vals, fill_value=fill_value)
        return ip(self.trg) 

def lin_interp(shape, xyd):
    # taken from https://github.com/hunse/kitti
    m, n = shape
    ij, d = xyd[:, 1::-1], xyd[:, 2]
    f = LinearNDInterpolator(ij, d, fill_value=0)
    # f = NearestNDInterpolator(ij, d)
    J, I = np.meshgrid(np.arange(n), np.arange(m))
    IJ = np.vstack([I.flatten(), J.flatten()]).T
    disparity = f(IJ).reshape(shape)
    return disparity 

def _initialize_distance_interpolator_flat(self, layer):
        self._layer = layer
        self.triangulate_layer_flat(layer=self._layer)

        self._interpolator = [None, None]
        for side in [0, 1]:
            self._interpolator[side] = LinearNDInterpolator(
                self.surf_triang[side], self.triangulation_points[side][:, 2]) 

def test_masked_input():
    """
    Now lets test how well this all works if we pass a masked input
    """

    # First set up 4 grid points and fill them randomly
    interpolation_points = {
        (0, 0): np.random.rand(2, 2),
        (0, 1): np.random.rand(2, 2),
        (1, 0): np.random.rand(2, 2),
        (1, 1): np.random.rand(2, 2),
    }

    # Create UnstructuredInterpolator and LinearNDInterpolator with these points
    interpolator = UnstructuredInterpolator(interpolation_points, remember_last=True)
    linear_nd = LinearNDInterpolator(
        list(interpolation_points.keys()), list(interpolation_points.values())
    )

    # Create some random coordinates in this space
    points = np.random.rand(10, 2)
    # And interpolate...
    interpolator(points)
    interpolated_points = interpolator(points)

    linear_nd_points = linear_nd(points)

    # Check everything agrees to a reasonable precision
    assert np.all(np.abs(interpolated_points - linear_nd_points) < 1e-10) 

def test_linear_nd():
    """
    In its simplest configuration this code should behave exactly the same as the scipy
    LinearNDInterpolator, so lets test that
    """

    # First set up 4 grid points and fill them randomly
    interpolation_points = {
        (0, 0): np.random.rand(2, 2),
        (0, 1): np.random.rand(2, 2),
        (1, 0): np.random.rand(2, 2),
        (1, 1): np.random.rand(2, 2),
    }

    # Create UnstructuredInterpolator and LinearNDInterpolator with these points
    interpolator = UnstructuredInterpolator(interpolation_points)
    linear_nd = LinearNDInterpolator(
        list(interpolation_points.keys()), list(interpolation_points.values())
    )

    # Create some random coordinates in this space
    points = np.random.rand(10, 2)
    # And interpolate...
    interpolated_points = interpolator(points)
    linear_nd_points = linear_nd(points)

    # Check everything agrees to a reasonable precision
    assert np.all(np.abs(interpolated_points - linear_nd_points) < 1e-10) 

def __init__(self, filename, verbose=1):
        """
        Initialisation of class to load templates from a file and create the interpolation
        objects

        Parameters
        ----------
        filename: string
            Location of Template file
        verbose: int
            Verbosity level,
            0 = no logging
            1 = File + interpolation point information
            2 = Detailed description of interpolation points
        """
        self.verbose = verbose
        if self.verbose:
            print("Loading lookup tables from", filename)

        grid, bins, template = self.parse_fits_table(filename)
        x_bins, y_bins = bins

        self.interpolator = interpolate.LinearNDInterpolator(
            grid, template, fill_value=0
        )
        self.nearest_interpolator = interpolate.NearestNDInterpolator(grid, template)

        self.grid_interp = interpolate.RegularGridInterpolator(
            (x_bins, y_bins),
            np.zeros([x_bins.shape[0], y_bins.shape[0]]),
            method="linear",
            bounds_error=False,
            fill_value=0,
        ) 

def build_interpolator(self, dem_proc):
        # Build an interpolator
        gc = dem_proc.elev.grid_coordinates
#       points = np.meshgrid(gc.x_axis, gc.y_axis)
#       points = np.column_stack([pts.ravel() for pts in points])
#       interp = spinterp.NearestNDInterpolator(points, dem_proc.data.ravel())
#       interp = spinterp.LinearNDInterpolator(points, np.ravel(dem_proc.data),
#                                               fill_value=np.nan)
        interp = spinterp.interpolate.RegularGridInterpolator(
            points=(gc.y_axis[::-1], gc.x_axis),
            values=dem_proc.data[::-1, :].astype(float),
            method='nearest', fill_value=np.nan, bounds_error=False)
        return interp 

def make_linear_interpolator_unstructured(field, grid=None, fill_value=np.nan):
	'''Make a linear interpolator for an unstructured grid.

	Parameters
	----------
	field : Field or array_like
		The field to interpolate.
	grid : Grid or None
		The grid of the field. If it is given, the grid of `field` is replaced by this grid.
	fill_value : scalar
		The value to use for points outside of the domain of the input field. Extrapolation is not supported.

	Returns
	-------
	Field generator
		The interpolator as a Field generator. The grid on which this field generator will be evaluated does
		not need to have any structure.
	'''
	if fill_value is None:
		raise ValueError('Extrapolation is not supported for a linear interpolator on an unstructured grid.')

	if grid is None:
		grid = field.grid
	else:
		field = Field(field, grid)

	interp = LinearNDInterpolator(grid.points, field, fill_value)

	def interpolator(evaluated_grid):
		res = interp(grid.points)
		return Field(res, evaluated_grid)

	return interpolator 

def get_greedy(self, v):
        """
        Compute optimal actions taking v as the value function.

        Parameters
        ----------
        v : array_like(float, ndim=1, length=len(π_grid))
            An approximate value function represented as a
            one-dimensional array.

        Returns
        -------
        policy : array_like(float, ndim=1, length=len(π_grid))
            The decision to accept or reject an offer where 1 indicates
            accept and 0 indicates reject

        """
        # == Simplify names == #
        f, g, β, c, q = self.f, self.g, self.β, self.c, self.q

        vf = LinearNDInterpolator(self.grid_points, v)
        N = len(v)
        policy = np.zeros(N, dtype=int)

        for i in range(N):
            w, π = self.grid_points[i, :]
            v1 = w / (1 - β)
            integrand = lambda m: vf(m, q(m, π)) * (π * f(m) +
                                                     (1 - π) * g(m))
            integral, error = fixed_quad(integrand, 0, self.w_max)
            v2 = c + β * integral
            policy[i] = v1 > v2  # Evaluates to 1 or 0

        return policy 

def bellman_operator(self, v):
        """

        The Bellman operator.  Including for comparison. Value function
        iteration is not recommended for this problem.  See the
        reservation wage operator below.

        Parameters
        ----------
        v : array_like(float, ndim=1, length=len(π_grid))
            An approximate value function represented as a
            one-dimensional array.

        Returns
        -------
        new_v : array_like(float, ndim=1, length=len(π_grid))
            The updated value function

        """
        # == Simplify names == #
        f, g, β, c, q = self.f, self.g, self.β, self.c, self.q

        vf = LinearNDInterpolator(self.grid_points, v)
        N = len(v)
        new_v = np.empty(N)

        for i in range(N):
            w, π = self.grid_points[i, :]
            v1 = w / (1 - β)
            integrand = lambda m: vf(m, q(m, π)) * (π * f(m)
                                                     + (1 - π) * g(m))
            integral, error = fixed_quad(integrand, 0, self.w_max)
            v2 = c + β * integral
            new_v[i] = max(v1, v2)

        return new_v 

def lin_interp(shape, xyd):
    # taken from https://github.com/hunse/kitti
    from scipy.interpolate import LinearNDInterpolator
    m, n = shape
    ij, d = xyd[:, 1::-1], xyd[:, 2]
    f = LinearNDInterpolator(ij, d, fill_value=0)
    J, I = np.meshgrid(np.arange(n), np.arange(m))
    IJ = np.vstack([I.flatten(), J.flatten()]).T
    disparity = f(IJ).reshape(shape)
    return disparity 

def define_tc_speed_model(
        normalized_VDI253_model, m1000_curve_factor, idle_engine_speed):
    """
    Define torque converter speed model.

    :param normalized_VDI253_model:
        Normalized VDI253 model function.
    :type normalized_VDI253_model: scipy.interpolate.LinearNDInterpolator

    :param m1000_curve_factor:
        Rescaling factor of m1000 curve [N*m/1e6].
    :type m1000_curve_factor: float

    :param idle_engine_speed:
        Idle engine speed and its standard deviation [RPM].
    :type idle_engine_speed: (float, float)

    :return:
        Torque converter speed model.
    :rtype: callable
    """

    # noinspection PyMissingOrEmptyDocstring,PyUnusedLocal
    def model(times, **kwargs):
        gbs, gbt = kwargs['gear_box_speeds_in'], kwargs['gear_box_torques_in']
        es = normalized_VDI253_model((gbs, gbt / m1000_curve_factor))
        return np.nan_to_num(es - np.maximum(gbs, idle_engine_speed[0]))

    return model 

def test_class_output():
    """
    The final test is to use the more useful functionality of interpolating between the
    outputs of a class member function. I will do this by interpolating between a
    number of numpy regulat grid interpolators and comparing to the output of the
    linear nd interpolator. Again this is a crazy use case, but is a good test.
    """

    x = np.linspace(0, 1, 11)
    # Create a bunch of random numbers to interpolate between
    rand_numbers = np.random.rand(4, 11, 11)

    # Create input for UnstructuredInterpolator
    interpolation_points = {
        (0, 0): RegularGridInterpolator((x, x), rand_numbers[0]),
        (0, 1): RegularGridInterpolator((x, x), rand_numbers[1]),
        (1, 0): RegularGridInterpolator((x, x), rand_numbers[2]),
        (1, 1): RegularGridInterpolator((x, x), rand_numbers[3]),
    }

    # Create some random points to evaluate our interpolators
    pts1 = np.random.rand(1, 2)
    pts2 = np.random.rand(10, 2)

    interpolator = UnstructuredInterpolator(interpolation_points)
    unsort_value = interpolator(pts1, pts2)

    interpolation_points = {
        (0, 0): rand_numbers[0],
        (0, 1): rand_numbers[1],
        (1, 0): rand_numbers[2],
        (1, 1): rand_numbers[3],
    }

    # Perform the same operation by interpolating the values of the full numpy array
    linear_nd = LinearNDInterpolator(
        list(interpolation_points.keys()), list(interpolation_points.values())
    )
    array_out = linear_nd(pts1)
    # Then interpolate on this grid
    reg_interpolator = RegularGridInterpolator((x, x), array_out[0])
    lin_nd_val = reg_interpolator(pts2)

    # Check they give the same answer
    assert np.all(np.abs(unsort_value - lin_nd_val) < 1e-10) 

def define_normalized_VDI253_model(
        m1000_curve_ratios, m1000_curve_norm_torques, idle_engine_speed,
        engine_max_speed, k_factor_curve):
    """
    Defines normalized VDI253 model function.

    :param m1000_curve_ratios:
        Speed ratios of m1000 curve [-].
    :type m1000_curve_ratios: numpy.array

    :param m1000_curve_norm_torques:
        Normalized torques of m1000 curve [-].
    :type m1000_curve_norm_torques: numpy.array

    :param idle_engine_speed:
        Idle engine speed and its standard deviation [RPM].
    :type idle_engine_speed: (float, float)

    :param engine_max_speed:
        Maximum allowed engine speed [RPM].
    :type engine_max_speed: float

    :param k_factor_curve:
        k factor curve.
    :type k_factor_curve: callable

    :return:
        Normalized VDI253 model function.
    :rtype: scipy.interpolate.LinearNDInterpolator
    """
    from scipy.interpolate import interp1d, LinearNDInterpolator
    maximum_ratio, idle = np.max(m1000_curve_ratios), idle_engine_speed[0]
    eng_s = np.linspace(idle, engine_max_speed, 100)
    gb_s = np.linspace(0, engine_max_speed * maximum_ratio, 250)
    x, z = np.meshgrid(gb_s, eng_s)
    r = x / z
    b = r <= maximum_ratio
    x, z, r = x[b], z[b], r[b]
    func = interp1d(m1000_curve_ratios, m1000_curve_norm_torques, kind='cubic')
    tout = func(r) * k_factor_curve(r) * z ** 2
    return LinearNDInterpolator((x, tout), z)


# noinspection PyPep8Naming 

def calibrate_m1000_curve_factor(
        full_load_curve, normalized_VDI253_model, clutch_phases,
        engine_speeds_out_hot, gear_box_speeds_in, gear_box_torques_in,
        clutch_tc_speeds_delta):
    """
    Calibrate the rescaling factor of m1000 curve [N*m/1e6].

    :param full_load_curve:
        Vehicle full load curve.
    :type full_load_curve: function

    :param normalized_VDI253_model:
        Normalized VDI253 model function.
    :type normalized_VDI253_model: scipy.interpolate.LinearNDInterpolator

    :param clutch_phases:
        When the clutch is active [-].
    :type clutch_phases: numpy.array

    :param engine_speeds_out_hot:
        Engine speed at hot condition [RPM].
    :type engine_speeds_out_hot: numpy.array

    :param gear_box_speeds_in:
        Gear box speed vector [RPM].
    :type gear_box_speeds_in: numpy.array

    :param gear_box_torques_in:
        Torque required vector [N*m].
    :type gear_box_torques_in: numpy.array

    :param clutch_tc_speeds_delta:
        Engine speed delta due to the clutch or torque converter [RPM].
    :type clutch_tc_speeds_delta: numpy.array

    :return:
        Rescaling factor of m1000 curve [N*m/1e6].
    :rtype: float
    """
    if clutch_phases.sum() <= 10:
        return sh.NONE
    from co2mpas.utils import mae
    from scipy.optimize import fmin

    # noinspection PyUnresolvedReferences
    es, gbs, gbt, predict, ds = (
        engine_speeds_out_hot[clutch_phases], gear_box_speeds_in[clutch_phases],
        gear_box_torques_in[clutch_phases], normalized_VDI253_model.predict,
        clutch_tc_speeds_delta[clutch_phases]

    )

    def _err(factor):
        e = mae(ds, np.nan_to_num(predict((gbs, gbt / factor)) - es))
        return np.float32(e)

    return fmin(_err, default_m1000_curve_factor(full_load_curve)) 

def minimize_triangle_surface_loss(ip):
    """Loss function that is similar to the distance loss function in the
    `~adaptive.Learner1D`. The loss is the area spanned by the 3D
    vectors of the vertices.

    Works with `~adaptive.Learner2D` only.

    Parameters
    ----------
    ip : `scipy.interpolate.LinearNDInterpolator` instance

    Returns
    -------
    losses : numpy.ndarray
        Loss per triangle in ``ip.tri``.

    Examples
    --------
    >>> from adaptive.learner.learner2D import minimize_triangle_surface_loss
    >>> def f(xy):
    ...     x, y = xy
    ...     return x**2 + y**2
    >>>
    >>> learner = adaptive.Learner2D(f, bounds=[(-1, -1), (1, 1)],
    ...     loss_per_triangle=minimize_triangle_surface_loss)
    >>>
    """
    tri = ip.tri
    points = tri.points[tri.vertices]
    values = ip.values[tri.vertices]
    values = values / (ip.values.ptp(axis=0).max() or 1)

    def _get_vectors(points):
        delta = points - points[:, -1, :][:, None, :]
        vectors = delta[:, :2, :]
        return vectors[:, 0, :], vectors[:, 1, :]

    a_xy, b_xy = _get_vectors(points)
    a_z, b_z = _get_vectors(values)

    a = np.hstack([a_xy, a_z])
    b = np.hstack([b_xy, b_z])

    return np.linalg.norm(np.cross(a, b) / 2, axis=1) 

def information_from_grid(self, theta_grid, fisherinformation_grid, option="smooth", inverse="exact"):
        """
        Loads a grid of coordinates and corresponding Fisher Information, which is then interpolated.
            
        Parameters
        ----------
        theta_grid : ndarray
            List if parameter points `theta` at which the Fisher information matrices `I_ij(theta)`
            is evaluated. Shape (n_gridpoints, n_dimension).
            
        fisherinformation_grid : ndarray
            List if Fisher information matrices `I_ij(theta)`. Shape (n_gridpoints, n_dimension, n_dimension).
            
        option : {"smooth", "linear"}
            Defines if the Fisher Information is interpolated smoothly using the function
            CloughTocher2DInterpolator() or piecewise linear using LinearNDInterpolator().
            Default = 'smooth'.
            
        inverse : {"exact", "interpolate"}
            Defines if the inverse Fisher Information is obtained by either first interpolating
            the Fisher Information and then inverting it ("exact") or by first inverting the grid
            of Fisher Informations and then interpolating the inverse ("interpolate"). Default = 'exact'.
        """

        self.infotype = "grid"
        self.inverse = inverse

        # load from file
        theta_grid = load_and_check(theta_grid)
        fisherinformation_grid = load_and_check(fisherinformation_grid)
        self.dimension = len(fisherinformation_grid[0])

        # Interpolate Information
        if option == "linear":
            self.infofunction = LinearNDInterpolator(points=theta_grid, values=np.array(fisherinformation_grid))
        elif option == "smooth":
            self.infofunction = CloughTocher2DInterpolator(points=theta_grid, values=np.array(fisherinformation_grid))
        else:
            RuntimeError("Option %s unknown", option)

        # Interpolate inverse information
        if self.inverse == "interpolate":
            inv_fisherinformation_grid = np.array([np.linalg.inv(info) for info in fisherinformation_grid])
            if option == "linear":
                self.infofunction_inv = LinearNDInterpolator(points=theta_grid, values=inv_fisherinformation_grid)
            elif option == "smooth":
                self.infofunction_inv = CloughTocher2DInterpolator(points=theta_grid, values=inv_fisherinformation_grid) 

def _make_interp_probe_count_for_dataset_nd_fn(probe_counts):
    """Generate and return a function that interpolates probe count for a dataset.

    This uses a function from scipy's interpolate package to operate on
    an arbitrary number of parameters.

    Args:
        probe_counts: dict giving number of probes for each dataset and
            choice of parameters

    Returns:
        function whose input is a dataset and values for arbitrary
        parameters. The function linearly interpolates the number of
        probes required in that dataset for those parameter values, based
        on the values (which were explicitly calculated) in probe_counts.
    """
    # Reset the memoized dict for a new call (this is useful for unit tests,
    # which may call this function multiple times with different inputs --
    # i.e., values in probe_counts)
    _interp_nd_fn_memoized = {}

    def interp_probe_count_for_dataset(dataset, param_vals):
        """
        Using the given probe counts at particular parameter values, interpolate
        the number of probes for 'dataset' and the parameter values given
        in 'param_vals', where each of these may be floats.
        """
        if dataset in _interp_nd_fn_memoized:
            nd_fn = _interp_nd_fn_memoized[dataset]
        else:
            points = []
            values = []
            for p in probe_counts[dataset].keys():
                points += [p]
                values += [probe_counts[dataset][p]]
            points = np.array(points)
            values = np.array(values)

            nd_fn = interpolate.LinearNDInterpolator(points, values,
                rescale=True)
            _interp_nd_fn_memoized[dataset] = nd_fn

        try:
            return nd_fn(np.array(param_vals))[0]
        except ValueError:
            raise ValueError(param_vals, dataset, probe_counts[dataset])

    return interp_probe_count_for_dataset