python - interpld - Cómo obtener la base spline utilizada por scipy.interpolate.splev

scipy interpolate extrapolate (2)

Necesito evaluar b-splines en python. Para hacerlo, escribí el siguiente código que funciona muy bien.

Darle 6 puntos de control y pedirle que evalúe 100k puntos en la curva es muy sencillo:

# Control points cv = np.array([[ 50., 25., 0.], [ 59., 12., 0.], [ 50., 10., 0.], [ 57., 2., 0.], [ 40., 4., 0.], [ 40., 14., 0.]]) n = 100000 # 100k Points degree = 3 # Curve degree points_scipy = scipy_bspline(cv,n,degree) #cProfile clocks this at 0.012 seconds

Aquí hay una trama de "points_scipy":

Ahora, para hacerlo más rápido, podría calcular una base para todos los 100k puntos en la curva, almacenar eso en la memoria, y cuando necesito dibujar una curva, todo lo que tendría que hacer es multiplicar las nuevas posiciones del punto de control con la base almacenada para obtener la nueva curva Para probar mi punto, escribí una función que usa el algoritmo de DeBoor para calcular mi base:

def basis(c, n, degree): """ bspline basis function c = number of control points. n = number of points on the curve. degree = curve degree """ # Create knot vector and a range of samples on the curve kv = np.array([0]*degree + range(c-degree+1) + [c-degree]*degree,dtype=''int'') # knot vector u = np.linspace(0,c-degree,n) # samples range # Cox - DeBoor recursive function to calculate basis def coxDeBoor(u, k, d): # Test for end conditions if (d == 0): if (kv[k] <= u and u < kv[k+1]): return 1 return 0 Den1 = kv[k+d] - kv[k] Den2 = 0 Eq1 = 0 Eq2 = 0 if Den1 > 0: Eq1 = ((u-kv[k]) / Den1) * coxDeBoor(u,k,(d-1)) try: Den2 = kv[k+d+1] - kv[k+1] if Den2 > 0: Eq2 = ((kv[k+d+1]-u) / Den2) * coxDeBoor(u,(k+1),(d-1)) except: pass return Eq1 + Eq2 # Compute basis for each point b = np.zeros((n,c)) for i in xrange(n): for k in xrange(c): b[i][k%c] += coxDeBoor(u[i],k,degree) b[n-1][-1] = 1 return b

Ahora use esto para calcular una nueva base, multiplique eso por los puntos de control y confirme que estamos obteniendo los mismos resultados que splev:

b = basis(len(cv),n,degree) #5600011 function calls (600011 primitive calls) in 10.975 seconds points_basis = np.dot(b,cv) #3 function calls in 0.002 seconds print np.allclose(points_basis,points_scipy) # Returns True

Mi función extremadamente lenta devolvió 100k valores de base en 11 segundos, pero dado que estos valores solo necesitan ser calculados una vez, el cálculo de los puntos en la curva terminó siendo 6 veces más rápido de esta manera que hacerlo a través de splev.

El hecho de que pude obtener exactamente los mismos resultados tanto de mi método como del esplendor me lleva a pensar que internamente el splend probablemente calcule una base como yo, excepto mucho más rápido.

Entonces mi objetivo es averiguar cómo calcular mi base rápidamente, almacenarla en la memoria y simplemente usar np.dot () para calcular los nuevos puntos en la curva, y mi pregunta es: ¿es posible utilizar spicy.interpolate para obtener el valores de base que (supongo) usa splev para calcular su resultado? Y si es así, ¿cómo?

[APÉNDICE]

Siguiendo el muy útil conocimiento de unutbu y ev-br sobre cómo scipy calcula la base de una spline, busqué el código fortran y escribí un equivalente a lo mejor de mis habilidades:

def fitpack_basis(c, n=100, d=3, rMinOffset=0, rMaxOffset=0): """ fitpack''s spline basis function c = number of control points. n = number of points on the curve. d = curve degree """ # Create knot vector kv = np.array([0]*d + range(c-d+1) + [c-d]*d, dtype=''int'') # Create sample range u = np.linspace(rMinOffset, rMaxOffset + c - d, n) # samples range # Create buffers b = np.zeros((n,c)) # basis bb = np.zeros((n,c)) # basis buffer left = np.clip(np.floor(u),0,c-d-1).astype(int) # left knot vector indices right = left+d+1 # right knot vector indices # Go! nrange = np.arange(n) b[nrange,left] = 1.0 for j in xrange(1, d+1): crange = np.arange(j)[:,None] bb[nrange,left+crange] = b[nrange,left+crange] b[nrange,left] = 0.0 for i in xrange(j): f = bb[nrange,left+i] / (kv[right+i] - kv[right+i-j]) b[nrange,left+i] = b[nrange,left+i] + f * (kv[right+i] - u) b[nrange,left+i+1] = f * (u - kv[right+i-j]) return b

Prueba contra la versión de unutbu de mi función de base original:

