polinomial - polinomio de newton en python

Speedup scipy griddata para múltiples interpolaciones entre dos grillas irregulares

Question

Tengo varios valores que están definidos en la misma cuadrícula irregular (x, y, z) que quiero interpolar en una nueva cuadrícula (x1, y1, z1) . es decir, tengo f(x, y, z), g(x, y, z), h(x, y, z) y quiero calcular f(x1, y1, z1), g(x1, y1, z1), h(x1, y1, z1) .

En este momento estoy haciendo esto usando scipy.interpolate.griddata y funciona bien. Sin embargo, debido a que tengo que realizar cada interpolación por separado y hay muchos puntos, es bastante lento, con una gran duplicación en el cálculo (es decir, encontrar qué puntos están más cerca, configurar las cuadrículas, etc.).

¿Hay alguna manera de acelerar el cálculo y reducir los cálculos duplicados? es decir, algo en la línea de definir las dos cuadrículas, luego cambiar los valores para la interpolación?

Answer 1

Hay varias cosas sucediendo cada vez que haces una llamada a scipy.interpolate.griddata :

1. Primero, se realiza una llamada a sp.spatial.qhull.Delaunay para triangular las coordenadas de la cuadrícula irregular.
2. Luego, para cada punto de la nueva cuadrícula, se busca la triangulación para encontrar en qué triángulo (en realidad, en qué simplex, cuál en su caso 3D estará en qué tetraedro) yace.
3. Se calculan las coordenadas baricéntricas de cada nuevo punto de cuadrícula con respecto a los vértices del simplex adjunto.
4. Se calculan valores interpolados para ese punto de la cuadrícula, utilizando las coordenadas baricéntricas y los valores de la función en los vértices del símplex simplificado.

Los primeros tres pasos son idénticos para todas sus interpolaciones, por lo que si pudiera almacenar, para cada nuevo punto de cuadrícula, los índices de los vértices del simplex adjunto y los pesos de la interpolación, minimizaría la cantidad de cálculos por mucho. Lamentablemente, esto no es fácil de hacer directamente con la funcionalidad disponible, aunque sí es posible:

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)

La función interp_weights realiza los cálculos de los primeros tres pasos que enumeré anteriormente. Luego, la función interpolate utiliza esos valores calculados para realizar el paso 4 muy rápido:

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

Entonces, primero, hace lo mismo que griddata , lo cual es bueno. En segundo lugar, la configuración de la interpolación, es decir, la computación de vtx y wts toma aproximadamente lo mismo que una llamada a griddata . Pero en tercer lugar, ahora puede interpolar para diferentes valores en la misma cuadrícula prácticamente en ningún momento.

Lo único que hace griddata que no se contempla aquí es asignar fill_value a los puntos que deben extrapolarse. Puede hacerlo revisando los puntos para los cuales al menos uno de los pesos es negativo, por ejemplo:

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

Answer 2

Muchas gracias a Jaime por su solución (incluso si no entiendo realmente cómo se realiza el cálculo baricéntrico ...)

Aquí encontrarás un ejemplo adaptado de su caso en 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))  

Es posible aplicar una transformación de imagen, como el mapeo de imágenes con una aceleración de desplazamiento

No puede usar la misma definición de función, ya que las nuevas coordenadas cambiarán en cada iteración, pero puede calcular la triangulación de una vez por todas.

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

En mi portátil, la velocidad es de entre 20 y 40 veces!

Espero que pueda ayudar a alguien.

Answer 3

Puede intentar usar Pandas , ya que proporciona estructuras de datos de alto rendimiento.

Es cierto que el método de interpolación es una envoltura de la interpolación escéptica, PERO tal vez con las estructuras mejoradas obtenga una mejor velocidad.

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

interpolate() rellena los valores de NaN en el conjunto de datos del Panel usando diferentes métodos . Espero que sea más rápido que Scipy.

Si no funciona , hay una manera de mejorar el rendimiento (en lugar de usar una versión de su código en paralelo): use Cython e implemente una pequeña rutina en C para usar dentro de su código de Python. Here tienes un ejemplo sobre esto.