Sin ganancias de velocidad de Cython

Question

Estoy tratando de definir una función que contenga un bucle interno para simular una integral.

El problema es la velocidad. Evaluar la función una vez puede tomar hasta 30 segundos en mi máquina. Como mi objetivo final es minimizar esta función, una velocidad extra sería agradable.

Como tal, he tratado de hacer que Cython funcione, pero debo cometer un grave error (¡probablemente muchos de ellos!). Siguiendo la documentación de Cython, intenté escribir mis variables. Después de hacerlo, el código es tan lento como Python puro. Esto parece extraño.

Aquí está mi código:

import numpy as np 
cimport cython 
cimport numpy as np 
import minuit 

data = np.genfromtxt('q6data.csv', usecols = np.arange(1, 24, 1), delimiter = ',') 

cdef int ns = 1000     # Number of simulation draws 
cdef int K  = 5     # Number of observed characteristics, including   constant 
cdef int J  = len(data[:,1])  # Number of products, including outside 
cdef double tol = 0.0001   # Inner GMM loop tolerance 
nu = np.random.normal(0, 1, (6, ns)) # ns random deviates 

@cython.boundscheck(False) 
@cython.wraparound(False) 


def S(np.ndarray[double, ndim=1] delta, double s1, double s2, double s3, double s4, double s5, double a): 
    """Computes the simulated integrals, one for each good. 
    Parameters: delta is an array of length J containing mean product specific utility levels 
    Returns: Numpy array with length J.""" 
    cdef np.ndarray[double, ndim=2] mu_ij = np.dot((data[:,2:7]*np.array([s1, s2, s3, s4, s5])), nu[1:K+1,:]) 
    cdef np.ndarray[double, ndim=2] mu_y = a * np.log(np.exp(data[:,21].reshape(J,1) + data[:,22].reshape(J,1)*nu[0,:].reshape(1, ns)) - data[:,7].reshape(J,1)) 
    cdef np.ndarray[double, ndim=2] V = delta.reshape(J,1) + mu_ij + mu_y 
    cdef np.ndarray[double, ndim=2] exp_vi = np.exp(V) 
    cdef np.ndarray[double, ndim=2] P_i = (1.0/np.sum(exp_vi[np.where(data[:,1] == 71)], 0)) * exp_vi[np.where(data[:,1] == 71)] 
    cdef int yrs = 19 
    cdef int yr 
    for yr in xrange(yrs): 
     P_yr = (1.0/np.sum(exp_vi[np.where(data[:,1]== (yr + 72))], 0)) * exp_vi[np.where(data[:,1] == (yr + 72))] 
     P_i = np.concatenate((P_i, P_yr)) 
    cdef np.ndarray[double, ndim=1] S = np.zeros(dtype = "d", shape = J) 
    cdef int j 
    for j in xrange(ns): 
     S += P_i[:,j] 
    return (1.0/ns) * S 

def d_infty(np.ndarray[double, ndim=1] x, np.ndarray[double, ndim=1] y): 
    """Sup norm.""" 
    return np.max(np.abs(x - y)) 

def T(np.ndarray[double, ndim=1] delta_exp, double s1, double s2, double s3, double s4, double s5, double a): 
    """The contraction operator. This function takes the parameters and the exponential 
    of the starting value of delta and returns the fixed point.""" 
    cdef int iter = 0 
    cdef int maxiter = 200 
    cdef int i 
    for i in xrange(maxiter): 
     delta1_exp = delta_exp * (data[:, 8]/S(np.log(delta_exp), s1, s2, s3, s4, s5, a))      
     print i 
     if d_infty(delta_exp, delta1_exp) < tol:          
      break 
     delta_exp = delta1_exp 
    return np.log(delta1_exp) 


def Q(double s1, double s2, double s3, double s4, double s5, double a): 
    """GMM objective function.""" 
    cdef np.ndarray[double, ndim=1] delta0_exp = np.exp(data[:,10])              
    cdef np.ndarray[double, ndim=1] delta1 = T(delta0_exp, s1, s2, s3, s4, s5, a) 
    delta1[np.where(data[:,10]==0)] = np.zeros(len(np.where(data[:,10]==0)))    
    cdef np.ndarray[double, ndim=1] xi = delta1 - (np.dot(data[:,2:7], np.linalg.lstsq(data[:,2:7], delta1)[0])) 
    cdef np.ndarray[double, ndim=2] g_J = xi.reshape(J,1) * data[:,11:21] 
    cdef np.ndarray[double, ndim=1] G_J = (1.0/J) * np.sum(g_J, 0) 
    return np.sqrt(np.dot(G_J, G_J)) 

he perfilado el código, y que parece ser la función S, el simulador integral, que está matando el rendimiento. De todos modos, habría esperado al menos algunas ganancias de velocidad al escribir mis variables. Como no produjo ganancias, me hacen creer que estoy cometiendo algunos errores fundamentales.

¿Alguien ve un error evidente en el código de Cython que podría conducir a este resultado?

Ah, y como soy muy nuevo en la programación, seguramente hay un montón de malos estilos y cosas que ralentizan el código. Si tienes tiempo, no dudes en aclararme estos puntos también.

Answer 1

El "error fundamental" es que espera un buen rendimiento en bucles largos de Python. Es un lenguaje interpretado, y el cambio entre implementaciones y codificación no hace nada con esto. Hay algunas bibliotecas numéricas de Python para computación rápida, principalmente en C. Por ejemplo, si ya usas numpy para matrices, ¿por qué no vas más allá y usas scipy para tus matemáticas avanzadas? Aumentará la velocidad de lectura y.