fb = fitpack_basis(c,n,d) #22 function calls in 0.044 seconds b = basis(c,n,d) #81 function calls (45 primitive calls) in 0.013 seconds ~5 times faster print np.allclose(b,fb) # Returns True

Mi función es 5 veces más lenta, pero aún relativamente rápida. Lo que me gusta de esto es que me permite usar rangos de muestra que van más allá de los límites, lo que es útil en mi aplicación. Por ejemplo:

print fitpack_basis(c,5,d,rMinOffset=-0.1,rMaxOffset=.2) [[ 1.331 -0.3468 0.0159 -0.0002 0. 0. ] [ 0.0208 0.4766 0.4391 0.0635 0. 0. ] [ 0. 0.0228 0.4398 0.4959 0.0416 0. ] [ 0. 0. 0.0407 0.3621 0.5444 0.0527] [ 0. 0. -0.0013 0.0673 -0.794 1.728 ]]

Entonces, por ese motivo, probablemente usaré fitpack_basis ya que es relativamente rápido. Pero me encantaría sugerencias para mejorar su rendimiento y espero estar más cerca de la versión de unutbu de la función base original que escribí.

Aquí hay una versión de coxDeBoor que es (en mi máquina) 840 veces más rápida que la original.

In [102]: %timeit basis(len(cv), 10**5, degree) 1 loops, best of 3: 24.5 s per loop In [104]: %timeit bspline_basis(len(cv), 10**5, degree) 10 loops, best of 3: 29.1 ms per loop

import numpy as np import scipy.interpolate as si def scipy_bspline(cv, n, degree): """ bspline basis function c = list of control points. n = number of points on the curve. degree = curve degree """ # Create a range of u values c = len(cv) kv = np.array( [0] * degree + range(c - degree + 1) + [c - degree] * degree, dtype=''int'') u = np.linspace(0, c - degree, n) # Calculate result arange = np.arange(n) points = np.zeros((n, cv.shape[1])) for i in xrange(cv.shape[1]): points[arange, i] = si.splev(u, (kv, cv[:, i], degree)) return points def memo(f): # Peter Norvig''s """Memoize the return value for each call to f(args). Then when called again with same args, we can just look it up.""" cache = {} def _f(*args): try: return cache[args] except KeyError: cache[args] = result = f(*args) return result except TypeError: # some element of args can''t be a dict key return f(*args) _f.cache = cache return _f def bspline_basis(c, n, degree): """ bspline basis function c = number of control points. n = number of points on the curve. degree = curve degree """ # Create knot vector and a range of samples on the curve kv = np.array([0] * degree + range(c - degree + 1) + [c - degree] * degree, dtype=''int'') # knot vector u = np.linspace(0, c - degree, n) # samples range # Cox - DeBoor recursive function to calculate basis @memo def coxDeBoor(k, d): # Test for end conditions if (d == 0): return ((u - kv[k] >= 0) & (u - kv[k + 1] < 0)).astype(int) denom1 = kv[k + d] - kv[k] term1 = 0 if denom1 > 0: term1 = ((u - kv[k]) / denom1) * coxDeBoor(k, d - 1) denom2 = kv[k + d + 1] - kv[k + 1] term2 = 0 if denom2 > 0: term2 = ((-(u - kv[k + d + 1]) / denom2) * coxDeBoor(k + 1, d - 1)) return term1 + term2 # Compute basis for each point b = np.column_stack([coxDeBoor(k, degree) for k in xrange(c)]) b[n - 1][-1] = 1 return b # Control points cv = np.array([[50., 25., 0.], [59., 12., 0.], [50., 10., 0.], [57., 2., 0.], [40., 4., 0.], [40., 14., 0.]]) n = 10 ** 6 degree = 3 # Curve degree points_scipy = scipy_bspline(cv, n, degree) b = bspline_basis(len(cv), n, degree) points_basis = np.dot(b, cv) print np.allclose(points_basis, points_scipy)

La mayoría de la aceleración se logra al hacer que coxDeBoor calcule un vector de resultados en lugar de un solo valor a la vez. Observe que u se elimina de los argumentos pasados a coxDeBoor . En cambio, el nuevo coxDeBoor(k, d) calcula lo que estaba antes de np.array([coxDeBoor(u[i], k, d) for i in xrange(n)]) .

Como la aritmética de matrices NumPy tiene la misma sintaxis que la aritmética escalar, se necesita muy poco código para cambiar. El único cambio sintáctico fue en la condición final:

if (d == 0): return ((u - kv[k] >= 0) & (u - kv[k + 1] < 0)).astype(int)

(u - kv[k] >= 0) y (u - kv[k + 1] < 0) son matrices booleanas. astype cambia los valores de la matriz a 0 y 1. Por lo tanto, antes de que se devolviera un único 0 o 1, ahora se devuelve una matriz completa de 0 y 1, uno para cada valor en u .

La memorización también puede mejorar el rendimiento, ya que las relaciones de recurrencia provocan que se coxDeBoor(k, d) por los mismos valores de coxDeBoor(k, d) más de una vez. La sintaxis del decorador

