Cómo obtener la base de spline utilizada por scipy.interpolate.splev

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

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 = cv.shape[0] kv = np.clip(np.arange(c+degree+1)-degree,0,c-degree) u = np.linspace(0,c-degree,n) # Calculate result return np.array(si.splev(u, (kv,cv.T,degree))).T 

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

 # 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 ttwig de “points_scipy”: introduzca la descripción de la imagen aquí

Ahora, para hacer esto más rápido, podría calcular una base para todos los 100k puntos de la curva, almacenarlo en la memoria, y cuando necesito dibujar una curva, todo lo que tendría que hacer es multiplicar las nuevas posiciones de puntos de control con la base almacenada para obtener la nueva curva Para probar mi punto, escribí una función que utiliza 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  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, utilicemos esto para calcular una nueva base, multiplíquelo por los puntos de control y confirme que obtenemos 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ó valores base de 100k en 11 segundos, pero como estos valores solo necesitan ser computados 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 de splev me lleva a creer que internamente splev probablemente calcula una base como yo, excepto que es mucho más rápido.

Así que mi objective es descubrir cómo calcular rápidamente mi base, almacenarla en la memoria y simplemente usar np.dot () para calcular los nuevos puntos en la curva, y mi pregunta es: ¿es posible usar spicy.interpolate para obtener el ¿Valores base que (supongo) utiliza Splev para calcular su resultado? Y si es así, ¿cómo?

[APÉNDICE]

Siguiendo la muy útil información de unutbu y ev-br sobre cómo Scipy calcula la base de una spline, busqué el código de 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) + [cd]*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,cd-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+ij]) b[nrange,left+i] = b[nrange,left+i] + f * (kv[right+i] - u) b[nrange,left+i+1] = f * (u - kv[right+ij]) return b 

Probando contra la versión de unutbu de mi función 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 cual es algo de uso 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 esa razón, probablemente usaré fitpack_basis ya que es relativamente rápido. Pero me encantarían las sugerencias para mejorar su rendimiento y, con suerte, estar más cerca de la versión de unutbu de la función base original que escribí.

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

Aquí hay una versión de fitpack_basis , que funciona tan rápido como bspline_basis . Las ideas principales utilizadas para boost la velocidad con Cython es declarar el tipo de cada variable y reescribir todos los usos de NumPy Fancy Indexing como bucles utilizando la indexación de enteros.

Consulta esta página para obtener instrucciones sobre cómo comstackr 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) + [cd]*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, cd-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]+ij]) 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]+ij]) 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 de fitpack_basis ejecute tan rápido (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 

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 parte 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 su lugar, el nuevo coxDeBoor(k, d) calcula lo que estaba antes de np.array([coxDeBoor(u[i], k, d) for i in xrange(n)]) .

Dado que la aritmética de matriz NumPy tiene la misma syntax que la aritmética escalar, se necesita muy poco código para cambiar. El único cambio sintáctico estaba 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. Entonces, antes de que se devolviera un solo 0 o 1, ahora se devuelve toda una matriz 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 hacen que coxDeBoor(k, d) se llame para los mismos valores de k y d más de una vez. La syntax 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 registre en una cache una asignación de (k, d) pares de argumentos a valores devueltos. Si se coxDeBoor(k, d) , se devolverá el valor del cache lugar de volver a calcularlo.


scipy_bspline aún es más rápido, pero al menos bspline_basis plus np.dot está en el estadio de béisbol, y puede ser útil si desea 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