Ускорение scipy griddata для множественных интерполяций между двумя нерегулярными сетками

У меня есть несколько значений, которые определены на одной и той же нерегулярной сетке (x, y, z) которую я хочу интерполировать на новую сетку (x1, y1, z1) . f(x, y, z), g(x, y, z), h(x, y, z) и я хочу рассчитать f(x1, y1, z1), g(x1, y1, z1), h(x1, y1, z1) .

На данный момент я делаю это с помощью scipy.interpolate.griddata и он работает хорошо. Однако, поскольку я должен выполнять каждую интерполяцию по отдельности и есть много точек, она довольно медленная, с большим дублированием в расчете (то есть, поиск ближайших точек, настройка сетки и т. Д.).

Есть ли способ ускорить расчет и уменьшить дублированные вычисления? т.е. что-то вдоль линий определения двух сеток, а затем изменение значений для интерполяции?

3 Solutions collect form web for “Ускорение scipy griddata для множественных интерполяций между двумя нерегулярными сетками”

При каждом вызове scipy.interpolate.griddata происходит несколько вещей:

  1. Во-первых, вызов sp.spatial.qhull.Dealunay сделан для триангуляции нерегулярных координат сетки.
  2. Затем для каждой точки новой сетки выполняется поиск триангуляции, чтобы найти треугольник (на самом деле, в котором симплекс, который в вашем трехмерном случае будет находиться в тетраэдре).
  3. Вычисляются барицентрические координаты каждой новой точки сетки относительно вершин окружающего симплекса.
  4. Интерполированные значения вычисляются для этой точки сетки, используя барицентрические координаты и значения функции в вершинах охватывающего симплекса.

Первые три шага идентичны для всех ваших интерполяций, поэтому, если бы вы могли сохранить для каждой новой точки сетки индексы вершин прилагаемого симплекса и веса для интерполяции, вы бы минимизировали количество вычислений на много. К сожалению, это не так просто сделать с доступной функциональностью, хотя это действительно возможно:

 import scipy.interpolate as spint import scipy.spatial.qhull as qhull import itertools def interp_weights(xyz, uvw): tri = qhull.Delaunay(xyz) simplex = tri.find_simplex(uvw) vertices = np.take(tri.simplices, simplex, axis=0) temp = np.take(tri.transform, simplex, axis=0) delta = uvw - temp[:, d] bary = np.einsum('njk,nk->nj', temp[:, :d, :], delta) return vertices, np.hstack((bary, 1 - bary.sum(axis=1, keepdims=True))) def interpolate(values, vtx, wts): return np.einsum('nj,nj->n', np.take(values, vtx), wts) 

Функция interp_weights выполняет вычисления для первых трех шагов, перечисленных выше. Затем interpolate функции использует эти исчисляемые значения, чтобы сделать шаг 4 очень быстро:

 m, n, d = 3.5e4, 3e3, 3 # make sure no new grid point is extrapolated bounding_cube = np.array(list(itertools.product([0, 1], repeat=d))) xyz = np.vstack((bounding_cube, np.random.rand(m - len(bounding_cube), d))) f = np.random.rand(m) g = np.random.rand(m) uvw = np.random.rand(n, d) In [2]: vtx, wts = interp_weights(xyz, uvw) In [3]: np.allclose(interpolate(f, vtx, wts), spint.griddata(xyz, f, uvw)) Out[3]: True In [4]: %timeit spint.griddata(xyz, f, uvw) 1 loops, best of 3: 2.81 s per loop In [5]: %timeit interp_weights(xyz, uvw) 1 loops, best of 3: 2.79 s per loop In [6]: %timeit interpolate(f, vtx, wts) 10000 loops, best of 3: 66.4 us per loop In [7]: %timeit interpolate(g, vtx, wts) 10000 loops, best of 3: 67 us per loop 

Итак, во-первых, он делает то же самое, что и griddata , что хорошо. Во-вторых, настройка интерполяции, то есть вычисление vtx и wts принимает примерно то же самое, что и вызов griddata . Но в-третьих, вы можете теперь интерполировать для разных значений в одной и той же сетке практически мгновенно.

Единственное, что делает griddata , которое здесь не предусмотрено, это присвоение fill_value точкам, которые нужно экстраполировать. Вы можете сделать это, проверив точки, для которых хотя бы один из весов отрицательный, например:

 def interpolate(values, vtx, wts, fill_value=np.nan): ret = np.einsum('nj,nj->n', np.take(values, vtx), wts) ret[np.any(wts < 0, axis=1)] = fill_value return ret 

Огромное спасибо Хайме за его решение (даже если я не понимаю, как делается барицентрическое вычисление …)

