Interpolacja 1D… siatek 2D

Question

Wiem, że może to być dość zagmatwane, więc daj mi znać, jeśli to wyjaśnienie wymaga edycji.

Powiedzmy, że mam dane wejściowe w tym formacie:

dla danego ciśnienia p_0-> siatka 2x2 temperatur ( T_0) odnoszących się do tej wartości ciśnienia

dla danego ciśnienia p_1-> siatka 2x2 temperatur ( T_1) odnoszących się do tej wartości ciśnienia

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])

Teraz otrzymuję siatkę 2x2 nowych wartości ciśnienia

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

i chciałbym uzyskać siatkę 2x2 interpolowanych wartości temperatur, korzystając z danych wejściowych.

Sposób, w jaki to robię, dotyczy każdego punktu siatki, buduję funkcję interpolacji, a następnie używam jej, aby uzyskać nową interpolowaną wartość temperatury dla tego punktu siatki:

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]])

Jak łatwo się domyślić, jest to dość powolne w przypadku dużych tablic i wydaje się być dość sprzeczne z odrętwiałym sposobem robienia rzeczy.

EDYCJA: Używam interp1drównież dlatego, że pozwala to również na ekstrapolację, co jest opcją, którą chciałbym zachować.

Answer 1

Możesz po prostu samodzielnie obliczyć interpolację. Tutaj zakładam, że masz więcej niż dwie Twartości i pniekoniecznie są one rozmieszczone w równych odstępach. Ponadto kod zakłada, że ​​masz kilka p_targetwartości, ale oczywiście działa tylko dla jednej.

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]]]

EDYCJA: Jeśli chcesz ekstrapolacji liniowej, po prostu nie przycinaj alphawartości:

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

Answer 2

Dodałem kilka parametrów do wymiarów; z twojego wyboru n_x = n_y = n_p = 2zależności nie były tak jasne.

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])

Niż słowo do twojego modelowania. Jeśli dobrze zrozumiałem, chcesz temperature_xy = fun_xy(pressure), oddzielną funkcję dla każdej współrzędnej na twojej siatce przestrzennej. Inną opcją może być uwzględnienie komponentów przestrzennych w połączonej funkcji temperature_xy = fun(pressure, x, y). W przypadku drugiego podejścia spójrz na scipy.interpolate.griddata .

Możesz zmienić ustawienie pierwszego podejścia, aby działało interp2d(). W tym przypadku pierwszy wymiar to ciśnienie, x=pressurea drugi wymiar to połączone wymiary przestrzenne y=product(x, y). Aby to zachowywało się jak n_x * n_yniezależne interpolacje wartości ciśnienia, po prostu używam tych samych wartości zastępczych 0, 1, 2 ... dla komponentów przestrzennych zarówno podczas tworzenia interpolacji, jak i podczas jej oceny. Ponieważ ocena interp2d()normalnie działa tylko na współrzędnych siatki, użyłem metody dostarczonej przez user6655984, aby ocenić funkcję tylko na określonym zestawie punktów.

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

Przy tych ustawieniach mam prawie czasową poprawę 10x. Ale jeśli używasz jeszcze większych wartości, takich jak n_x=n_y=1000użycie pamięci w tym niestandardowym podejściu interp2d, staje się zbyt duże i wygrywa podejście iteracyjne.

# 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

Mając tę ​​wiedzę, możesz zapętlić dużą 1000x1000siatkę i przetwarzać 100x100elementy sekwencyjnie, a skończy się to około 3 sekund zamiast 30 sekund.

def interpolate2d_flat(p, p_target_flat, T_flat):
    n_p, n_xy = T_flat.shape
    f2 = interp2d(x=p, y=np.arange(n_xy), z=T_flat.T)
    return evaluate_interp2d(f=f2, x=p_target_flat, y=np.arange(n_xy))


n_splits = n_x * n_y // 1000  # So each patch has size n_p*1000, can be changed 

# Flatten and split the spatial dimensions
T_flat_s = np.array_split(T.reshape(n_p, n_x*n_y), n_splits, axis=1)
p_target_flat_s = np.array_split(p_target.ravel(), n_splits, axis=0)

# Loop over the patches
T_new_flat = np.concatenate([interpolate2d_flat(p=p, p_target_flat=ptf, T_flat=Tf)
                             for (ptf, Tf) in zip(p_target_flat_s, T_flat_s)])
T_new2 = T_new_flat.reshape(n_x, n_y)