@memo def coxDeBoor(k, d): ...

es equivalente a

def coxDeBoor(k, d): ... coxDeBoor = memo(coxDeBoor)

y el decorador de memo hace que coxDeBoor grabe en un cache un mapeo de (k, d) pares de argumentos a valores devueltos. Si se llama de nuevo a coxDeBoor(k, d) , se devolverá el valor de la cache lugar de volver a calcular el valor.

scipy_bspline es aún más rápido, pero al menos bspline_basis plus np.dot está en el estadio, y puede ser útil si quieres reutilizar b con muchos puntos de control, cv .

In [109]: %timeit scipy_bspline(cv, 10**5, degree) 100 loops, best of 3: 17.2 ms per loop In [104]: %timeit b = bspline_basis(len(cv), 10**5, degree) 10 loops, best of 3: 29.1 ms per loop In [111]: %timeit points_basis = np.dot(b, cv) 100 loops, best of 3: 2.46 ms per loop

fitpack_basis usa un bucle doble que modifica iterativamente elementos en bb y b . No veo una forma de usar NumPy para vectorizar estos bucles, ya que el valor de bb b en cada etapa de la iteración depende de los valores de las iteraciones anteriores. En situaciones como esta, a veces se puede usar Cython para mejorar el rendimiento de los bucles.

Aquí hay una versión Cython fitpack_basis de fitpack_basis , que se ejecuta tan rápido como bspline_basis . Las ideas principales utilizadas para aumentar la velocidad con Cython es declarar el tipo de cada variable y volver a escribir todos los usos de la indexación de fantasía de NumPy como bucles utilizando la indexación de entero simple.

Consulte esta página para obtener instrucciones sobre cómo crear el código y ejecutarlo desde Python.

import numpy as np cimport numpy as np cimport cython ctypedef np.float64_t DTYPE_f ctypedef np.int64_t DTYPE_i @cython.boundscheck(False) @cython.wraparound(False) @cython.nonecheck(False) def cython_fitpack_basis(int c, int n=100, int d=3, double rMinOffset=0, double rMaxOffset=0): """ fitpack''s spline basis function c = number of control points. n = number of points on the curve. d = curve degree """ cdef Py_ssize_t i, j, k, l cdef double f # Create knot vector cdef np.ndarray[DTYPE_i, ndim=1] kv = np.array( [0]*d + range(c-d+1) + [c-d]*d, dtype=np.int64) # Create sample range cdef np.ndarray[DTYPE_f, ndim=1] u = np.linspace( rMinOffset, rMaxOffset + c - d, n) # basis cdef np.ndarray[DTYPE_f, ndim=2] b = np.zeros((n,c)) # basis buffer cdef np.ndarray[DTYPE_f, ndim=2] bb = np.zeros((n,c)) # left knot vector indices cdef np.ndarray[DTYPE_i, ndim=1] left = np.clip(np.floor(u), 0, c-d-1).astype(np.int64) # right knot vector indices cdef np.ndarray[DTYPE_i, ndim=1] right = left+d+1 for k in range(n): b[k, left[k]] = 1.0 for j in range(1, d+1): for l in range(j): for k in range(n): bb[k, left[k] + l] = b[k, left[k] + l] b[k, left[k]] = 0.0 for i in range(j): for k in range(n): f = bb[k, left[k]+i] / (kv[right[k]+i] - kv[right[k]+i-j]) b[k, left[k]+i] = b[k, left[k]+i] + f * (kv[right[k]+i] - u[k]) b[k, left[k]+i+1] = f * (u[k] - kv[right[k]+i-j]) return b

Usando este código de tiempo para comparar su rendimiento,

import timeit import numpy as np import cython_bspline as CB import numpy_bspline as NB c = 6 n = 10**5 d = 3 fb = NB.fitpack_basis(c, n, d) bb = NB.bspline_basis(c, n, d) cfb = CB.cython_fitpack_basis(c,n,d) assert np.allclose(bb, fb) assert np.allclose(cfb, fb) # print(NB.fitpack_basis(c,5,d,rMinOffset=-0.1,rMaxOffset=.2)) timing = dict() timing[''NB.fitpack_basis''] = timeit.timeit( stmt=''NB.fitpack_basis(c, n, d)'', setup=''from __main__ import NB, c, n, d'', number=10) timing[''NB.bspline_basis''] = timeit.timeit( stmt=''NB.bspline_basis(c, n, d)'', setup=''from __main__ import NB, c, n, d'', number=10) timing[''CB.cython_fitpack_basis''] = timeit.timeit( stmt=''CB.cython_fitpack_basis(c, n, d)'', setup=''from __main__ import CB, c, n, d'', number=10) for func_name, t in timing.items(): print "{:>25}: {:.4f}".format(func_name, t)

parece que Cython puede hacer que el código fitpack_basis ejecute tan rápido como (y quizás un poco más rápido) que bspline_basis :

NB.bspline_basis: 0.3322 CB.cython_fitpack_basis: 0.2939 NB.fitpack_basis: 0.9182