加速Python/Cython cod的技巧

Question

我写了一个Python程序并对它进行了Cythonize。Cython（30%）的加速效果并不令人满意。当然可以通过改变代码结构或代码的方式来优化它。任何帮助都将不胜感激。 程序基本上是采用数字高程模型（DEM）栅格图和同一形状的过量水地图。对于多余水图中的每个像素，它搜索四个相邻像素，并确定该像素是否低于周围的相邻像素，是否具有相同的级别或高于它们。基于此，它会增加像素处的水位，或者将多余的水分配给海拔较低的邻居。这个密码一直持续到所有多余的水都散落在地上。这是Cython版本的代码。在

import numpy as np
cimport numpy as np

cdef unravel(np.ndarray[np.int_t, ndim = 1] idx,int shape0, int shape1):
    return idx//shape0, idx%shape1

cdef find_lower_neighbours(int i, int j, np.ndarray[np.double_t, ndim=2] water_level, double friction_head_loss):

    cdef double current_water_el, minlevels, deltav_total, deltav_min
    current_water_el = water_level[i,j]
    cdef np.ndarray[np.double_t, ndim = 1] levels = np.zeros(4, dtype = np.double)

    levels[:] = water_level[i - 1, j], water_level[i, j + 1], water_level[i + 1, j], water_level[i, j - 1]
    minlevels = levels.min()

    if current_water_el - minlevels < 0:
        return 0, minlevels, 0

    elif np.absolute(current_water_el - minlevels) < 0.0001:
        res = np.where(np.absolute(levels - current_water_el) < 0.0001)
        return 1, res[0], 0

    else:
        levels = current_water_el - levels
        low_values_flags = levels < 0
        levels[low_values_flags] = 0

        deltav_total = np.sum(levels)
        deltav_min = levels[levels > 0].min()
        return 2, levels / (deltav_total + deltav_min), deltav_min / (deltav_total + deltav_min)


cpdef np.ndarray[np.double_t, ndim=2] new_algorithm( np.ndarray[np.double_t, ndim=2] DEM, np.ndarray[np.double_t, ndim=2] extra_volume_map, double nodata, double pixel_area, double friction_head_loss , int von_neuman):

    cdef int terminate = 0
    cdef int iteration = 1

    cdef double sum_extra_volume_map


    index_dic_von_neuman = [[-1, 0], [0, 1], [1, 0], [0, -1]]
    cdef int DEMshape0 = DEM.shape[0]
    cdef int DEMshape1 = DEM.shape[1]

    cdef np.ndarray[np.double_t, ndim = 2] water_levels
    water_levels = np.copy(DEM)

    cdef np.ndarray[np.double_t, ndim = 2] temp_water_levels
    temp_water_levels = np.copy(water_levels)

    cdef np.ndarray[np.double_t, ndim = 2] temp_extra_volume_map
    temp_extra_volume_map = np.copy(extra_volume_map)


    cdef np.ndarray[np.int_t, ndim = 2] wetcells
    wetcells = np.zeros((DEMshape0,DEMshape1), dtype= np.int)
    wetcells[extra_volume_map > 0] = 1

    cdef np.ndarray[np.int_t, ndim = 2] temp_wetcells
    temp_wetcells = np.zeros((DEMshape0,DEMshape1), dtype= np.int )

    cdef double min_excess = friction_head_loss * pixel_area
    cdef np.ndarray[np.int_t, ndim = 1] fdx
    cdef int i, j , k, condition, have_any_dry_cells_in_neghbors
    cdef double water_level_difference, n , volume_to_each_neghbour, weight, w0

    if von_neuman == 0:
        index_dic = index_dic_von_neuman


    while terminate != 1:
        fdx = np.flatnonzero(extra_volume_map > min_excess)
        extra_volume_locations = unravel(fdx, DEMshape1,DEMshape1)
        if not extra_volume_locations[0].size:
            terminate = 1
            return water_levels

        for item in zip(*extra_volume_locations):
            i = item[0]
            j = item[1]

            if DEM[i,j] == nodata:
                print "warningggggg", i, j
                temp_extra_volume_map[i,j] = 0.
            else:
                condition, wi, w0 = find_lower_neighbours(i, j, water_levels, friction_head_loss)

                if condition == 0:
                    water_level_difference = wi - water_levels[i,j]
                    temp_water_levels[i,j] = wi
                    temp_extra_volume_map[i,j] -= water_level_difference * pixel_area

                elif condition == 1:
                    n = len(wi)
                    volume_to_each_neghbour = (extra_volume_map[i, j] - friction_head_loss * pixel_area * .02)/ (n * 1. )
                    for itemm in wi:
                        temp_extra_volume_map[i + index_dic[itemm][0], j + index_dic[itemm][1]] += volume_to_each_neghbour
                        temp_wetcells[i + index_dic[itemm][0], j + index_dic[itemm][1]] += 1
                    temp_extra_volume_map[i,j] -= extra_volume_map[i,j]
                    temp_water_levels[i,j] += friction_head_loss * .02

                elif condition == 2:
                    temp_wetcells[i,j] += 1
                    have_any_dry_cells_in_neghbors = 0 # means "no"
                    for k, weight in enumerate(wi):
                        if weight > 0:
                            temp_extra_volume_map[i + index_dic[k][0], j + index_dic[k][1]] += weight * (extra_volume_map[i,j])
                            if wetcells[i + index_dic[k][0], j + index_dic[k][1]] == 0:
                                have_any_dry_cells_in_neghbors = 1
                            temp_wetcells[i + index_dic[k][0], j + index_dic[k][1]] += 1
                    if have_any_dry_cells_in_neghbors == 1:
                        temp_water_levels[i,j] += friction_head_loss
                        temp_extra_volume_map[i, j] -= (1. - w0) * (extra_volume_map[i, j])
                    else:
                        temp_extra_volume_map[i, j] -= (1. - w0) * (extra_volume_map[i, j])


        wetcells = np.copy(temp_wetcells)
        water_levels = np.copy(temp_water_levels)
        extra_volume_map = np.copy(temp_extra_volume_map)

        iteration += 1
        if iteration*1. %500. == 0.:
            sum_extra_volume_map = np.sum(extra_volume_map)
            print "iteration", iteration, "volume left =",  sum_extra_volume_map


    print "finished at iteration =", iteration - 1

    return water_levels

以下是分析结果：

^{pr2}$

Answer 1

您首先需要profile您的程序，并精确地测量所需时间（哪些函数消耗的CPU时间最多）。在

您可以改进程序中的算法并降低某些函数的time complexity。在

你也可以手动> EME>重写你的代码中重的（资源需求，例如CPU消耗）部分，例如在（手写）C中的{a3}（也许还有C++中的某些部分或任何容易从C调用的语言，并且与Python和C的内存模型兼容）。不要指望任何工具（例如Cython）能自动有效地做到这一点。在

你的程序足够小了。考虑花几个星期（或几个月）重写它的某些部分，甚至重新设计和重写它。我猜想静态类型编译的语言（如Ocaml、Go、D、C++）可以提供一些性能改进。也许您可以考虑将代码完全重新设计为并发应用程序（multi-threaded，或OpenCL，MPI，…）

将现有的Python原型看作是开始理解您正在处理的问题的一种方法，但是重新设计并重写您的代码（可能是用另一种语言）。花点时间阅读关于这个主题的现有文献并与专家讨论。在

顺便说一句，你的用户可以接受不到一分钟的执行时间。那么你就不需要重写代码了。请记住，您的开发时间也有一些成本！有No Silver Bullet。。。在