计算分子在顶点处的推力或拉力

Question

代码的基本思想是计算顶点处的推力或拉力，给定“推致分子”和“拉致分子”在多边形顶点的数目。

因此，代码主要执行算术(包括它调用的函数calculate_edge_force)，并且有几个if-statements。函数被多次调用，因为它是计算导数所必需的，后来被scipy.odeint用于集成。

下面是我用普通Python编写的代码(请注意，import numpy as np放在文件的较高位置，函数位于类中，因此引用为self.<something>)：

def calculate_point_forces(self, PushMolecules, PushMolecules_iPlus1, PushMolecules_iMinus1, PullMolecules, PullMolecules_iPlus1, PullMolecules_iMinus1, point_index, num_nodes, polygon, polygon_compression, point_interpolygon_contact, push_vector, pull_vector, a1, a2):

    max_push_force = self.max_push_force
    max_pull_force = self.max_pull_force
    tangent_factor = self.tangent_factor
    rac_force_exponent = self.rac_force_exponent
    rho_force_exponent = self.rho_force_exponent

    avg_PushMolecules_plus1 = (PushMolecules + PushMolecules_iPlus1)*0.5
    avg_PullMolecules_plus1 = (PullMolecules + PullMolecules_iPlus1)*0.5

    avg_PushMolecules_minus1 = (PushMolecules + PushMolecules_iMinus1)*0.5
    avg_PullMolecules_minus1 = (PullMolecules + PullMolecules_iMinus1)*0.5

    if self.enable_compression_stiffening == True:
        if polygon_compression < 1:
            F_internal = ((1-polygon_compression)+self.cytosolic_baseline)*self.cytosolic_pressure
        else:
            F_internal = 0
    else:
        F_internal = 0

    if PushMolecules > PullMolecules:
        Fpush = int(not point_interpolygon_contact)*max_push_force*(1 - (1/(1 + (PushMolecules/a2)**rac_force_exponent)))*push_vector
        Fpull = 0
        Fe_pull = 0
    elif PushMolecules < PullMolecules:
        pull_magnitude = max_pull_force*(1 - (1/(1 + (PullMolecules/a1)**rho_force_exponent)))
        Fpull = (1-tangent_factor)*pull_magnitude*pull_vector
        Fe_pull = 0.5*tangent_factor*pull_magnitude
        Fpush = 0
    else:
        Fpush = 0
        Fpull = 0
        Fe_pull = 0

    Fe_pull_minus1 = 0
    Fe_pull_plus1 = 0
    if PushMolecules_iMinus1 < PullMolecules_iMinus1:
        pull_magnitude_iMinus1 = max_pull_force*(1 - (1/(1 + (PullMolecules_iMinus1/a1)**rho_force_exponent)))
        Fe_pull_minus1 = 0.5*tangent_factor*pull_magnitude_iMinus1
    if PushMolecules_iPlus1 < PullMolecules_iPlus1:
        pull_magnitude_iPlus1 = max_pull_force*(1 - (1/(1 + (PullMolecules_iPlus1/a1)**rho_force_exponent)))
        Fe_pull_plus1 = 0.5*tangent_factor*pull_magnitude_iPlus1

    Fe_pull  = np.max([Fe_pull_minus1, Fe_pull_plus1, Fe_pull])

    Fe_plus, Fe_plus_ = self.calculate_edge_force(
        avg_PushMolecules_plus1, avg_PullMolecules_plus1, point_index,
        (point_index+1) % num_nodes, polygon)
    Fe_minus, Fe_minus_ = self.calculate_edge_force(
        avg_PushMolecules_minus1, avg_PullMolecules_minus1, point_index,
        (point_index-1) % num_nodes, polygon)

    F = Fpush + Fpull + Fe_plus + Fe_minus + Fe_pull*Fe_plus_ + Fe_pull*Fe_minus_ + push_vector*F_internal

    return F

当我分析这一点时，我得到函数本身需要大约3秒的时间。其中，1秒用于调用calculate_edge_force，1秒用于np.max。

如果我“天真地”把float和int放到我能做的任何地方(我不知道如何键入NumPy数组，例如push_vector或pull_vector，所以我只为像PushMolecules、max_push_force或point_index这样的东西(我知道这将是float或int)这样做)，我实际上会得到一个非常温和的速度减速(大约0.2秒)！

显然，这不是我的代码“细胞化”的方法。我该怎么做呢？

Answer 1

首先，您可能应该键入numpy数组，但是要注意，无论您使用np.ndarray表示法还是内存视图表示法，只会以c的速度进行有限的操作，即索引、“shape”和其他一些操作。使用内存视图，您只能使用这些快速操作，使用ndarray类型，您也可以访问所有python操作，但它们仍然以python的速度进行，而不是c。

在这种情况下，您的代码实际上没有使用任何数组函数，这些数组函数将受益于正确的数组类型，因此只需键入数组就不会获得任何好处。调用python函数的开销很大，因此不使用numpy数组函数，如果numpy数组很短，如果功能使用数组索引(fast)直接作为cython代码再现，如果函数多次被调用，您可能会从中受益。分析人员应该告诉你这些事。