Answer 2

Definitivamente puede acelerar su código utilizando más capacidades de Numpy.

Por ejemplo:

cdef np.ndarray[double, ndim=1] S = np.zeros(dtype = "d", shape = J) 
cdef int j 
for j in xrange(ns): 
    S += P_i[:,j] 

sería mucho más rápido y más legible

S = P_i.sum(axis=1) 

También repetir algunos cálculos, que por lo tanto tienen dos veces más tiempo de lo necesario. Por ejemplo,

np.where(data[:,1]==(yr + 72)) 

podría calcularse solo una vez y almacenarse en una variable que pueda volver a utilizar.

También realiza una gran cantidad de remodelación y división: podría ayudar a tener sus variables en un formato más simple desde el principio. Si es posible, su código sería mucho más claro, y las optimizaciones podrían ser mucho más obvias.

Answer 3

Cython no ofrece aumentos de rendimiento automáticos, debe conocer sus partes internas y verificar el código C generado.

En particular, si desea mejorar los rendimientos de los bucles, debe evitar invocar funciones de Python en ellos, lo cual sucede mucho en este caso (todas las llamadas np. son llamadas Python, slicing y probablemente otras) .

Consulte this page para obtener pautas generales sobre la optimización del rendimiento con Cython (el modificador -a es realmente útil al optimizar) y this one para conocer las características específicas al optimizar el código numpy.

Answer 4

¿Le ayudará un perfilador a determinar qué pieza es lenta? Me gusta correr programas utilizando el generador de perfiles biblioteca estándar:

python -O -m cProfile -o profile.out MYAPP.py 

y luego ver el resultado de que en el 'RunSnakeRun' GUI:

runsnake profile.out 

Un RunSnakeRun se puede instalar desde aquí: http://www.vrplumber.com/programming/runsnakerun/

Answer 5

Cython puede producir un archivo hTML para ayudar con esto:

cython -a MODULE.py 

Esto muestra cada línea de código fuente de color blanco a través de varios tonos de amarillo. Cuanto más oscuro es el color amarillo, más comportamiento dinámico de Python se sigue realizando en esa línea. Para cada línea que contenga algo de amarillo, debe agregar más declaraciones de tipado estáticas.

Cuando hago esto, me gusta dividir partes de mi código fuente con las que estoy teniendo problemas en muchas líneas separadas, una para cada expresión u operador, para obtener la vista más granular.

Sin esto, es fácil pasar por alto algunas declaraciones de tipo estático de variables, llamadas a funciones u operadores. (Por ejemplo, el operador de indexación x [y] sigue siendo una operación de Python completamente dinámica a menos que usted declare lo contrario)

Answer 6

Tomando el consejo dado aquí, he dedicado más tiempo a perfilar el código anterior. Para poder limpiar un poco las cosas, definí

He perfilado el código un poco más y tengo una mejor idea de qué piezas son las más lentas. Definí adicionalmente

X = data[:, 2:7] 
m_y = data[:, 21].reshape(J,1) 
sigma_y = 1.0 
price = data[:, 7].reshape(J, 1) 
shares_data = data[:,8] 

Las siguientes líneas ocupan la mayor parte del tiempo total.

mu_ij = np.dot((X*np.array([s1, s2, s3, s4, s5])), nu[1:K+1,:]) 
mu_y = a * np.log(np.exp(m_y + sigma_y*nu[0,:].reshape(1,ns)) - price) 
V = delta.reshape(J,1) + mu_ij + mu_y 
exp_vi = np.exp(V) 
P_i = (1.0/np.sum(exp_vi[np.where(data[:,1]==71)], 0)) * exp_vi[np.where(data[:,1]==71)] 
for yr in xarange(19): 
    P_yr = (1.0/np.sum(exp_vi[np.where(data[:,1]==yr)], 0)) * exp_vi[np.where(data[:,1]==yr)] 
P_i = np.concatenate((P_i, P_yr)) 

Me da la impresión de que esta es una forma demasiado engorrosa para lograr mi objetivo. Esperaba que alguien pudiera darme algunos consejos sobre cómo acelerar estas líneas. Tal vez hay capacidades Numpy que me estoy perdiendo? Si este problema no está lo suficientemente bien especificado para que usted sea útil, me complacería brindarle más detalles sobre el contexto de mi problema. ¡Gracias!

Answer 7

dividir los datos única, después de su comentario el 28 de noviembre:

import sys 
import time 
import numpy as np 

def splitdata(data, n, start=1971): 
    """ split data into n pieces, where col 1 == start .. start + n """ 
     # not fancy, fast enough for small n 
    split = n * [None] 
    for j in range(n): 
     split[j] = data[ data[:,1] == start + j ] 
    return split # [ arrays: col1 0, col1 1 ... ] 

#........................................................................... 
N = 2237 
ncol = 21 
start = 1971 
n = 20 
seed = 1 
exec "\n".join(sys.argv[1:]) # run this.py N= ... 
np.set_printoptions(2, threshold=100, suppress=True) # .2f 
np.random.seed(seed) 

print "N=%d ncol=%d n=%d" % (N, ncol, n) 
data = np.random.uniform(start, start + n, (N,ncol)) 
data[:,1] = data[:,1].round() 
t0 = time.time() 

split = splitdata(data, n, start) # n pieces 

print "time: %.0f us splitdata" % ((time.time() - t0) * 1e6) 
for y, yeardata in enumerate(split): 
    print "%d %d %g" % (start + y, len(yeardata), yeardata[:,0].sum()) 

->

time: 27632 us splitdata # old mac ppc 
1971 69 136638 
1972 138 273292 
...