Interpolasi 1D… dari kisi 2D

Question

Saya tahu ini bisa sangat membingungkan, jadi beri tahu saya jika penjelasan ini perlu diedit.

Katakanlah saya memiliki data masukan dalam format ini:

untuk tekanan tertentu p_0-> 2x2 kisi suhu ( T_0) yang mengacu pada nilai tekanan ini

untuk tekanan tertentu p_1-> 2x2 kisi suhu ( T_1) yang mengacu pada nilai tekanan ini

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

Sekarang, saya diberi kisi 2x2 dengan nilai tekanan baru

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

dan saya ingin mendapatkan kisi 2x2 dari nilai suhu yang diinterpolasi, menggunakan data masukan.

Cara saya melakukan ini adalah untuk setiap titik kisi, saya membangun fungsi interpolasi dan kemudian saya menggunakannya untuk mendapatkan nilai suhu interpolasi baru untuk titik kisi tersebut:

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

Karena mudah ditebak, ini cukup lambat untuk array besar, dan tampaknya cukup bertentangan dengan cara numpy dalam melakukan sesuatu.

EDIT: Saya menggunakan interp1djuga karena memungkinkan untuk ekstrapolasi juga, yang merupakan opsi yang ingin saya simpan.

Answer 1

Anda bisa menghitung sendiri interpolasi. Di sini saya berasumsi Anda memiliki lebih dari dua Tnilai dan itu ptidak selalu berjarak sama. Juga, kode mengasumsikan bahwa Anda memiliki beberapa p_targetnilai, tetapi jelas berfungsi hanya untuk satu nilai.

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

EDIT: Jika Anda ingin ekstrapolasi linier, cukup jangan potong alphanilainya:

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

Answer 2

Saya menambahkan beberapa parameter untuk dimensi; dari pilihan Anda n_x = n_y = n_p = 2, dependensinya tidak begitu jelas.

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

Dari pada kata untuk pemodelan Anda. Jika saya mengerti benar yang Anda inginkan temperature_xy = fun_xy(pressure), fungsi terpisah untuk setiap koordinat di kisi spasial Anda. Pilihan lain mungkin memasukkan komponen spasial dalam fungsi gabungan temperature_xy = fun(pressure, x, y). Untuk pendekatan kedua, lihat scipy.interpolate.griddata .

Anda dapat mengatur ulang pendekatan pertama untuk membuatnya berhasil interp2d(). Untuk ini, dimensi pertama adalah tekanan x=pressuredan dimensi kedua mewakili dimensi spasial gabungan y=product(x, y). Untuk membuat ini berperilaku sebagai n_x * n_yinterpolasi independen dari nilai tekanan, saya hanya menggunakan nilai dummy yang sama 0, 1, 2 ... untuk komponen spasial baik saat membuat interpolasi maupun saat mengevaluasinya. Karena evaluasi interp2d()normaly hanya bekerja pada koordinat grid, saya menggunakan metode yang disediakan oleh user6655984 untuk mengevaluasi fungsi hanya pada satu set titik tertentu.

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

Dengan pengaturan tersebut saya mendapatkan peningkatan waktu hampir 10x. Tetapi jika Anda menggunakan nilai yang lebih besar lagi, seperti n_x=n_y=1000penggunaan memori dari pendekatan interp2d kustom ini tumbuh terlalu besar dan pendekatan iteratif Anda menang.

# 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

Dengan pengetahuan ini Anda dapat melakukan loop pada 1000x1000grid besar dan memproses 100x100potongan secara berurutan, kemudian Anda akan berakhir sekitar 3 detik, bukan 30 detik.

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)