Zastrzeżenie: piszę ten post głównie z myślą o składni i ogólnym zachowaniu. Nie jestem zaznajomiony z aspektem pamięci i procesora opisanych metod i kieruję tę odpowiedź do tych, którzy mają dość małe zestawy danych, tak że jakość interpolacji może być głównym aspektem do rozważenia. Zdaję sobie sprawę, że pracując z bardzo dużymi zbiorami danych, skuteczniejsze metody (a mianowicie griddata
i Rbf
) mogą nie być wykonalne.
Porównam trzy rodzaje wielowymiarowych metod interpolacji ( interp2d
/ splajny griddata
i Rbf
). Poddam je dwóm rodzajom zadań interpolacyjnych i dwóm rodzajom funkcji bazowych (punkty, z których mają być interpolowane). Konkretne przykłady zademonstrują dwuwymiarową interpolację, ale realne metody można zastosować w dowolnych wymiarach. Każda metoda zapewnia różne rodzaje interpolacji; we wszystkich przypadkach użyję interpolacji sześciennej (lub czegoś zbliżonego do 1 ). Ważne jest, aby pamiętać, że za każdym razem, gdy używasz interpolacji, wprowadzasz odchylenie w porównaniu z surowymi danymi, a określone metody mają wpływ na artefakty, które uzyskasz. Zawsze bądź tego świadomy i interpoluj odpowiedzialnie.
Będą dwa zadania interpolacji
- upsampling (dane wejściowe na siatce prostokątnej, dane wyjściowe na siatce gęstszej)
- interpolacja rozproszonych danych do zwykłej siatki
Dwie funkcje (w domenie [x,y] in [-1,1]x[-1,1]
) będą
- płynna i przyjazna funkcja
cos(pi*x)*sin(pi*y)
:; zakres w[-1, 1]
- funkcja zła (a zwłaszcza nieciągła):
x*y/(x^2+y^2)
o wartości 0,5 w pobliżu źródła; zakres w[-0.5, 0.5]
Oto jak wyglądają:
Najpierw pokażę, jak zachowują się te trzy metody w tych czterech testach, a następnie szczegółowo opiszę składnię wszystkich trzech. Jeśli wiesz, czego powinieneś oczekiwać od metody, możesz nie chcieć tracić czasu na naukę jej składni (patrzenie na siebie interp2d
).
Dane testowe
Dla jasności, oto kod, za pomocą którego wygenerowałem dane wejściowe. Chociaż w tym konkretnym przypadku jestem oczywiście świadomy funkcji leżącej u podstaw danych, użyję tego tylko do wygenerowania danych wejściowych dla metod interpolacji. Używam numpy dla wygody (i głównie do generowania danych), ale samo scipy też by wystarczyło.
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)
Płynna funkcja i upsampling
Zacznijmy od najłatwiejszego zadania. Oto jak działa upsampling z siatki kształtu [6,7]
do jednej z [20,21]
płynnych funkcji testu:
Mimo że jest to proste zadanie, istnieją już subtelne różnice między wynikami. Na pierwszy rzut oka wszystkie trzy wyniki są rozsądne. Warto zwrócić uwagę na dwie cechy, opierając się na naszej wcześniejszej wiedzy na temat podstawowej funkcji: środkowy przypadek griddata
najbardziej zniekształca dane. Zwróć uwagę na y==-1
granicę wykresu (najbliżej x
etykiety): funkcja powinna być ściśle zerowa (ponieważ y==-1
jest linią węzłową dla funkcji gładkiej), ale tak nie jest griddata
. Zwróć także uwagę na x==-1
granicę wykresów (z tyłu, po lewej): podstawowa funkcja ma lokalne maksimum (implikujące zerowy gradient w pobliżu granicy) przy [-1, -0.5]
, ale griddata
wynik pokazuje wyraźnie niezerowy gradient w tym regionie. Efekt jest subtelny, ale mimo wszystko jest to błąd. (WiernośćRbf
jest jeszcze lepszy z domyślnym wyborem funkcji radialnych, dubbingowanych multiquadratic
).
Zła funkcja i upsampling
Nieco trudniejszym zadaniem jest wykonanie upsamplingu na naszej złej funkcji:
Między tymi trzema metodami zaczynają się pojawiać wyraźne różnice. Patrząc na wykresy powierzchni, na wyjściu z interp2d
(zwróć uwagę na dwa garby po prawej stronie wykreślonej powierzchni) widać wyraźne fałszywe ekstrema . Choć griddata
i Rbf
wydaje się produkować podobne wyniki na pierwszy rzut oka, ten ostatni wydaje się produkować głębsze minimum blisko [0.4, -0.4]
, że jest nieobecny z funkcją podstawową.
Jest jednak jeden kluczowy aspekt, który Rbf
jest o wiele lepszy: szanuje symetrię podstawowej funkcji (co jest oczywiście możliwe dzięki symetrii siatki próbki). Wyjście z griddata
łamie symetrię punktów próbkowania, co jest już słabo widoczne w gładkiej obudowie.
Gładka funkcja i rozproszone dane
Najczęściej chce się dokonać interpolacji na danych rozproszonych. Z tego powodu oczekuję, że te testy będą ważniejsze. Jak pokazano powyżej, punkty próbkowania wybrano pseudojednorodnie w interesującej nas dziedzinie. W realistycznych scenariuszach możesz mieć dodatkowy szum przy każdym pomiarze i powinieneś rozważyć, czy na początku ma sens interpolacja surowych danych.
Wyjście dla funkcji płynnej:
Teraz jest już trochę horroru. Przyciąłem dane wyjściowe od interp2d
do między [-1, 1]
wyłącznie na potrzeby kreślenia, aby zachować przynajmniej minimalną ilość informacji. Oczywiste jest, że chociaż niektóre z podstawowych kształtów są obecne, istnieją ogromne hałaśliwe obszary, w których metoda całkowicie się psuje. Drugi przypadek griddata
całkiem ładnie odtwarza kształt, ale zwróć uwagę na białe obszary na granicy działki konturowej. Wynika to z faktu, że griddata
działa tylko wewnątrz wypukłego kadłuba punktów danych wejściowych (innymi słowy, nie dokonuje żadnej ekstrapolacji ). Zachowałem domyślną wartość NaN dla punktów wyjściowych leżących poza wypukłym kadłubem. 2 Biorąc pod uwagę te cechy, Rbf
wydaje się działać najlepiej.
Zła funkcja i rozproszone dane
I moment, na który wszyscy czekaliśmy:
To żadna niespodzianka, że się interp2d
poddaje. W rzeczywistości podczas rozmowy interp2d
powinieneś spodziewać się przyjaznych RuntimeWarning
narzekań na niemożność zbudowania splajnu. Jeśli chodzi o pozostałe dwie metody, Rbf
wydaje się , że daje najlepsze wyniki, nawet w pobliżu granic dziedziny, w której wynik jest ekstrapolowany.
Pozwólcie więc, że powiem kilka słów na temat tych trzech metod, w kolejności malejącej preferencji (tak, aby najgorsze było najmniej prawdopodobne).
scipy.interpolate.Rbf
Rbf
Klasa oznacza „radialnych funkcji bazowych”. Szczerze mówiąc, nigdy nie rozważałem tego podejścia, dopóki nie zacząłem szukać informacji na temat tego postu, ale jestem prawie pewien, że będę ich używać w przyszłości.
Podobnie jak w przypadku metod opartych na splajnach (patrz dalej), użycie odbywa się w dwóch krokach: najpierw tworzy się wywoływalną Rbf
instancję klasy na podstawie danych wejściowych, a następnie wywołuje ten obiekt dla danej siatki wyjściowej, aby uzyskać interpolowany wynik. Przykład z testu płynnego upsamplingu:
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
Zauważ, że zarówno punkty wejściowe, jak i wyjściowe były w tym przypadku tablicami 2d, a dane wyjściowe z_dense_smooth_rbf
mają taki sam kształt jak x_dense
i y_dense
bez żadnego wysiłku. Należy również pamiętać, że Rbf
obsługuje dowolne wymiary interpolacji.
Więc, scipy.interpolate.Rbf
- generuje grzeczne dane wyjściowe nawet dla szalonych danych wejściowych
- obsługuje interpolację w wyższych wymiarach
- ekstrapoluje poza wypukłą powłokę punktów wejściowych (oczywiście ekstrapolacja jest zawsze ryzykowna i generalnie nie należy w ogóle na niej polegać)
- tworzy interpolator jako pierwszy krok, więc ocena go w różnych punktach wyjściowych wymaga mniej dodatkowego wysiłku
- może mieć punkty wyjściowe o dowolnym kształcie (w przeciwieństwie do ograniczania do prostokątnych siatek, patrz dalej)
- skłonne do zachowania symetrii danych wejściowych
- obsługuje wiele rodzajów funkcji radialnych dla słowa kluczowego
function
: multiquadric
, inverse
, gaussian
, linear
, cubic
, quintic
, thin_plate
i zdefiniowany przez użytkownika dowolna
scipy.interpolate.griddata
Mój poprzedni faworyt griddata
to ogólny konik roboczy do interpolacji w dowolnych wymiarach. Nie wykonuje ekstrapolacji poza ustaleniem jednej zadanej wartości dla punktów znajdujących się poza wypukłym kadłubem punktów węzłowych, ale ponieważ ekstrapolacja jest bardzo zmienną i niebezpieczną rzeczą, niekoniecznie jest to oszustwo. Przykład użycia:
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
Zwróć uwagę na nieco niezdarną składnię. Punkty wejściowe muszą być określone w tablicy kształtu [N, D]
w D
wymiarach. W tym celu musimy najpierw spłaszczyć nasze tablice współrzędnych 2D (używając ravel
), a następnie połączyć tablice i transponować wynik. Można to zrobić na wiele sposobów, ale wszystkie wydają się być nieporęczne. Dane wejściowe z
również muszą zostać spłaszczone. Mamy trochę więcej swobody, jeśli chodzi o punkty wyjściowe: z jakiegoś powodu można je również określić jako krotkę wielowymiarowych tablic. Zwróć uwagę, że help
of griddata
jest mylący, ponieważ sugeruje, że to samo dotyczy punktów wejściowych (przynajmniej dla wersji 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.
W skrócie, scipy.interpolate.griddata
- generuje grzeczne dane wyjściowe nawet dla szalonych danych wejściowych
- obsługuje interpolację w wyższych wymiarach
- nie dokonuje ekstrapolacji, można ustawić pojedynczą wartość dla wyjścia poza wypukłą otoczką punktów wejściowych (patrz
fill_value
)
- oblicza interpolowane wartości w jednym wywołaniu, więc sondowanie wielu zestawów punktów wyjściowych zaczyna się od zera
- może mieć punkty wyjściowe o dowolnym kształcie
- obsługuje interpolację najbliższego sąsiada i interpolację liniową w dowolnych wymiarach sześciennych w 1d i 2d. Stosuje się odpowiednio interpolację najbliższego sąsiada i interpolację liniową
NearestNDInterpolator
oraz LinearNDInterpolator
pod maską. Interpolacja sześcienna 1d wykorzystuje splajn, interpolacja sześcienna 2d wykorzystuje CloughTocher2DInterpolator
do konstruowania ciągłego różniczkowalnego interpolatora odcinkowo-sześciennego.
- może naruszyć symetrię danych wejściowych
scipy.interpolate.interp2d
/scipy.interpolate.bisplrep
Jedynym powodem, dla którego omawiam interp2d
i jego krewnych, jest to, że ma zwodniczą nazwę i ludzie prawdopodobnie będą próbować go używać. Alert spoilera: nie używaj go (od wersji Scipy 0.17.0). Jest już bardziej wyjątkowy niż poprzednie przedmioty, ponieważ jest specjalnie używany do interpolacji dwuwymiarowej, ale podejrzewam, że jest to zdecydowanie najczęstszy przypadek interpolacji wielowymiarowej.
Jeśli chodzi o składnię, interp2d
jest podobna do Rbf
tego, że najpierw wymaga skonstruowania instancji interpolacji, którą można wywołać w celu podania rzeczywistych wartości interpolowanych. Jest jednak pewien haczyk: punkty wyjściowe muszą znajdować się na prostokątnej siatce, więc dane wejściowe do wywołania interpolatora muszą być wektorami 1d, które obejmują siatkę wyjściową, jakby z 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
Jednym z najczęstszych błędów podczas używania interp2d
jest umieszczanie pełnych siatek 2d w wywołaniu interpolacyjnym, co prowadzi do gwałtownego zużycia pamięci i miejmy nadzieję, do pośpiechu MemoryError
.
Największym problemem interp2d
jest to, że często nie działa. Aby to zrozumieć, musimy zajrzeć pod maskę. Okazuje się, że interp2d
jest to opakowanie dla funkcji niższego poziomu bisplrep
+ bisplev
, które z kolei są opakowaniami dla procedur FITPACK (napisanych w języku Fortran). Równoważne wywołanie z poprzedniego przykładu byłoby
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
A teraz sprawa interp2d
: (w wersji scipy 0.17.0) jest miły komentarzinterpolate/interpolate.py
dla 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)
i rzeczywiście w interpolate/fitpack.py
, bisplrep
jest trochę konfiguracji i ostatecznie
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)
I to wszystko. Podstawowe procedury interp2d
nie są tak naprawdę przeznaczone do wykonywania interpolacji. Mogą wystarczyć do wystarczająco dobrze wychowanych danych, ale w realnych okolicznościach prawdopodobnie będziesz chciał użyć czegoś innego.
Podsumowując, interpolate.interp2d
- może prowadzić do artefaktów nawet w przypadku dobrze temperowanych danych
- jest przeznaczony specjalnie do problemów z dwiema zmiennymi (chociaż istnieje ograniczenie
interpn
dla punktów wejściowych zdefiniowanych w siatce)
- dokonuje ekstrapolacji
- tworzy interpolator jako pierwszy krok, więc ocena go w różnych punktach wyjściowych wymaga mniej dodatkowego wysiłku
- może generować wyjście tylko przez prostokątną siatkę, dla wyjścia rozproszonego musiałbyś wywołać interpolator w pętli
- obsługuje interpolację liniową, sześcienną i kwintyczną
- może naruszyć symetrię danych wejściowych
1 Jestem prawie pewien, że cubic
i linear
rodzaj funkcji bazowych Rbf
nie odpowiadają dokładnie innym interpolatorom o tej samej nazwie.
2 Te NaN są również powodem, dla którego wykres powierzchni wydaje się tak dziwny: matplotlib historycznie miał trudności z wykreślaniem złożonych obiektów 3D z odpowiednimi informacjami o głębokości. Wartości NaN w danych dezorientują mechanizm renderujący, więc części powierzchni, które powinny znajdować się z tyłu, są wykreślane z przodu. Jest to problem dotyczący wizualizacji, a nie interpolacji.