Может ли scipy RegularGridInterpolator возвращать значения и градиенты одним вызовом?

python scipy interpolation gradient linear-interpolation

611 просмотра

1 ответ

Я использую scipy.interpolate.RegularGridInterpolatorс method='linear'. Получить интерполированные значения просто (см. Пример на https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html ). Что такое хороший способ получить градиенты вместе с интерполированными значениями?

Одна возможность состоит в том, чтобы вызывать интерполятор несколько раз и вычислять градиент «вручную» с использованием конечных разностей. Это кажется расточительным, учитывая, что каждый вызов интерполятора, вероятно, уже вычисляет градиент под капотом. Это правильно? Если так, как я могу изменить RegularGridInterpolator, чтобы он возвращал как значение интерполированной функции, так и ее градиент?

Чтобы было ясно, меня не интересует «истинный» градиент функции, которую я интерполирую - просто градиент линейного приближения, например, градиент my_interpolating_functionв примере на https://docs.scipy.org /doc/scipy-0.16.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html .


Вот пример. У меня есть функция f, я строю линейный интерполятор f_interp, и меня интересует градиент f_interp(в отличие от градиента f). Я могу вычислить это, используя конечные различия, но есть ли лучший способ? Я предполагаю, RegularGridInterpolatorчто уже вычисляет градиент под капотом - и делает это быстро. Как я могу изменить его, чтобы он возвращал градиент вместе с интерполированными значениями?

import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate as interp

def f(x, y, z):
    return 0.01*x**2 + 0.05*x**3 + 5*x*y + 3*x*y*z + 0.1*x*(y**2)*z + 9*x*z**2 - 2*y

x_grid = np.linspace(0.0, 10.0, 20)
y_grid = np.linspace(-10.0, 10.0, 20)
z_grid = np.linspace(0.0, 20.0, 20)

mesh = np.meshgrid(x_grid, y_grid, z_grid, indexing="ij")

f_on_grid = f(mesh[0], mesh[1], mesh[2])
assert np.isclose(f_on_grid[0, 0, 0], f(x_grid[0], y_grid[0], z_grid[0]))  # Sanity check

grid = (x_grid, y_grid, z_grid)
f_interp = interp.RegularGridInterpolator(grid, f_on_grid, method="linear",
                                          bounds_error=False, fill_value=None)

dense_x = np.linspace(0.0, 20.0, 400)
plt.plot(dense_x, [f_interp([x, 1.0, 2.0])[0] for x in dense_x], label="y=1.0, z=2.0")
plt.plot(dense_x, [f_interp([x, 1.0, 4.0])[0] for x in dense_x], label="y=1.0, z=4.0")
plt.legend()
plt.show()

f_interp([0.05, 1.0, 2.0])  # Linearly interpolated value, distinct from f(0.05, 1.0, 2.0)

## Suppose I want to compute both f_interp and its gradient at point_of_interest
point_of_interest = np.array([0.23, 1.67, 5.88])
f_interp(point_of_interest)  # Function value -- how can I get f_interp to also return gradient?

## First gradient calculation using np.gradient and a mesh around point_of_interest +- delta
delta = 0.10
delta_mesh = np.meshgrid(*([-delta, 0.0, delta], ) * 3, indexing="ij")
delta_mesh_long = np.column_stack((delta_mesh[0].flatten(),
                                   delta_mesh[1].flatten(),
                                   delta_mesh[2].flatten()))
assert delta_mesh_long.shape[1] == 3
point_plus_delta_mesh = delta_mesh_long + point_of_interest.reshape((1, 3))
values_for_gradient = f_interp(point_plus_delta_mesh).reshape(delta_mesh[0].shape)
gradient = [x[1, 1, 1] for x in np.gradient(values_for_gradient, delta)]
gradient  # Roughly [353.1, 3.8, 25.2]

## Second gradient calculation using finite differences, should give similar result
gradient = np.zeros((3, ))
for idx in [0, 1, 2]:
    point_right = np.copy(point_of_interest)
    point_right[idx] += delta
    point_left = np.copy(point_of_interest)
    point_left[idx] -= delta
    gradient[idx] = (f_interp(point_right)[0] - f_interp(point_left)[0]) / (2*delta)

gradient  # Roughly [353.1, 3.8, 25.2]

Вот изображение f и f_interp. Меня интересует градиент f_interp (сплошные линии):

f и f_interp

Автор: Adrian Источник Размещён: 08.11.2019 11:32

Ответы (1)


1 плюс

Решение

Нет.

Вот что scipy.interpolate.RegularGridInterpolatorделает под капотом:

class CartesianGrid(object):
    """
    Linear Multivariate Cartesian Grid interpolation in arbitrary dimensions
    This is a regular grid with equal spacing.
    """
    def __init__(self, limits, values):
        self.values = values
        self.limits = limits

    def __call__(self, *coords):
        # transform coords into pixel values
        coords = numpy.asarray(coords)
        coords = [(c - lo) * (n - 1) / (hi - lo) for (lo, hi), c, n in zip(self.limits, coords, self.values.shape)]

        return scipy.ndimage.map_coordinates(self.values, coords, 
            cval=numpy.nan, order=1)

https://github.com/JohannesBuchner/regulargrid/blob/master/regulargrid/cartesiangrid.py

Используется scipy.ndimage.map_coordinatesдля линейной интерполяции. coordsсодержит местоположение в пиксельных координатах. Вы должны иметь возможность использовать эти веса, а также нижнее и верхнее значения в каждом измерении, чтобы выяснить, насколько крутой возрастает интерполяция.

Однако градиент также зависит от значений угловых точек.

Вы можете найти математику здесь: https://en.wikipedia.org/wiki/Trilinear_interpolation

Автор: j13r Размещён: 28.11.2016 02:54
Вопросы из категории :
32x32