《Python高性能编程》——2.3 计算完整的Julia集合

本节书摘来自异步社区《Python高性能编程》一书中的第2章，第2.3节，作者[美] 戈雷利克 （Micha Gorelick），胡世杰，徐旭彬 译，更多章节内容可以访问云栖社区“异步社区”公众号查看。

2.3　计算完整的Julia集合

我们在本节分解Julia集合的生成代码。我们将在本章以各种方法分析它。如例2-1所示，在模块的一开始，我们导入time模块作为我们的第一种分析手段并定义一些坐标常量。

例2-1　定义空间坐标的全局常量

"""Julia set generator without optional PIL-based image drawing"""
import time

# area of complex space to investigate
x1, x2, y1, y2 = -1.8, 1.8, -1.8, 1.8
c_real, c_imag = -0.62772, -.42193

为了生成图像，我们创建两列输入数据。第一列zs（复数z的坐标），第二列cs（复数的初始条件）。这两列数据都不会发生变化，且cs列表可以被优化成一个c常量。建立两列输入的理由是之后当我们在本章分析RAM使用情况时可以有一些看上去合理的数据来进行分析。

为了创建zs列表和cs列表，我们需要知道每个z的坐标。例2-2中，我们用xcoord和ycoord以及指定的x_step和y_step创建了这些坐标。这样一个冗长的设定方式有助于将代码移植到其他工具（如numpy）或Python环境上，因为它将一切都明确定义下来，方便除错。

例2-2　建立坐标列表作为计算函数的输入

def calc_pure_python(desired_width, max_iterations):
    """Create a list of complex coordinates (zs) and complex
    parameters (cs), build Julia set, and display"""
    x_step = (float(x2 - x1) / float(desired_width))
    y_step = (float(y1 - y2) / float(desired_width))
    x = []
    y = []
    ycoord = y2
    while ycoord > y1:
        y.append(ycoord)
        ycoord += y_step
    xcoord = x1
    while xcoord < x2:
        x.append(xcoord)
        xcoord += x_step
    # Build a list of coordinates and the initial condition for each cell.
    # Note that our initial condition is a constant and could easily be removed;
    # we use it to simulate a real-world scenario with several inputs to
    # our function.
    zs = []
    cs = []
    for ycoord in y:
        for xcoord in x:
            zs.append(complex(xcoord, ycoord))
            cs.append(complex(c_real, c_imag))

    print "Length of x:", len(x)
    print "Total elements:", len(zs)
    start_time = time.time()
    output = calculate_z_serial_purepython(max_iterations, zs, cs)
    end_time = time.time()
    secs = end_time - start_time
    print calculate_z_serial_purepython.func_name + " took", secs, "seconds"

    # This sum is expected for a 1000^2 grid with 300 iterations.
    # It catches minor errors we might introduce when we're
    # working on a fixed set of inputs.
    assert sum(output) == 33219980

创建了zs列表和cs列表后，我们输出一些有关列表长度的信息并用calculate_z_serial_purepython计算output列表。最后我们对output的内容求和并断言它符合预期的输出值。Ian靠它来确保本书此处不至于出现错误。

既然代码已经确定，我们只需要对所有计算出的值求和就可以验证函数是否如我们预期的那样工作。这叫完整性检查——当我们对计算代码进行修改时，很有必要进行这样的检查来确保我们没有破坏算法的正确性。理想情况下，我们会使用单元测试且对该问题的多个不同配置进行测试。

接下来，在例2-3中，我们定义calculate_z_serial_purepython函数，它扩展了我们之前讨论的算法。注意，我们同时还在函数开头定义了一个output列表，其长度跟zs和cs列表相同。你可能会奇怪为什么我们使用range而不是xrange——这里就先这样，在2.9节，我们会告诉你range有多浪费。

例2-3　我们的CPU密集型计算函数

def calculate_z_serial_purepython(maxiter, zs, cs):
    """Calculate output list using Julia update rule"""
    output = [0] * len(zs)
    for i in range(len(zs)):
        n = 0
        z = zs[i]
        c = cs[i]
        while abs(z) < 2 and n < maxiter:
            z = z * z + c
            n += 1
        output[i] = n
    return output

现在我们可以在例2-4中调用整个计算逻辑。只要将它包入一个__main__检查块，我们就可以安全地为某些性能分析手段导入这个模块而不是直接开始计算。注意，将输出画成图像的方法我们这里没有显示出来。

例2-4　我们代码的main函数

if __name__ == "__main__":
    # Calculate the Julia set using a pure Python solution with
    # reasonable defaults for a laptop
    calc_pure_python(desired_width=1000, max_iterations=300)

一旦我们运行这段代码，就能够看到一些关于问题复杂度的输出：

# running the above produces:
Length of x: 1000
Total elements: 1000000
calculate_z_serial_purepython took 12.3479790688 seconds

图2-1所示的伪灰阶图中，高对比的颜色变化给我们一个该函数运行快慢的直观感受。在此处的图2-3中，我们有一个线性的灰阶图：黑色是算得快的地方而白色是计算开销大的地方。通过对同一数据的两种表达，我们能够看到线性图上丢失了很多细节。有时，当我们调查一个函数的开销时，多种表达方式有助于我们查清问题。