1D интерполяция… 2D сеток

Sep 01 2020

Я знаю, что это может сбивать с толку, поэтому дайте мне знать, если это объяснение требует некоторого редактирования.

Допустим, у меня есть входные данные в таком формате:

для данного давления p_0-> сетка температур 2x2 ( T_0), относящихся к этому значению давления

для данного давления p_1-> сетка температур 2x2 ( T_1), относящихся к этому значению давления

p_0 = 0
T_0 = np.array([[1, 4], [3, 2]])

p_1 = 1
T_1 = np.array([[1, 6], [4, 4]])

p = np.array([p_0, p_1])
T = np.array([T_0, T_1])

Теперь мне дана сетка 2x2 новых значений давления.

p_target = np.array([[0.1, 0.4], [0.3, 0.2]])

и я хотел бы получить сетку 2x2 интерполированных значений температуры, используя входные данные.

Я делаю это для каждой точки сетки, я создаю функцию интерполяции, а затем использую ее для получения нового интерполированного значения температуры для этой точки сетки:

from scipy.interpolate import interp1d

T_new = np.empty(p_target.shape)

for ix,iy in np.ndindex(p_target.shape):
    f = interp1d(p, T[:,ix,iy])
    T_new[ix,iy] = f(p_target[ix,iy])

T_new

array([[1. , 4.8],
       [3.3, 2.4]])

Как легко догадаться, это довольно медленно для больших массивов, и, похоже, это совершенно не соответствует принципу numpy.

РЕДАКТИРОВАТЬ: Я использую interp1dтакже, потому что он также позволяет экстраполяцию, которую я хотел бы сохранить.

Ответы

2 jdehesa Sep 02 2020 at 18:06

Вы можете просто вычислить интерполяцию самостоятельно. Здесь я предполагаю, что у вас есть более двух Tзначений, и это pне обязательно равномерно. Кроме того, код предполагает, что у вас есть несколько p_targetзначений, но, очевидно, работает только для одного.

import numpy as np

p_0 = 0
T_0 = np.array([[1., 4.], [3., 2.]])
p_1 = 1
T_1 = np.array([[1., 6.], [4., 4.]])
p = np.array([p_0, p_1])
T = np.array([T_0, T_1])
p_target = np.array([[0.1, 0.4], [0.3, 0.2]])
# Assume you may have several of p_target values
p_target = np.expand_dims(p_target, 0)

# Find the base index for each interpolated value (assume p is sorted)
idx_0 = (np.searchsorted(p, p_target) - 1).clip(0, len(p) - 2)
# And the next index
idx_1 = idx_0 + 1
# Get p values for each interpolated value
a = p[idx_0]
b = p[idx_1]
# Compute interpolation factor
alpha = ((p_target - a) / (b - a)).clip(0, 1)
# Get interpolation values
v_0 = np.take_along_axis(T, idx_0, axis=0)
v_1 = np.take_along_axis(T, idx_1, axis=0)
# Compute interpolation
out = (1 - alpha) * v_0 + alpha * v_1
print(out)
# [[[1.  4.8]
#   [3.3 2.4]]]

РЕДАКТИРОВАТЬ: если вам нужна линейная экстраполяция, просто не обрезайте alphaзначения:

alpha = ((p_target - a) / (b - a))

1 scleronomic Sep 01 2020 at 20:22

Я добавил параметры для размеров; по вашему выбору n_x = n_y = n_p = 2, зависимости были не так ясны.

from scipy.interpolate import interp1d, interp2d, dfitpack

n_x = 30
n_y = 40
n_p = 50
T = np.random.random((n_p, n_x, n_y)) * 100
p = np.random.random(n_p)
p[np.argmin(p)] = 0
p[np.argmax(p)] = 1
p_target = np.random.random((n_x, n_y))

T_new = np.empty(p_target.shape)

for ix, iy in np.ndindex(p_target.shape):
    f = interp1d(p, T[:, ix, iy])
    T_new[ix, iy] = f(p_target[ix, iy])

Чем слово для вашего моделирования. Если я правильно понял, что вы хотите temperature_xy = fun_xy(pressure), отдельная функция для каждой координаты в вашей пространственной сетке. Другой вариант - включить пространственные компоненты в комбинированную функцию temperature_xy = fun(pressure, x, y). Для второго подхода посмотрите scipy.interpolate.griddata .

Вы можете изменить первый подход, чтобы он работал interp2d(). Для этого первое измерение - это давление, x=pressureа второе измерение представляет собой комбинированные пространственные измерения y=product(x, y). Чтобы это выглядело как n_x * n_yнезависимая интерполяция значений давления, я просто использую одни и те же фиктивные значения 0, 1, 2 ... для пространственных компонентов как при создании интерполяции, так и при ее оценке. Поскольку оценка interp2d()нормальности работает только с координатами сетки, я использовал метод, предоставленный user6655984, чтобы оценить функцию только по определенному набору точек.

def evaluate_interp2d(f, x, y):
    """https://stackoverflow.com/a/47233198/7570817"""
    return dfitpack.bispeu(f.tck[0], f.tck[1], f.tck[2], f.tck[3], f.tck[4], x, y)[0]

f2 = interp2d(x=p, y=np.arange(n_x*n_y), z=T.reshape(n_p, n_x*n_y).T)

T_new2 = evaluate_interp2d(f=f2, x=p_target.ravel(), y=np.arange(n_x*n_y))
T_new2 = T_new2.reshape(n_x, n_y)

print(np.allclose(T_new, T_new2))
# True

С этими настройками я получаю улучшение почти 10x. Но если вы используете еще большие значения, например, n_x=n_y=1000использование памяти этим настраиваемым подходом interp2d становится слишком большим, и вы выигрываете итеративный подход.

# np=50
#    nx*ny      1e2      1e4      1e5      1e6
# interp1d  0.0056s  0.3420s  3.4133s  33.390s
# interp2d  0.0004s  0.0388s  2.0954s  191.66s

Обладая этими знаниями, вы могли бы перебирать большую 1000x1000сетку и 100x100последовательно обрабатывать детали, и в итоге у вас получилось бы около 3 секунд вместо 30 секунд.