在上面的代码中，我可以看到一些低挂的水果，取这一行，被提到是一个主要的罪犯：

Fe_pull  = np.max([Fe_pull_minus1, Fe_pull_plus1, Fe_pull])

这里发生的情况是，np.max是对python模块的python函数的调用，该列表将被构造为python列表，并且由于python列表不能容纳c floats，所以浮点数将被转换为python。然后，这些代码都将以python的速度运行。然后，np.max将解压缩所有这些内容，np.max的结果将是一个python，它被转换为c浮点数。这会很慢的！事实上，它甚至可能比python慢，因为它涉及创建额外的临时python对象。

您要做的是重写它以消除python函数调用(从而消除所有中间python对象)。碰巧的是，如果您在类型化变量上使用内置的python 'max‘，cython就足够聪明地完全消除python调用：

Fe_pull = max(Fe_pull_minus1, Fe_pull_plus1, Fe_pull)

这将自动成为计算3个变量最大值的最佳内联c代码(Cython实际上像if/else语句一样展开它)

下一个低挂水果，另一个被指控的违法者是：

Fe_plus, Fe_plus_ = self.calculate_edge_force( ...

在某些情况下，Cython优化了元组解压缩，但在函数调用的情况下并非如此。所以将要发生的事情是cython将创建一个中间元组( python对象)，将内容打包到tuple中(同样，转换为python对象，因为元组不能容纳c类型)，然后解压元组并将python对象转换为它们的c类型。同样，这应该转化为易于转换为c的内容。不幸的是，在这种情况下，没有更快、更干净的语法，最好是通过创建“对”结构并手动打包和解压缩，比如：

#At the top level
cdef struct float_pair:
    float a
    float b

# In the class method definition
cdef floar_pair calculate_edge_force(...
    ...
    cdef float_pair result
    result.a = ...
    result.b = ...
    return result


# In your function   
    cdef float_pair result = self.calculate_edge_force( ...
    Fe_plus = result.a
    Fe_plus_ = result.b

在这里，使用结构绝不是唯一的选择，但它将比使用一个小数组更快。您还可以使用C++的“对”模板，这个模板与结构一样快，但它需要引入C++。在这个问题上，你也可以通过指针传递浮标。这些代码在速度和优雅方面都没有明显的优势，但它们都会编译成干净的c代码。

我想，如果您根除了python调用中的那些大恶棍，您的速度就会有显著的提高。之后，在.pyx文件上使用cython -a就成了一个问题。生成的html文件以黄色/橙色显示“未优化”行，您可以单击一行查看生成的c代码。从这里，您可以了解如何消除python调用--如果值得的话。

Answer 2

当你有这样的代码

cdef Type x, y, z
z = some_python_function(x, y)

Cython不会创建一个C类型，然后将其转换为Python类型。这样做后，它会将答案转换回来。这两个转换步骤都可能是缓慢的，并且可能导致您提到的0.2秒的减速。

看一下您的代码，很难判断是什么类型的东西。有一点要注意的是如果

avg_PushMolecules_plus1 = (PushMolecules + PushMolecules_iPlus1)*0.5

正在对Numpy数组进行操作，键入这些将不会给您带来任何速度改进。相反，你应该看看

• 打开数组表达式并键入它们，然后使用Cython进行编译。
• 打开数组表达式并使用Numba的autojit
• 使用类似于numexpr的方法来优化计算

甚至只是重写了这个

avg_PushMolecules_plus1 = PushMolecules + PushMolecules_iPlus1
acg_PushMolecules_plus1 *= 0.5

可能有一个小的积极作用，通过去除中间体。

现在，以下是重要的部分：

在函数的内部，您不会做很重的循环。就目前情况而言，这一职能没有任何需要加快的地方。

值得做的一件事是使用line_profiler获取每行的时间。这将向您展示哪些表达式最需要优化。

你提到numpy.max是一个潜在的候选人。优化此调用的一种方法是使用纯Cython编写它：

cdef double Fe_pull_minus1, Fe_pull_plus1, Fe_pull
...
cdef double maximum = Fe_pull_minus1
if Fe_pull_plus1 > maximum: maximum = Fe_pull_plus1
if Fe_pull > maximum: maximum = Fe_pull

虽然这听起来很可疑，但这样的电话可能会对时间产生重大影响--当然不会超过几微秒！

慢的部分在哪里并不明显，因为

• 我不知道什么，如果有的话，是一个数组，
• 我不知道这是多久一次，所以我不知道什么样的优化是合适的，
• 我没有时间。

您可以使用line_profiler解决所有这些困惑点。

安装它并放在文件的顶部：

import line_profiler
profiler = line_profiler.LineProfiler()

然后在你运行任何东西之前，写：

profiler.enable()
profiler.add_function(calculate_point_forces)

最后，在脚本结束时，编写：

profiler.print_stats()

然后，输出应该包括函数的基于行的时间。