我正在尝试定义一个包含用于模拟积分的内部循环的函数。
问题在于速度。在我的机器上,计算一次函数可能需要30秒。因为我的最终目标是最小化这个函数,所以一些额外的速度会很好。
因此,我试图让Cython正常工作,但我肯定犯了一个严重的错误(很可能是许多错误!)。按照Cython文档,我已经尝试输入我的变量。这样做之后,代码就和纯Python一样慢了。这看起来很奇怪。
下面是我的代码:
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))我已经分析了代码,似乎是函数S,积分模拟器,正在扼杀性能。无论如何,我希望通过输入变量至少能获得一些速度上的提升。由于它没有产生任何收益,我被引导相信我正在犯一些根本性的错误。
有没有人看到导致这个结果的Cython代码中的一个明显的错误?
哦,因为我对编程非常陌生,所以肯定有很多糟糕的风格和东西会减慢代码的速度。如果你有时间,也可以在这些问题上给我提个明白。
发布于 2010-11-24 18:00:38
Cython可以生成一个html文件来帮助实现这一点:
cython -a MODULE.py这显示了源代码的每一行都通过各种不同的黄色着色为白色。黄色越深,在该行上执行的Python行为就越动态。对于包含一些黄色的每一行,您需要添加更多静态类型声明。
当我这样做时,我喜欢将我遇到问题的部分源代码拆分到许多单独的行中,每个表达式或运算符一个行,以获得最细粒度的视图。
如果没有这个,很容易忽略变量、函数调用或运算符的一些静态类型声明。(例如,索引操作符xy仍然是一个完全动态的Python操作,除非您另行声明)
发布于 2010-11-24 16:57:31
Cython不提供自动的性能增益,您必须了解它的内部结构并检查生成的C代码。
特别是,如果您想要提高循环性能,就必须避免在循环中调用Python函数,在这种情况下,您碰巧做了很多事情(所有的np.调用都是Python调用、分片,可能还有其他事情)。
有关使用Cython进行性能优化的一般指南( -a开关在优化时确实很方便),请参阅this page;有关优化numpy代码时的具体细节,请参阅this one。
发布于 2010-11-24 16:32:11
你可以通过使用更多的Numpy的功能来加速你的代码。
例如:
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]将会更快更易读,因为
S = P_i.sum(axis=1)您还会重复一些计算,因此花费的时间会比必要的多两倍。例如
np.where(data[:,1]==(yr + 72))只能计算一次,并存储在可以重用的变量中。
您还可以执行大量的整形和切片:从一开始就让变量采用更简单的格式可能会有所帮助。如果可能的话,你的代码会更清晰,优化也会更明显。
https://stackoverflow.com/questions/4264229
复制相似问题