Interpolação 1D ... de grades 2D

Question

Eu sei que isso pode ser bastante confuso, então, por favor, me avise se esta explicação precisar de alguma edição

Digamos que eu tenha dados de entrada neste formato:

para uma dada pressão p_0-> grade 2x2 de temperaturas ( T_0) que se referem a este valor de pressão

para uma dada pressão p_1-> grade 2x2 de temperaturas ( T_1) que se referem a este valor de pressão

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

Agora, recebo uma grade 2x2 de novos valores de pressão

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

e eu gostaria de obter uma grade 2x2 de valores de temperatura interpolados, usando os dados de entrada.

A forma como estou fazendo isso é para cada ponto da grade, eu construo uma função de interpolação e, em seguida, a utilizo para obter o novo valor de temperatura interpolada para esse ponto da grade:

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

Como é fácil de adivinhar, isso é bastante lento para matrizes grandes e parece ser totalmente contra a maneira entediante de fazer as coisas.

EDIT: Estou usando interp1dtambém porque permite extrapolação, que é uma opção que gostaria de manter.

Answer 1

Você pode apenas calcular a interpolação sozinho. Aqui, presumo que você tenha mais de dois Tvalores e isso pnão é necessariamente espaçado de maneira uniforme. Além disso, o código assume que você tem vários p_targetvalores, mas obviamente funciona para apenas um.

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

EDITAR: Se você quiser extrapolação linear, simplesmente não corte os alphavalores:

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

Answer 2

Eu adicionei alguns parâmetros para as dimensões; por sua escolha de n_x = n_y = n_p = 2, as dependências não eram tão claras.

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

Do que uma palavra para sua modelagem. Se entendi bem o que você quer temperature_xy = fun_xy(pressure), uma função separada para cada coordenada em sua grade espacial. Outra opção pode ser incluir os componentes espaciais em uma função combinada temperature_xy = fun(pressure, x, y). Para a segunda abordagem, dê uma olhada em scipy.interpolate.griddata .

Você pode reorganizar a primeira abordagem para fazê-la funcionar interp2d(). Para isso, a primeira dimensão é a pressão x=pressuree a segunda dimensão representa as dimensões espaciais combinadas y=product(x, y). Para fazer com que isso se comporte como n_x * n_yinterpolações independentes dos valores de pressão, eu apenas uso os mesmos valores fictícios 0, 1, 2 ... para os componentes espaciais tanto ao criar a interpolação quanto ao avaliá-la. Como a avaliação de interp2d()normalmente funciona apenas em coordenadas de grade, usei o método fornecido pelo usuário 6655984 para avaliar a função apenas em um conjunto específico de pontos.

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

Com essas configurações, obtenho uma melhoria de tempo de quase 10x. Mas se você usar valores ainda maiores, como n_x=n_y=1000o uso de memória desta abordagem interp2d personalizada aumenta muito e sua abordagem iterativa vence.

# 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

Com esse conhecimento, você poderia fazer um loop sobre uma grande 1000x1000grade e processar 100x100peças sequencialmente, então você terminaria em cerca de 3 segundos em vez de 30 segundos.

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)