提高循环性能的速度

Question

我正在尝试使用 pyMC 为我的马尔可夫链蒙特卡罗代码构建一个示例。因此使用模型的采样参数，每次输出都是通过从 nfw class 调用getLensing instance 构建的，并与观察到的数据进行比较。我的问题是我的代码在计算模型参数时非常慢。例如，我有 24000 个数据点，然后对于每个数据点，我都有一个概率分布 - 例如obj_pdf- 我在内部循环中将其边缘化（集成）。所以每次至少需要一个小时来计算模型的所有输出。

import numpy as np
z=np.arange(0,1.5,0.001)
z_h=0.15
for j in range(pos.shape[0]):
    value1=0;value2=0
    pdf=obj_pdf[j,:]/sum(obj_pdf[j,:])
    for i in range(len(z)):
        if (z[i]>z_h) :
        g1,g2=nfw.getLensing( pos[j,:], z[i])
            value1+=g1*pdf[i]
            value2+=g2*pdf[i]
    if (j<1):
       value=np.array([value1,value2])
    else:
       value=np.vstack((value, np.array([value1,value2]))) 

因此，如果我想重新采样输入参数，例如 100000 次，则需要几个月的时间来进行 MCMC 计算。有什么聪明的方法可以加速我的代码和循环吗？ 我是否需要使用numpy.vectorize 之类的东西，或者它仍然不会提高我的代码速度？ cython 怎么样，它会提高代码的性能吗？如果有帮助，它是如何工作的？

我跑了python -m cProfile mycode.py 看看是什么原因导致我的代码变慢，结果如下：

    12071    0.004    0.000    0.004    0.000 {min}
        2    0.000    0.000    0.000    0.000 {next}
        1    0.000    0.000    0.000    0.000 {numexpr.interpreter._set_num_threads}
        8    0.002    0.000    0.002    0.000 {numpy.core.multiarray.arange}
132424695  312.210    0.000  312.210    0.000 {numpy.core.multiarray.array}
    73498    3.933    0.000    3.933    0.000 {numpy.core.multiarray.concatenate}
 99151506  201.497    0.000  201.497    0.000 {numpy.core.multiarray.copyto}
 99151500  164.303    0.000  164.303    0.000 {numpy.core.multiarray.empty_like}
       28    0.000    0.000    0.000    0.000 {numpy.core.multiarray.empty}
        2    0.000    0.000    0.000    0.000 {numpy.core.multiarray.set_string_function}
        1    0.000    0.000    0.000    0.000 {numpy.core.multiarray.set_typeDict}
        1    0.000    0.000    0.000    0.000 {numpy.core.multiarray.where}
       14    0.000    0.000    0.000    0.000 {numpy.core.multiarray.zeros}
       14    0.000    0.000    0.000    0.000 {numpy.core.umath.geterrobj}
        7    0.000    0.000    0.000    0.000 {numpy.core.umath.seterrobj}
      270    0.000    0.000    0.000    0.000 {numpy.lib._compiled_base.add_docstring}
        6    0.011    0.002    0.011    0.002 {open}
        1    0.000    0.000    0.000    0.000 {operator.div}
        2    0.000    0.000    0.000    0.000 {operator.mul}
     1918    0.000    0.000    0.000    0.000 {ord}
        2    0.000    0.000    0.000    0.000 {posix.WEXITSTATUS}
        2    0.000    0.000    0.000    0.000 {posix.WIFEXITED}
        1    0.000    0.000    0.000    0.000 {posix.WIFSIGNALED}
        9    0.002    0.000    0.002    0.000 {posix.access}
        3    0.000    0.000    0.000    0.000 {posix.close}
        5    0.002    0.000    0.002    0.000 {posix.fdopen}
        1    0.002    0.002    0.002    0.002 {posix.fork}
        4    0.000    0.000    0.000    0.000 {posix.getcwd}
        6    0.000    0.000    0.000    0.000 {posix.getpid}
        1    0.000    0.000    0.000    0.000 {posix.getuid}
        1    0.000    0.000    0.000    0.000 {posix.listdir}
        6    0.000    0.000    0.000    0.000 {posix.lstat}
        4    0.043    0.011    0.043    0.011 {posix.open}
        2    0.000    0.000    0.000    0.000 {posix.pipe}
        2    0.004    0.002    0.004    0.002 {posix.popen}
        1    0.007    0.007    0.007    0.007 {posix.read}
      205    0.059    0.000    0.059    0.000 {posix.stat}
        3    0.000    0.000    0.000    0.000 {posix.sysconf}
        2    0.000    0.000    0.000    0.000 {posix.uname}
        4    0.004    0.001    0.004    0.001 {posix.unlink}
        3    0.000    0.000    0.000    0.000 {posix.urandom}
        1    0.000    0.000    0.000    0.000 {posix.waitpid}
        1    0.000    0.000    0.000    0.000 {pow}
     1522    0.004    0.000    0.004    0.000 {range}
       73    0.000    0.000    0.000    0.000 {repr}
 99151501 2102.879    0.000 6380.906    0.000 {scipy.integrate._quadpack._qagse}
     1776    0.002    0.000    0.002    0.000 {setattr}
       32    0.000    0.000    0.000    0.000 {sorted}
    24500   18.861    0.001   18.861    0.001 {sum}
      184    0.000    0.000    0.000    0.000 {sys._getframe}
        1    0.000    0.000    0.000    0.000 {sys.getfilesystemencoding}
        2    0.000    0.000    0.000    0.000 {sys.settrace}
        1    0.000    0.000    0.000    0.000 {tables.utilsextension._broken_hdf5_long_double}
        1    0.000    0.000    0.000    0.000 {tables.utilsextension.blosc_compressor_list}
        2    0.000    0.000    0.000    0.000 {tables.utilsextension.get_hdf5_version}
        1    0.000    0.000    0.000    0.000 {tables.utilsextension.get_pytables_version}
        2    0.000    0.000    0.000    0.000 {tables.utilsextension.which_lib_version}
       27    0.000    0.000    0.000    0.000 {thread.allocate_lock}
        6    0.000    0.000    0.000    0.000 {thread.get_ident}
        4    0.000    0.000    0.000    0.000 {thread.start_new_thread}
        1    0.000    0.000    0.000    0.000 {time.localtime}
        2    0.000    0.000    0.000    0.000 {time.time}
      105    0.000    0.000    0.000    0.000 {unichr}
      229    0.000    0.000    0.000    0.000 {vars}
    49300    2.127    0.000    2.127    0.000 {zip}

Answer 1

这是一些代码。如果时间从 60 分钟变为 59 分钟，我会感到惊讶。

import numpy as np
z_h=0.15
z=np.arange(z_h, 1.5,0.001) #start the range from what you need (not exactly
z=z[1:] # needed because you said if (z[i]>z_h), range gives (z[i]>=z_h)

value=np.array([])

for j in range(pos.shape[0]):
    value1=0;value2=0
    pdf=obj_pdf[j,:]/sum(obj_pdf[j,:])
    posj=pos[j,:] #precalculate 
    for i,zi in enumerate(z): #use enumerate if you need value and index
        g1,g2=nfw.getLensing( posj, zi)
        value1+=g1*pdf[i]
        value2+=g2*pdf[i]
    value=np.append(value, np.array([value1,value2])) # use a proper append function

与其他人一样，我认为 getLensing 正在占用您的 CPU 周期。

根据this 的第一个答案，np.vectorize 不会加速你的函数。