Penafian: Saya kebanyakan menulis posting ini dengan pertimbangan sintaksis dan perilaku umum. Saya tidak terbiasa dengan memori dan aspek CPU dari metode yang dijelaskan, dan jawaban ini saya tujukan pada mereka yang memiliki kumpulan data yang cukup kecil, sehingga kualitas interpolasi dapat menjadi aspek utama yang perlu dipertimbangkan. Saya sadar bahwa ketika bekerja dengan kumpulan data yang sangat besar, metode yang berkinerja lebih baik (yaitu griddata
dan Rbf
) mungkin tidak dapat dilakukan.
Saya akan membandingkan tiga jenis metode interpolasi multi-dimensi ( interp2d
/ splines, griddata
dan Rbf
). Saya akan membebankan mereka pada dua jenis tugas interpolasi dan dua jenis fungsi dasar (poin yang akan diinterpolasi). Contoh spesifik akan menunjukkan interpolasi dua dimensi, tetapi metode yang layak dapat diterapkan dalam dimensi yang berubah-ubah. Setiap metode menyediakan berbagai macam interpolasi; dalam semua kasus saya akan menggunakan interpolasi kubik (atau sesuatu yang mendekati 1 ). Penting untuk diperhatikan bahwa setiap kali Anda menggunakan interpolasi, Anda menimbulkan bias dibandingkan dengan data mentah Anda, dan metode spesifik yang digunakan memengaruhi artefak yang akan Anda hadapi. Selalu sadari hal ini, dan lakukan interpolasi secara bertanggung jawab.
Dua tugas interpolasi akan menjadi
- upsampling (data masukan berada pada kisi persegi panjang, data keluaran berada pada kisi yang lebih padat)
- interpolasi data yang tersebar ke kisi biasa
Kedua fungsi (di atas domain [x,y] in [-1,1]x[-1,1]
) akan menjadi
- a halus dan fungsi ramah:
cos(pi*x)*sin(pi*y)
; berkisar dalam[-1, 1]
- fungsi jahat (dan khususnya, tidak berkelanjutan):
x*y/(x^2+y^2)
dengan nilai 0,5 di dekat asalnya; berkisar dalam[-0.5, 0.5]
Begini tampilannya:
Pertama-tama saya akan mendemonstrasikan bagaimana ketiga metode berperilaku di bawah keempat pengujian ini, kemudian saya akan merinci sintaks dari ketiganya. Jika Anda tahu apa yang diharapkan dari suatu metode, Anda mungkin tidak ingin membuang waktu mempelajari sintaksnya (melihat Anda, interp2d
).
Uji data
Demi kesederhanaan, berikut adalah kode yang saya gunakan untuk menghasilkan data masukan. Sementara dalam kasus khusus ini saya jelas mengetahui fungsi yang mendasari data, saya hanya akan menggunakan ini untuk menghasilkan input untuk metode interpolasi. Saya menggunakan numpy untuk kenyamanan (dan sebagian besar untuk menghasilkan data), tetapi scipy saja sudah cukup.
import numpy as np
import scipy.interpolate as interp
# auxiliary function for mesh generation
def gimme_mesh(n):
minval = -1
maxval = 1
# produce an asymmetric shape in order to catch issues with transpositions
return np.meshgrid(np.linspace(minval,maxval,n), np.linspace(minval,maxval,n+1))
# set up underlying test functions, vectorized
def fun_smooth(x, y):
return np.cos(np.pi*x)*np.sin(np.pi*y)
def fun_evil(x, y):
# watch out for singular origin; function has no unique limit there
return np.where(x**2+y**2>1e-10, x*y/(x**2+y**2), 0.5)
# sparse input mesh, 6x7 in shape
N_sparse = 6
x_sparse,y_sparse = gimme_mesh(N_sparse)
z_sparse_smooth = fun_smooth(x_sparse, y_sparse)
z_sparse_evil = fun_evil(x_sparse, y_sparse)
# scattered input points, 10^2 altogether (shape (100,))
N_scattered = 10
x_scattered,y_scattered = np.random.rand(2,N_scattered**2)*2 - 1
z_scattered_smooth = fun_smooth(x_scattered, y_scattered)
z_scattered_evil = fun_evil(x_scattered, y_scattered)
# dense output mesh, 20x21 in shape
N_dense = 20
x_dense,y_dense = gimme_mesh(N_dense)
Fungsi halus dan upsampling
Mari kita mulai dengan tugas termudah. Inilah cara upampling dari bentuk jaring [6,7]
menjadi salah satu hasil [20,21]
untuk fungsi uji yang mulus:
Meskipun ini adalah tugas sederhana, sudah ada perbedaan kecil di antara hasilnya. Sekilas ketiga keluaran ini masuk akal. Ada dua fitur yang perlu diperhatikan, berdasarkan pengetahuan kami sebelumnya tentang fungsi yang mendasarinya: kasus tengah griddata
yang paling banyak mendistorsi data. Perhatikan y==-1
batas plot (paling dekat dengan x
label): fungsinya harus benar-benar nol (karena y==-1
merupakan garis nodal untuk fungsi smooth), namun tidak demikian griddata
. Perhatikan juga x==-1
batas plot (di belakang, ke kiri): fungsi yang mendasarinya memiliki maksimum lokal (menyiratkan gradien nol di dekat batas) di [-1, -0.5]
, namun griddata
keluarannya menunjukkan dengan jelas gradien bukan nol di wilayah ini. Efeknya halus, tapi tetap saja bias. (KesetiaanRbf
bahkan lebih baik dengan pilihan default fungsi radial, dijuluki multiquadratic
.)
Fungsi jahat dan upsampling
Tugas yang sedikit lebih sulit adalah melakukan upampling pada fungsi jahat kita:
Perbedaan yang jelas mulai terlihat di antara ketiga metode tersebut. Melihat plot permukaan, ada ekstrema palsu yang jelas muncul di keluaran dari interp2d
(perhatikan dua punuk di sisi kanan permukaan yang diplot). Meskipun sekilas griddata
dan Rbf
tampaknya menghasilkan hasil yang serupa, yang terakhir tampaknya menghasilkan nilai minimum yang lebih dalam di dekat [0.4, -0.4]
yang tidak ada pada fungsi yang mendasarinya.
Namun, ada satu aspek penting yang Rbf
jauh lebih unggul: ia menghormati kesimetrian fungsi yang mendasarinya (yang tentu saja juga dimungkinkan oleh simetri jaring sampel). Output dari griddata
memecah simetri titik sampel, yang sudah terlihat lemah dalam kasus halus.
Fungsi halus dan data yang tersebar
Paling sering seseorang ingin melakukan interpolasi pada data yang tersebar. Untuk alasan ini saya berharap tes ini menjadi lebih penting. Seperti yang ditunjukkan di atas, titik sampel dipilih secara pseudo-seragam dalam domain yang diminati. Dalam skenario realistis, Anda mungkin memiliki gangguan tambahan dengan setiap pengukuran, dan Anda harus mempertimbangkan apakah masuk akal untuk menginterpolasi data mentah Anda.
Keluaran untuk fungsi mulus:
Sekarang sudah ada sedikit pertunjukan horor yang sedang berlangsung. Saya memotong keluaran dari interp2d
ke antara [-1, 1]
secara eksklusif untuk merencanakan, untuk menjaga setidaknya sejumlah informasi minimal. Jelas bahwa sementara beberapa bentuk yang mendasarinya ada, ada daerah berisik yang sangat besar tempat metode ini benar-benar rusak. Kasus kedua griddata
mereproduksi bentuk dengan cukup baik, tetapi perhatikan daerah putih di perbatasan plot kontur. Hal ini disebabkan oleh fakta bahwa griddata
hanya bekerja di dalam cembung titik data masukan (dengan kata lain, tidak melakukan ekstrapolasi ). Saya menyimpan nilai NaN default untuk titik output yang terletak di luar cembung. 2 Mempertimbangkan fitur-fitur ini, Rbf
tampaknya berkinerja terbaik.
Fungsi jahat dan data yang tersebar
Dan saat yang kita semua tunggu-tunggu:
Tidaklah mengherankan jika interp2d
menyerah. Bahkan, selama panggilan ke interp2d
Anda harus mengharapkan beberapa RuntimeWarning
keluhan ramah tentang ketidakmungkinan spline akan dibangun. Adapun dua metode lainnya, Rbf
tampaknya menghasilkan keluaran terbaik, bahkan dekat perbatasan domain tempat hasil diekstrapolasi.
Jadi izinkan saya menjelaskan beberapa kata tentang ketiga metode tersebut, dalam mengurangi urutan preferensi (sehingga yang terburuk adalah yang paling kecil kemungkinannya untuk dibaca oleh siapa pun).
scipy.interpolate.Rbf
The Rbf
kelas singkatan dari "radial fungsi dasar". Sejujurnya saya tidak pernah mempertimbangkan pendekatan ini sampai saya mulai meneliti untuk posting ini, tetapi saya cukup yakin saya akan menggunakan ini di masa depan.
Sama seperti metode berbasis spline (lihat nanti), penggunaan datang dalam dua langkah: pertama membuat Rbf
instance kelas yang dapat dipanggil berdasarkan data masukan, dan kemudian memanggil objek ini untuk mesh keluaran yang diberikan untuk mendapatkan hasil interpolasi. Contoh dari uji upampling halus:
import scipy.interpolate as interp
zfun_smooth_rbf = interp.Rbf(x_sparse, y_sparse, z_sparse_smooth, function='cubic', smooth=0) # default smooth=0 for interpolation
z_dense_smooth_rbf = zfun_smooth_rbf(x_dense, y_dense) # not really a function, but a callable class instance
Perhatikan bahwa titik masukan dan keluaran adalah array 2d dalam kasus ini, dan keluarannya z_dense_smooth_rbf
memiliki bentuk yang sama x_dense
dan y_dense
tanpa upaya apa pun. Perhatikan juga yang Rbf
mendukung dimensi arbitrer untuk interpolasi.
Begitu, scipy.interpolate.Rbf
- menghasilkan keluaran yang berperilaku baik bahkan untuk data masukan gila
- mendukung interpolasi dalam dimensi yang lebih tinggi
- mengekstrapolasi di luar lambung cembung dari titik-titik masukan (tentu saja ekstrapolasi selalu merupakan pertaruhan, dan Anda biasanya tidak harus bergantung padanya sama sekali)
- membuat interpolator sebagai langkah pertama, jadi mengevaluasinya di berbagai titik keluaran akan mengurangi upaya tambahan
- dapat memiliki titik keluaran dalam bentuk sembarang (sebagai lawan yang dibatasi untuk mata jaring persegi panjang, lihat nanti)
- cenderung mempertahankan simetri data masukan
- mendukung beberapa macam fungsi radial untuk kata kunci
function
: multiquadric
, inverse
, gaussian
, linear
, cubic
, quintic
, thin_plate
dan user-defined sewenang-wenang
scipy.interpolate.griddata
Favorit saya sebelumnya griddata
,, adalah pekerja keras umum untuk interpolasi dalam dimensi sewenang-wenang. Itu tidak melakukan ekstrapolasi di luar pengaturan nilai preset tunggal untuk titik-titik di luar lambung cembung dari titik-titik nodal, tetapi karena ekstrapolasi adalah hal yang sangat berubah-ubah dan berbahaya, ini belum tentu sebuah penipuan. Contoh penggunaan:
z_dense_smooth_griddata = interp.griddata(np.array([x_sparse.ravel(),y_sparse.ravel()]).T,
z_sparse_smooth.ravel(),
(x_dense,y_dense), method='cubic') # default method is linear
Perhatikan sintaks yang sedikit kludgy. Titik masukan harus ditentukan dalam larik bentuk [N, D]
dalam D
dimensi. Untuk ini, pertama-tama kita harus meratakan array koordinat 2d kita (menggunakan ravel
), lalu menggabungkan array dan mengubah urutan hasilnya. Ada beberapa cara untuk melakukan ini, tetapi semuanya tampak besar. z
Data masukan juga harus diratakan. Kami memiliki sedikit lebih banyak kebebasan dalam hal poin keluaran: untuk beberapa alasan ini juga dapat ditentukan sebagai tupel array multidimensi. Perhatikan bahwa help
dari griddata
menyesatkan, karena menunjukkan bahwa hal yang sama juga berlaku untuk titik masukan (setidaknya untuk versi 0.17.0):
griddata(points, values, xi, method='linear', fill_value=nan, rescale=False)
Interpolate unstructured D-dimensional data.
Parameters
----------
points : ndarray of floats, shape (n, D)
Data point coordinates. Can either be an array of
shape (n, D), or a tuple of `ndim` arrays.
values : ndarray of float or complex, shape (n,)
Data values.
xi : ndarray of float, shape (M, D)
Points at which to interpolate data.
Pendeknya, scipy.interpolate.griddata
- menghasilkan keluaran yang berperilaku baik bahkan untuk data masukan gila
- mendukung interpolasi dalam dimensi yang lebih tinggi
- tidak melakukan ekstrapolasi, satu nilai dapat diatur untuk keluaran di luar lambung cembung dari titik masukan (lihat
fill_value
)
- menghitung nilai interpolasi dalam satu panggilan, jadi menyelidiki beberapa set titik keluaran dimulai dari awal
- dapat memiliki titik keluaran dengan bentuk yang berubah-ubah
- mendukung terdekat-tetangga dan interpolasi linier dalam dimensi sewenang-wenang, kubik dalam 1d dan 2d. Tetangga terdekat dan penggunaan interpolasi linier
NearestNDInterpolator
dan di LinearNDInterpolator
bawah kap, masing-masing. Interpolasi kubik 1d menggunakan spline, interpolasi kubik 2d digunakan CloughTocher2DInterpolator
untuk membangun interpolator sepotong-sepotong kubik yang dapat dibedakan secara kontinyu.
- mungkin melanggar simetri data masukan
scipy.interpolate.interp2d
/scipy.interpolate.bisplrep
Satu-satunya alasan saya berdiskusi interp2d
dan kerabatnya adalah karena nama itu menipu, dan orang-orang cenderung mencoba menggunakannya. Peringatan spoiler: jangan gunakan (mulai scipy versi 0.17.0). Ini sudah lebih istimewa daripada subjek sebelumnya karena secara khusus digunakan untuk interpolasi dua dimensi, tapi saya menduga ini adalah kasus paling umum untuk interpolasi multivariat.
Sejauh sintaks berjalan, interp2d
mirip dengan Rbf
yang pertama perlu membangun sebuah contoh interpolasi, yang dapat dipanggil untuk memberikan nilai interpolasi aktual. Namun ada kekurangannya: titik keluaran harus ditempatkan pada jaring persegi panjang, jadi masukan yang masuk ke panggilan ke interpolator harus berupa vektor 1d yang menjangkau kisi keluaran, seolah-olah dari numpy.meshgrid
:
# reminder: x_sparse and y_sparse are of shape [6, 7] from numpy.meshgrid
zfun_smooth_interp2d = interp.interp2d(x_sparse, y_sparse, z_sparse_smooth, kind='cubic') # default kind is 'linear'
# reminder: x_dense and y_dense are of shape [20, 21] from numpy.meshgrid
xvec = x_dense[0,:] # 1d array of unique x values, 20 elements
yvec = y_dense[:,0] # 1d array of unique y values, 21 elements
z_dense_smooth_interp2d = zfun_smooth_interp2d(xvec,yvec) # output is [20, 21]-shaped array
Salah satu kesalahan paling umum saat menggunakan interp2d
adalah memasukkan mesh 2d penuh Anda ke dalam panggilan interpolasi, yang mengarah pada konsumsi memori yang eksplosif, dan mudah-mudahan menjadi terburu-buru MemoryError
.
Sekarang, masalah terbesarnya interp2d
adalah sering tidak berhasil. Untuk memahami ini, kita harus melihat di balik kap mesin. Ternyata itu interp2d
adalah pembungkus untuk fungsi tingkat yang lebih rendah bisplrep
+ bisplev
, yang pada gilirannya menjadi pembungkus untuk rutinitas FITPACK (ditulis dalam Fortran). Panggilan yang setara dengan contoh sebelumnya adalah
kind = 'cubic'
if kind=='linear':
kx=ky=1
elif kind=='cubic':
kx=ky=3
elif kind=='quintic':
kx=ky=5
# bisplrep constructs a spline representation, bisplev evaluates the spline at given points
bisp_smooth = interp.bisplrep(x_sparse.ravel(),y_sparse.ravel(),z_sparse_smooth.ravel(),kx=kx,ky=ky,s=0)
z_dense_smooth_bisplrep = interp.bisplev(xvec,yvec,bisp_smooth).T # note the transpose
Sekarang, inilah hal tentang interp2d
: (dalam scipy versi 0.17.0) ada komentar yanginterpolate/interpolate.py
bagus untuk interp2d
:
if not rectangular_grid:
# TODO: surfit is really not meant for interpolation!
self.tck = fitpack.bisplrep(x, y, z, kx=kx, ky=ky, s=0.0)
dan memang interpolate/fitpack.py
, di bisplrep
sana ada beberapa penyiapan dan akhirnya
tx, ty, c, o = _fitpack._surfit(x, y, z, w, xb, xe, yb, ye, kx, ky,
task, s, eps, tx, ty, nxest, nyest,
wrk, lwrk1, lwrk2)
Dan itu dia. Rutinitas yang mendasarinya interp2d
sebenarnya tidak dimaksudkan untuk melakukan interpolasi. Mereka mungkin cukup untuk data berperilaku cukup baik, tetapi dalam keadaan realistis Anda mungkin ingin menggunakan sesuatu yang lain.
Sekadar menyimpulkan, interpolate.interp2d
- dapat menyebabkan artefak bahkan dengan data yang ditempa dengan baik
- secara khusus untuk masalah bivariat (meskipun ada batasan
interpn
untuk titik masukan yang ditentukan pada kisi)
- melakukan ekstrapolasi
- membuat interpolator sebagai langkah pertama, jadi mengevaluasinya di berbagai titik keluaran akan mengurangi upaya tambahan
- hanya dapat menghasilkan keluaran melalui kisi persegi panjang, untuk keluaran yang tersebar Anda harus memanggil interpolator dalam satu lingkaran
- mendukung interpolasi linier, kubik dan kuintik
- mungkin melanggar simetri data masukan
1 Saya cukup yakin bahwa fungsi dasar cubic
dan linear
jenis Rbf
tidak sama persis dengan interpolator lain dengan nama yang sama.
2 NaN ini juga merupakan alasan mengapa plot permukaan tampak begitu aneh: matplotlib secara historis mengalami kesulitan dalam merencanakan objek 3d yang kompleks dengan informasi kedalaman yang tepat. Nilai NaN dalam data membingungkan perender, sehingga bagian permukaan yang seharusnya berada di belakang diplot berada di depan. Ini adalah masalah dengan visualisasi, dan bukan interpolasi.