Здесь вы найдете пример, адаптированный из его случая в 2D:

 import scipy.interpolate as spint import scipy.spatial.qhull as qhull import numpy as np def interp_weights(xy, uv,d=2): tri = qhull.Delaunay(xy) simplex = tri.find_simplex(uv) vertices = np.take(tri.simplices, simplex, axis=0) temp = np.take(tri.transform, simplex, axis=0) delta = uv - temp[:, d] bary = np.einsum('njk,nk->nj', temp[:, :d, :], delta) return vertices, np.hstack((bary, 1 - bary.sum(axis=1, keepdims=True))) def interpolate(values, vtx, wts): return np.einsum('nj,nj->n', np.take(values, vtx), wts) m, n = 101,201 mi, ni = 1001,2001 [Y,X]=np.meshgrid(np.linspace(0,1,n),np.linspace(0,2,m)) [Yi,Xi]=np.meshgrid(np.linspace(0,1,ni),np.linspace(0,2,mi)) xy=np.zeros([X.shape[0]*X.shape[1],2]) xy[:,0]=Y.flatten() xy[:,1]=X.flatten() uv=np.zeros([Xi.shape[0]*Xi.shape[1],2]) uv[:,0]=Yi.flatten() uv[:,1]=Xi.flatten() values=np.cos(2*X)*np.cos(2*Y) #Computed once and for all ! vtx, wts = interp_weights(xy, uv) valuesi=interpolate(values.flatten(), vtx, wts) valuesi=valuesi.reshape(Xi.shape[0],Xi.shape[1]) print "interpolation error: ",np.mean(valuesi-np.cos(2*Xi)*np.cos(2*Yi)) print "interpolation uncertainty: ",np.std(valuesi-np.cos(2*Xi)*np.cos(2*Yi)) 

Можно применить преобразование изображения, такое как картирование с ускорением udge

Вы не можете использовать то же определение функции, что и новые координаты, будут меняться на каждой итерации, но вы можете вычислить триангуляцию Один раз для всех.

 import scipy.interpolate as spint import scipy.spatial.qhull as qhull import numpy as np import time # Definition of the fast interpolation process. May be the Tirangulation process can be removed !! def interp_tri(xy): tri = qhull.Delaunay(xy) return tri def interpolate(values, tri,uv,d=2): simplex = tri.find_simplex(uv) vertices = np.take(tri.simplices, simplex, axis=0) temp = np.take(tri.transform, simplex, axis=0) delta = uv- temp[:, d] bary = np.einsum('njk,nk->nj', temp[:, :d, :], delta) return np.einsum('nj,nj->n', np.take(values, vertices), np.hstack((bary, 1.0 - bary.sum(axis=1, keepdims=True)))) m, n = 101,201 mi, ni = 101,201 [Y,X]=np.meshgrid(np.linspace(0,1,n),np.linspace(0,2,m)) [Yi,Xi]=np.meshgrid(np.linspace(0,1,ni),np.linspace(0,2,mi)) xy=np.zeros([X.shape[0]*X.shape[1],2]) xy[:,1]=Y.flatten() xy[:,0]=X.flatten() uv=np.zeros([Xi.shape[0]*Xi.shape[1],2]) # creation of a displacement field uv[:,1]=0.5*Yi.flatten()+0.4 uv[:,0]=1.5*Xi.flatten()-0.7 values=np.zeros_like(X) values[50:70,90:150]=100. #Computed once and for all ! tri = interp_tri(xy) t0=time.time() for i in range(0,100): values_interp_Qhull=interpolate(values.flatten(),tri,uv,2).reshape(Xi.shape[0],Xi.shape[1]) t_q=(time.time()-t0)/100 t0=time.time() values_interp_griddata=spint.griddata(xy,values.flatten(),uv,fill_value=0).reshape(values.shape[0],values.shape[1]) t_g=time.time()-t0 print "Speed-up:", t_g/t_q print "Mean error: ",(values_interp_Qhull-values_interp_griddata).mean() print "Standard deviation: ",(values_interp_Qhull-values_interp_griddata).std() 

На моем ноутбуке ускорение составляет от 20 до 40 раз!

Надежда, которая может помочь кому-то

Вы можете попробовать использовать Pandas , поскольку он обеспечивает высокопроизводительные структуры данных.

Это правда, что метод интерполяции является оберткой scipy-интерполяции, но, возможно, с улучшенными структурами вы получаете лучшую скорость.

 import pandas as pd; wp = pd.Panel(randn(2, 5, 4)); wp.interpolate(); 

interpolate() заполняет значения NaN в наборе данных панели с использованием разных методов . Надеюсь, что это быстрее, чем Scipy.

Если это не сработает , есть один способ повысить производительность (вместо использования параллельной версии вашего кода): используйте Cython и реализуйте небольшую процедуру на C для использования внутри вашего кода Python. Здесь у вас есть пример об этом.

  • ошибка неупорядоченных типов при импорте sklearn
  • 2D Gaussian Fit для интенсивностей в определенных координатах в Python
  • Эффективный метод расчета плотности нерегулярно разнесенных точек
  • Определение логнормального распределения с использованием Scipy vs Matlab
  • Запуск scipy.integrate.ode в многопроцессорном пуле приводит к огромному результату
  • эффективный способ получить максимум каждой строки для большой разреженной матрицы
  • Самый простой способ решения математических уравнений в Python
  • Поиск длины кубического B-сплайна
  • scipy.interpolate.griddata: вырезать z-значение и получить область внутри
  • Установка scipy для python 2.7
  • Чтение файла и построение CDF в Python
  • Python - лучший язык программирования в мире.