Python-高性能编程-全-

Python 高性能编程（全）

原文：zh.annas-archive.org/md5/7289a66790cbfec1c8bb1a2cb27941dd

译者：飞龙

协议：CC BY-NC-SA 4.0

前言

Python 是一种以其简单性、优雅性和出色的社区支持而闻名的编程语言。得益于大量高质量第三方库的令人印象深刻，Python 被用于许多领域。

在性能关键的应用程序中，通常首选低级语言，如 C、C++和 Fortran。用这些语言编写的程序性能极好，但编写和维护起来都很困难。

Python 是一种更容易处理的编程语言，可以用来快速编写复杂的应用程序。由于其与 C 语言的紧密集成，Python 能够避免动态语言相关的性能下降。你可以使用快速的 C 扩展来编写性能关键代码，同时保留 Python 在其余应用程序中的所有便利性。

在这本书中，你将学习如何逐步使用基本和高级技术找到并加速程序中的慢速部分。

这本书的风格是实用的；每个概念都通过示例进行解释和说明。这本书还讨论了常见的错误，并教授如何避免它们。本书使用的工具相当流行且经过实战检验；你可以确信它们在未来将保持相关性和良好的支持。

这本书从基础知识开始，在此基础上构建，因此，我建议你按顺序阅读章节。

不要忘了享受乐趣！

这本书涵盖的内容

第一章, 基准测试和性能分析 展示了你如何找到程序中需要优化的部分。我们将使用不同用例的工具，并解释如何分析和解释性能分析统计数据。

第二章, 使用 NumPy 进行快速数组操作 是关于 NumPy 包的指南。NumPy 是 Python 中进行数组计算的框架。它提供了一个干净且简洁的 API，以及高效的数组操作。

第三章, 使用 Cython 进行 C 性能优化 是关于 Cython 的教程：一种作为 Python 和 C 之间桥梁的语言。Cython 可以使用 Python 语法的超集来编写代码，并将其编译为高效的 C 扩展。

第四章, 并行处理 是对并行编程的介绍。在本章中，你将学习并行编程与串行编程的不同之处，以及如何并行化简单问题。我们还将解释如何使用多进程、IPython.parallel 和 cython.parallel 为多核编写代码。

这本书你需要的东西

这本书需要安装 Python。除非另有说明，否则示例适用于 Python 2.7 和 Python 3.3。

在这本书中，我们将利用一些流行的 Python 包：

• NumPy（版本 1.7.1 或更高版本）：此软件包可以从官方网站下载（www.scipy.org/scipylib/download.html），并在大多数 Linux 发行版中可用

• Cython（版本 0.19.1 或更高版本）：安装说明可在官方网站上找到（docs.cython.org/src/quickstart/install.html）；请注意，您还需要一个 C 编译器，例如 GCC（GNU 编译器集合），来编译您的 C 扩展

• IPython（版本 0.13.2 或更高版本）：安装说明可在官方网站上找到（ipython.org/install.html）

本书是在 Ubuntu 13.10 操作系统上编写和测试的。示例代码很可能在 Mac OS X 上运行，只需进行少量或无需修改。

我对 Windows 用户的建议是安装 Anaconda Python 发行版（store.continuum.io/cshop/anaconda/），它包含适合科学编程的完整环境。

一个方便的替代方案是使用免费的wakari.io服务：一个基于云的 Linux 和 Python 环境，包括所需的软件包及其工具和实用程序。无需设置。

在第一章中，“基准测试和性能分析”，我们将使用 KCachegrind (sourceforge.net/projects/kcachegrind/)，它适用于 Linux。KCachegrind 还有一个 Windows 版本——QcacheGrind，它也可以从源代码在 Mac OS X 上安装。

本书的目标读者

本书面向中级到高级的 Python 程序员，他们开发的是性能关键型应用程序。由于大多数示例都来自科学应用，本书非常适合希望加快其数值代码的科学家和工程师。

然而，本书的范围很广，概念可以应用于任何领域。由于本书涵盖了基本和高级主题，因此包含了不同 Python 熟练程度程序员的实用信息。

习惯用法

在本书中，您将找到多种文本样式，用于区分不同类型的信息。以下是一些这些样式的示例及其含义的解释。

文本中的代码单词、数据库表名、文件夹名、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 用户名都按照以下格式显示：“plot函数包含在matplotlib中，可以将我们的粒子显示在笛卡尔网格上的点，而FuncAnimation类可以动画化我们的粒子随时间的变化。”

代码块按照以下方式设置：

from matplotlib import pyplot as plt
from matplotlib import animation

def visualize(simulator):

    X = [p.x for p in simulator.particles]
    Y = [p.y for p in simulator.particles]

当我们希望引起您对代码块中特定部分的注意时，相关的行或项目将以粗体显示：

In [1]: import purepy
In [2]: %timeit purepy.loop()
100 loops, best of 3: 8.26 ms per loop
In [3]: %timeit purepy.comprehension()
100 loops, best of 3: 5.39 ms per loop
In [4]: %timeit purepy.generator()
100 loops, best of 3: 5.07 ms per loop

任何命令行输入或输出都按照以下格式编写：

$ time python simul.py # Performance Tuned
real    0m0.756s
user    0m0.714s
sys    0m0.036s

新术语和重要词汇将以粗体显示。您在屏幕上、菜单或对话框中看到的单词，例如，将以文本中的这种形式出现：“您可以通过双击矩形来导航到调用图或调用者地图选项卡。”

注意

警告或重要注意事项将以这样的框显示。

小贴士

小贴士和技巧将以这样的形式出现。

读者反馈

我们始终欢迎读者的反馈。让我们知道您对这本书的看法——您喜欢什么或可能不喜欢什么。读者反馈对我们开发您真正从中受益的标题非常重要。

如要向我们发送一般反馈，请发送电子邮件至<feedback@packtpub.com>，并在邮件主题中提及书籍标题。

如果你在某个领域有专业知识，并且对撰写或为书籍做出贡献感兴趣，请参阅我们的作者指南，网址为www.packtpub.com/authors。

客户支持

现在您已经是 Packt 书籍的骄傲拥有者，我们有一些事情可以帮助您充分利用您的购买。

下载示例代码

您可以从您在www.packtpub.com的账户下载所有已购买的 Packt 书籍的示例代码文件。如果您在其他地方购买了这本书，您可以访问www.packtpub.com/support，并注册以将文件直接通过电子邮件发送给您。

错误清单

尽管我们已经尽一切努力确保我们内容的准确性，但错误仍然可能发生。如果您在我们的某本书中发现错误——可能是文本或代码中的错误——如果您能向我们报告这一点，我们将不胜感激。通过这样做，您可以节省其他读者的挫败感，并帮助我们改进此书的后续版本。如果您发现任何错误清单，请通过访问www.packtpub.com/submit-errata，选择您的书籍，点击错误清单提交表单链接，并输入您的错误清单详情。一旦您的错误清单得到验证，您的提交将被接受，错误清单将被上传到我们的网站，或添加到该标题的错误清单部分。任何现有的错误清单都可以通过从www.packtpub.com/support选择您的标题来查看。

侵权

互联网上版权材料的侵权是一个跨所有媒体的持续问题。在 Packt，我们非常重视我们版权和许可证的保护。如果您在互联网上发现我们作品的任何非法副本，无论形式如何，请立即向我们提供位置地址或网站名称，以便我们可以追究补救措施。

请通过<copyright@packtpub.com>与我们联系，并提供涉嫌侵权材料的链接。

我们感谢您帮助我们保护作者，并确保我们能够为您提供有价值的内容。

询问

如果你在本书的任何方面遇到问题，可以通过<questions@packtpub.com>联系我们，我们将尽力解决。

第一章：基准测试和分析

在加快代码速度时，识别程序中的慢速部分是最重要的任务。在大多数情况下，瓶颈只占程序的一小部分。通过专门解决这些关键点，你可以专注于需要改进的部分，而无需在微优化上浪费时间。

分析是允许我们定位瓶颈的技术。分析器是一个运行代码并观察每个函数运行所需时间的程序，检测程序中的慢速部分。Python 提供了几个工具来帮助我们找到这些瓶颈并导航性能指标。在本章中，我们将学习如何使用标准cProfile模块、line_profiler和memory_profiler。我们还将学习如何使用程序KCachegrind来解释分析结果。

你可能还想评估你程序的总体执行时间，看看你的更改是如何影响的。我们将学习如何编写基准测试以及如何准确计时你的程序。

设计你的应用程序

当你设计一个性能密集型程序时，第一步是编写代码时不要考虑优化；引用唐纳德·克努特的话：

过早的优化是万恶之源。

在早期开发阶段，程序的设计可能会迅速变化，需要你重写和重新组织大量代码。通过在不考虑优化的情况下测试不同的原型，你可以更多地了解你的程序，这将帮助你做出更好的设计决策。

在优化代码时你应该记住的咒语如下：

• 让它运行：我们必须让软件处于工作状态，并确保它产生正确的结果。这一阶段旨在探索我们试图解决的问题，并在早期阶段发现主要的设计问题。

• 让它正确：我们想确保程序的设计是稳固的。在尝试任何性能优化之前应该进行重构。这实际上有助于将应用程序分解成独立且一致的单元，这些单元更容易维护。

• 让它变得更快：一旦我们的程序运行良好并且有良好的设计，我们希望优化程序中不够快的部分。如果内存使用构成问题，我们可能还想优化内存使用。

在本节中，我们将分析一个测试应用程序——粒子模拟器。模拟器是一个程序，它根据我们将建立的一系列法则，随着时间的推移演变一些粒子。这些粒子可以是抽象实体，也可以对应于物理对象。例如，它们可以是桌子上移动的台球、气体中的分子、在空间中移动的恒星、烟雾粒子、室内的流体等等。

这些模拟在物理学、化学和天文学等领域非常有用，用于模拟物理系统的程序通常是性能密集型的。为了研究现实系统，通常需要模拟尽可能多的物体。

在我们的第一个例子中，我们将模拟一个包含粒子围绕中心点以不同速度不断旋转的系统，就像时钟的指针一样。

运行我们的模拟所需的信息将是粒子的起始位置、速度和旋转方向。从这些元素中，我们必须计算粒子在下一个时间瞬间的位置。

圆周运动的基本特征是粒子始终垂直于连接粒子和中心的连线方向移动，如前图所示。为了移动粒子，我们只需通过在运动方向上采取一系列非常小的步骤来改变其位置，如下图所示：

我们将首先以面向对象的方式设计应用程序。根据我们的要求，自然有一个通用的 Particle 类，它简单地存储粒子的位置 (x，y) 和其角速度：

class Particle:
    def __init__(self, x, y, ang_speed):
        self.x = x
        self.y = y
        self.ang_speed = ang_speed

另一个类，称为 ParticleSimulator，将封装我们的运动定律，并负责随时间改变粒子的位置。__init__ 方法将存储 Particle 实例的列表，而 evolve 方法将根据我们的定律改变粒子的位置。

我们希望粒子围绕点 (x，y) 旋转，在这里，它等于 (0, 0)，以恒定速度。粒子的方向始终垂直于从中心（参见图章的第一幅图）的方向。为了找到这个向量

（对应于 Python 变量 v_x 和 v_y）使用以下公式就足够了：

如果我们让其中一个粒子移动，经过一定的时间 dt，它将沿着圆形路径移动，到达另一个位置。为了让粒子沿着这条轨迹移动，我们必须将时间间隔 dt 分成非常小的时间步长，其中粒子沿着圆周切线方向移动。最终结果只是圆形运动的近似，实际上它类似于多边形。时间步长应该非常小，否则粒子的轨迹会迅速发散，如下图所示：

以更简化的方式，为了计算时间 dt 时的粒子位置，我们必须执行以下步骤：

1. 计算运动方向：v_x，v_y。

2. 计算位移 (d_x，d_y)，它是时间和速度的乘积，并遵循运动方向。

3. 重复步骤 1 和 2 足够长的时间步，以覆盖总时间dt。

以下代码显示了完整的ParticleSimulator实现：

class ParticleSimulator:

    def __init__(self, particles):
        self.particles = particles

    def evolve(self, dt):
        timestep = 0.00001
        nsteps = int(dt/timestep)

        for i in range(nsteps):
            for p in self.particles:

 # 1\. calculate the direction
 norm = (p.x**2 + p.y**2)**0.5
 v_x = (-p.y)/norm
 v_y = p.x/norm

 # 2\. calculate the displacement
 d_x = timestep * p.ang_speed * v_x
 d_y = timestep * p.ang_speed * v_y

 p.x += d_x
 p.y += d_y
 # 3\. repeat for all the time steps

我们可以使用matplotlib库来可视化我们的粒子。这个库不包括在 Python 标准库中。要安装它，您可以按照官方文档中的说明进行操作：

matplotlib.org/users/installing.html

注意

或者，您可以使用包含matplotlib和本书中使用的其他大多数第三方包的 Anaconda Python 发行版（store.continuum.io/cshop/anaconda/)。Anaconda 是免费的，适用于 Linux、Windows 和 Mac。

matplotlib中的plot函数可以将我们的粒子显示在笛卡尔网格上的点，而FuncAnimation类可以动画化粒子随时间的变化。

visualize函数通过获取粒子模拟器并在动画图中显示轨迹来实现这一点。

visualize函数的结构如下：

• 设置坐标轴并使用绘图函数将粒子显示为点

• 编写一个初始化函数（init）和一个更新函数（animate），使用line.set_data方法改变数据点的x，y坐标

• 创建一个FuncAnimation实例，传递函数和一些参数

• 使用plt.show()运行动画

  可视化函数的完整实现如下：

  from matplotlib import pyplot as plt
  from matplotlib import animation
  
  def visualize(simulator):
  
      X = [p.x for p in simulator.particles]
      Y = [p.y for p in simulator.particles]
  
      fig = plt.figure()
      ax = plt.subplot(111, aspect='equal')
      line, = ax.plot(X, Y, 'ro')
  
      # Axis limits
      plt.xlim(-1, 1)
      plt.ylim(-1, 1)
  
      # It will be run when the animation starts
      def init():
          line.set_data([], [])
          return line,
  
      def animate(i):
          # We let the particle evolve for 0.1 time units
          simulator.evolve(0.01)
          X = [p.x for p in simulator.particles]
          Y = [p.y for p in simulator.particles]
  
          line.set_data(X, Y)
          return line,
  
      # Call the animate function each 10 ms
      anim = animation.FuncAnimation(fig, animate, init_func=init, blit=True,# Efficient animation
                                     interval=10)
      plt.show()
  

最后，我们定义了一个小的测试函数——test_visualize——它动画化了一个由三个不同方向旋转的粒子组成的系统。注意，第三个粒子比其他粒子快三倍完成一圈：

def test_visualize():
    particles = [Particle( 0.3,  0.5, +1),
                 Particle( 0.0, -0.5, -1),
                 Particle(-0.1, -0.4, +3)]

    simulator = ParticleSimulator(particles)
    visualize(simulator)

if __name__ == '__main__':
    test_visualize()

编写测试和基准

现在我们有一个工作的模拟器，我们可以开始测量我们的性能并调整我们的代码，以便我们的模拟器可以处理尽可能多的粒子。这个过程的第一步是编写一个测试和一个基准。

我们需要一个测试来检查模拟产生的结果是否正确。在优化过程中，我们将重写代码以尝试不同的解决方案；这样做可能会轻易引入错误。维护一个坚实的测试套件对于避免在损坏的代码上浪费时间至关重要。

我们的测试将使用三个粒子，并让系统演化 0.1 个时间单位。然后，我们将我们的结果与参考实现的结果进行比较，直到一定的精度：

def test():
    particles = [Particle( 0.3,  0.5, +1),
                 Particle( 0.0, -0.5, -1),
                 Particle(-0.1, -0.4, +3)]

    simulator = ParticleSimulator(particles)

    simulator.evolve(0.1)

    p0, p1, p2 = particles

    def fequal(a, b):
        return abs(a - b) < 1e-5

    assert fequal(p0.x, 0.2102698450356825)
    assert fequal(p0.y, 0.5438635787296997)

    assert fequal(p1.x, -0.0993347660567358)
    assert fequal(p1.y, -0.4900342888538049)

    assert fequal(p2.x,  0.1913585038252641)
    assert fequal(p2.y, -0.3652272210744360)

if __name__ == '__main__':
    test()

我们还希望编写一个基准，可以测量我们应用程序的性能。这将提供有关我们相对于先前实现的改进程度的指示。

在我们的基准测试中，我们实例化了 100 个具有随机坐标和角速度的Particle对象，并将它们输入到ParticleSimulator类中。然后，我们让系统演化 0.1 个时间单位：

from random import uniform

def benchmark():
    particles = [Particle(uniform(-1.0, 1.0),
                          uniform(-1.0, 1.0),
                          uniform(-1.0, 1.0))
                  for i in range(1000)]

    simulator = ParticleSimulator(particles)
    simulator.evolve(0.1)

if __name__ == '__main__':
    benchmark()

测试基准的时间

你可以通过使用 Unix 的 time 命令从命令行轻松测量任何进程的执行时间：

$ time python simul.py
real    0m1.051s
user    0m1.022s
sys    0m0.028s

注意

time 命令在 Windows 上不可用，但可以在从官方网站 www.cygwin.com/ 下载的 cygwin 命令行界面中找到。

默认情况下，time 命令显示三个指标：

• real：从开始到结束运行进程的实际时间，就像用秒表测量的人一样

• user：所有 CPU 在计算过程中累计花费的时间

• sys：所有 CPU 在与系统相关的任务（如内存分配）中累计花费的时间

注意，有时 user + sys 可能会大于 real，因为多个处理器可能并行工作。

提示

time 还提供了几个格式化选项；要了解概述，你可以查看其手册（使用 man time 命令）。如果你想要所有可用指标的总览，可以使用 -v 选项。

Unix 的 time 命令是衡量你程序的好方法。为了获得更精确的测量，基准测试应该运行足够长的时间（以秒为单位），使得过程的设置和拆除与执行时间相比变得很小。user 指标适合作为 CPU 性能的监控器，因为 real 指标还包括在其他进程或等待 I/O 操作中花费的时间。

另一个用于计时 Python 脚本的有用程序是 timeit 模块。该模块将代码片段在循环中运行 n 次，并测量所花费的时间。然后，它重复此操作 r 次（默认情况下，r 的值为 3）并取最佳运行结果。由于此过程，timeit 适用于准确计时独立的小语句。

timeit 模块可以用作 Python 模块，从命令行或从 IPython 中使用。

IPython 是一个专为交互式使用设计的 Python 命令行界面。它增强了 tab 自动完成和许多时间、性能分析和调试代码的实用工具。我们将利用这个 shell 在整本书中尝试代码片段。IPython 命令行界面接受 魔法命令——以 % 符号开头的语句，这些语句增强了 shell 的特殊行为。以 %% 开头的命令称为 单元魔法命令，这些命令可以应用于多行代码片段（称为 单元）。

IPython 在大多数 Linux 发行版中可用，并包含在 Anaconda 中。你可以按照官方文档中的安装说明进行操作：

ipython.org/install.html

提示

你可以将 IPython 用作常规的 Python 命令行界面（ipython），但它也提供基于 Qt 的版本（ipython qtconsole）以及强大的基于浏览器的界面（ipython notebook）。

在 IPython 和命令行接口中，可以使用 -n 和 -r 选项指定循环或重复的次数，否则它们将自动确定。当从命令行调用 timeit 时，你也可以提供一个在循环执行语句之前运行的设置代码。

在以下代码中，我们展示了如何从 IPython、命令行以及作为 Python 模块使用 timeit：

# IPython Interface
$ ipython
In [1]: from simul import benchmark
In [2]: %timeit benchmark()
1 loops, best of 3: 782 ms per loop

# Command Line Interface
$ python -m timeit -s 'from simul import benchmark' 'benchmark()'10 loops, best of 3: 826 msec per loop

# Python Interface
# put this function into the simul.py script

import timeit
result = timeit.timeit('benchmark()',
                                   setup='from __main__ import benchmark', number=10)
# result is the time (in seconds) to run the whole loop

result = timeit.repeat('benchmark()', setup='from __main__ import benchmark', number=10, repeat=3)
# result is a list containing the time of each repetition (repeat=3 in this case)

注意，虽然命令行和 IPython 接口会自动确定一个合理的 n 值，但 Python 接口需要你显式地将其作为 number 参数传递。

使用 cProfile 寻找瓶颈

在评估程序的执行时间后，我们准备识别需要性能调整的代码部分。这些部分通常与程序的大小相比相当小。

历史上，Python 的标准库中有三个不同的配置文件模块：

• profile 模块：此模块是用纯 Python 编写的，并为程序执行添加了显著的开销。它在标准库中的存在主要是由于其可扩展性。

• hotshot 模块：一个旨在最小化配置文件开销的 C 模块。Python 社区不建议使用它，并且在 Python 3 中不可用。

• cProfile 模块：主要的配置文件模块，其接口类似于 profile。它具有较小的开销，适合作为通用配置文件。

我们将展示如何以两种不同的方式使用 cProfile 模块：

• 从命令行

• 从 IPython

为了使用 cProfile，不需要对代码进行任何更改，可以直接在现有的 Python 脚本或函数上执行。

你可以使用以下方式从命令行使用 cProfile：

$ python -m cProfile simul.py

这将打印出包含多个配置文件指标的冗长输出。你可以使用 -s 选项按某个指标排序输出：

$ python -m cProfile -s tottime simul.py

你可以通过传递 -o 选项将输出文件保存为 stats 模块和其他工具可读的格式：

$ python -m cProfile -o prof.out simul.py

你还可以从 IPython 中交互式地配置文件。%prun 魔法命令允许你使用 cProfile 配置文件：

In [1]: from simul import benchmark
In [2]: %prun benchmark()
         707 function calls in 0.793 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.792    0.792    0.792    0.792 simul.py:12(evolve)
        1    0.000    0.000    0.000    0.000 simul.py:100(<listcomp>)
      300    0.000    0.000    0.000    0.000 random.py:331(uniform)
      100    0.000    0.000    0.000    0.000 simul.py:2(__init__)
        1    0.000    0.000    0.793    0.793 {built-in method exec}
      300    0.000    0.000    0.000    0.000 {method 'random' of '_random.Random' objects}
        1    0.000    0.000    0.793    0.793 simul.py:99(benchmark)
        1    0.000    0.000    0.793    0.793 <string>:1(<module>)
        1    0.000    0.000    0.000    0.000 simul.py:9(__init__)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

cProfile 的输出分为五列：

• ncalls：函数被调用的次数。

• tottime：在函数中花费的总时间，不考虑对其他函数的调用。

• cumtime：在函数中花费的时间，包括其他函数调用。

• percall：函数单次调用的花费时间——可以通过将总时间或累积时间除以调用次数来获得。

• filename:lineno：文件名和相应的行号。当调用 C 扩展模块时，此信息不存在。

最重要的指标是 tottime，即函数体中实际花费的时间，不包括子调用。在我们的案例中，大部分时间都花在了 evolve 函数上。我们可以想象，循环是代码中需要性能调整的部分。

对于具有大量调用和子调用的程序，以文本方式分析数据可能会令人望而却步。一些图形工具通过改进交互式界面来辅助任务。

KCachegrind 是一个 GUI（图形用户界面），用于分析不同程序的剖析输出。

注意

KCachegrind 可在 Ubuntu 13.10 官方仓库中找到。Qt 端口，QCacheGrind 可以从以下网页下载到 Windows：

sourceforge.net/projects/qcachegrindwin/

Mac 用户可以通过 Mac Ports (www.macports.org/)编译 QCacheGrind，按照此链接博客文章中的说明进行操作：

blogs.perl.org/users/rurban/2013/04/install-kachegrind-on-macosx-with-ports.html

KCachegrind 不能直接读取由cProfile产生的输出文件。幸运的是，pyprof2calltree第三方 Python 模块能够将cProfile输出文件转换为 KCachegrind 可读的格式。

您可以从源代码安装pyprof2calltree(pypi.python.org/pypi/pyprof2calltree/)或从 Python 包索引(pypi.python.org/)安装。

为了最好地展示 KCachegrind 的功能，我们将使用一个具有更多样化结构的另一个示例。我们定义一个递归函数factorial，以及两个使用factorial的其他函数，它们是taylor_exp和taylor_sin。它们代表exp(x)和sin(x)的泰勒近似的多项式系数：

def factorial(n):
    if n == 0:
        return 1.0
    else:
        return float(n) * factorial(n-1)

def taylor_exp(n):
    return [1.0/factorial(i) for i in range(n)]

def taylor_sin(n):
    res = []
    for i in range(n):
        if i % 2 == 1:
           res.append((-1)**((i-1)/2)/float(factorial(i)))
        else:
           res.append(0.0)
    return res

def benchmark():
    taylor_exp(500)
    taylor_sin(500)

if __name__ == '__main__':
    benchmark()

我们需要首先生成cProfile输出文件：

$ python -m cProfile -o prof.out taylor.py

然后，我们可以使用pyprof2calltree转换输出文件并启动 KCachegrind：

$ pyprof2calltree -i prof.out -o prof.calltree
$ kcachegrind prof.calltree # or qcachegrind prof.calltree

上一张图片是 KCachegrind 用户界面的截图。在左侧，我们有一个与cProfile相当输出的结果。实际的列名略有不同：Incl.对应于cProfile模块的cumtime；Self对应于tottime。通过在菜单栏上点击Relative按钮，数值以百分比形式给出。通过点击列标题，您可以按相应属性排序。

在右上角，点击Callee Map标签包含函数成本的图表。在图表中，每个函数由一个矩形表示，函数花费的时间百分比与矩形的面积成正比。矩形可以包含表示对其他函数子调用的子矩形。在这种情况下，我们可以很容易地看到有两个矩形表示factorial函数。左侧的对应于taylor_exp的调用，右侧的对应于taylor_sin的调用。

在右下角，您可以通过点击调用图选项卡来显示另一个图表——调用图。调用图是函数之间调用关系的图形表示：每个方块代表一个函数，箭头表示调用关系。例如，taylor_exp调用（一个列表推导），它调用factorial 500 次；taylor_sin调用factorial 250 次。KCachegrind 还可以检测递归调用：factorial调用自身 187250 次。

您可以通过双击矩形导航到调用图或调用者映射选项卡；界面将相应更新，显示时间属性相对于所选函数。例如，双击taylor_exp将导致图表更改，仅显示taylor_exp对总成本的贡献。

注意

Gprof2Dot (code.google.com/p/jrfonseca/wiki/Gprof2Dot) 是另一个流行的工具，用于生成调用图。从支持的剖析器产生的输出文件开始，它将生成一个.dot图，表示调用图。

使用 line_profiler 逐行分析

现在我们知道了哪个函数需要优化，我们可以使用line_profiler模块，它以逐行的方式显示时间消耗。这在难以确定哪些语句成本高昂的情况下非常有用。line_profiler模块是一个第三方模块，可在 Python 包索引上找到，并可通过其网站上的说明进行安装：

pythonhosted.org/line_profiler/

为了使用line_profiler，我们需要将@profile装饰器应用到我们打算监控的函数上。请注意，您不需要从另一个模块导入profile函数，因为它在运行分析脚本kernprof.py时会被注入到全局命名空间中。为了生成程序的配置文件输出，我们需要将@profile装饰器添加到evolve函数上：

@profile
def evolve:
    # code

脚本kernprof.py将生成一个输出文件，并将分析结果打印到标准输出。我们应该使用两个选项来运行脚本：

• -l 选项用于使用line_profiler函数

• -v 选项用于立即在屏幕上打印结果

  $ kernprof.py -l -v simul.py
  
  

还可以在 IPython shell 中运行剖析器以进行交互式编辑。您应该首先加载line_profiler扩展，这将提供lprun魔法命令。通过使用该命令，您可以避免添加@profile装饰器。

In [1]: %load_ext line_profiler
In [2]: from simul import benchmark, ParticleSimulator
In [3]: %lprun -f ParticleSimulator.evolve benchmark()

Timer unit: 1e-06 s

File: simul.py
Function: evolve at line 12
Total time: 5.31684 s

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    12                                               def evolve(self, dt):
    13         1            9      9.0      0.0          timestep = 0.00001
    14         1            4      4.0      0.0          nsteps = int(dt/timestep)
    15                                                   
    16     10001         5837      0.6      0.1          for i in range(nsteps):
    17   1010000       517504      0.5      9.7              for p 
in self.particles:
    18                                           
    19   1000000       963498      1.0     18.1                  norm = (p.x**2 + p.y**2)**0.5
    20   1000000       621063      0.6     11.7                  v_x = (-p.y)/norm
    21   1000000       577882      0.6     10.9                  v_y = p.x/norm
    22                                                           
    23   1000000       672811      0.7     12.7                  d_x = timestep * p.ang_speed * v_x
    24   1000000       685092      0.7     12.9                  d_y = timestep * p.ang_speed * v_y
    25                                           
    26   1000000       650802      0.7     12.2                  p.x += d_x
    27   1000000       622337      0.6     11.7                  p.y += d_y

输出非常直观，分为几列：

• 行号：运行的行号

• 命中次数：该行运行的次数

• 时间：该行的执行时间，以微秒为单位（时间）

• 每次命中：时间除以命中次数

• % 时间：执行该行所花费的总时间的比例

• 行内容：相应行的源代码

通过查看百分比列，我们可以大致了解时间花费在哪里。在这种情况下，for循环体中有几个语句，每个语句的成本大约为 10-20%。

优化我们的代码

现在我们已经确切地知道了时间是如何花费的，我们可以修改代码并评估性能的变化。

有几种不同的方法可以调整我们的纯 Python 代码。通常产生最显著结果的方法是改变算法。在这种情况下，与其计算速度并添加小步骤，不如用半径r和角度alpha（而不是x和y）表示运动方程，然后使用以下方程计算圆上的点：

x = r * cos(alpha)
y = r * sin(alpha)

另一种方法是通过最小化指令的数量。例如，我们可以预先计算因子timestep * p.ang_speed，这个因子不会随时间变化。我们可以交换循环顺序（首先迭代粒子，然后迭代时间步），并将因子的计算放在循环外的粒子中。

行内分析还显示，即使是简单的赋值操作也可能花费相当多的时间。例如，以下语句占用了超过 10%的总时间：

 v_x = (-p.y)/norm

因此，优化循环的一种方法是通过减少赋值操作的数量。为了做到这一点，我们可以牺牲可读性，通过将表达式重写为单个稍微复杂一些的语句来避免中间变量（注意，右侧在赋值给变量之前被完全评估）：

p.x, p.y = p.x - t_x_ang*p.y/norm, p.y + t_x_ang * p.x/norm

这导致以下代码：

    def evolve_fast(self, dt):
        timestep = 0.00001
        nsteps = int(dt/timestep)

        # Loop order is changed
        for p in self.particles:
            t_x_ang = timestep * p.ang_speed
            for i in range(nsteps):
                norm = (p.x**2 + p.y**2)**0.5
                p.x, p.y = (p.x - t_x_ang *  p.y/norm,p.y + t_x_ang * p.x/norm)

应用更改后，我们应该通过运行我们的测试来确保结果仍然相同。然后我们可以使用我们的基准来比较执行时间：

$ time python simul.py # Performance Tuned
real    0m0.756s
user    0m0.714s
sys    0m0.036s

$ time python simul.py # Original
real    0m0.863s
user    0m0.831s
sys    0m0.028s

通过对纯 Python 进行操作，我们只获得了速度的微小提升。

dis模块

有时候，很难评估 Python 语句将执行多少操作。在本节中，我们将探索 Python 内部机制来估计 Python 语句的性能。Python 代码被转换为中间表示——称为字节码——由 Python 虚拟机执行。

为了帮助检查代码如何转换为字节码，我们可以使用 Python 模块dis（反汇编）。它的使用非常简单，只需在ParticleSimulator.evolve方法上调用函数dis.dis即可：

import dis
from simul import ParticleSimulator
dis.dis(ParticleSimulator.evolve)

这将为每一行生成一个字节码指令列表。例如，语句v_x = (-p.y)/norm在以下指令集中展开：

20         85 LOAD_FAST                5 (p)
             88 LOAD_ATTR                4 (y)
             91 UNARY_NEGATIVE       
             92 LOAD_FAST                6 (norm)
             95 BINARY_TRUE_DIVIDE   
             96 STORE_FAST               7 (v_x)

LOAD_FAST将变量p的引用加载到栈上，LOAD_ATTR加载栈顶元素的y属性。其他指令（UNARY_NEGATIVE和BINARY_TRUE_DIVIDE）在栈顶元素上执行算术运算。最后，结果存储在v_x中（STORE_FAST）。

通过分析完整的 dis 输出，我们可以看到循环的第一个版本产生了 51 条字节码指令，而第二个版本被转换成了 35 条指令。

dis 模块有助于发现语句是如何被转换的，并主要作为探索和学习 Python 字节码表示的工具。

为了进一步提高我们的性能，我们可以继续尝试找出其他方法来减少指令的数量。然而，很明显，这种方法有一些限制，可能不是完成这项工作的正确工具。在下一章中，我们将看到如何借助 NumPy 加快这类计算。

使用 memory_profiler 分析内存使用情况

在某些情况下，内存使用可能成为一个问题。例如，如果我们想要处理大量的粒子，由于创建了大量的 Particle 实例，我们将会有内存开销。

memory_profiler 模块以一种类似于 line_profiler 的方式总结进程的内存使用情况。

注意

memory_profiler 包也可在 Python 包索引中找到。你还应该安装 psutil 模块（code.google.com/p/psutil/）作为可选依赖项，这将使 memory_profiler 运行得更快。

就像 line_profiler 一样，memory_profiler 也需要源代码的仪器化，通过在我们要监控的函数上放置 @profile 装饰器。在我们的例子中，我们想要分析 benchmark 函数。

我们可以稍微修改 benchmark 以实例化大量的 Particle 实例（100000 个）并减少模拟时间：

def benchmark_memory():
    particles = [Particle(uniform(-1.0, 1.0),
                          uniform(-1.0, 1.0),
                          uniform(-1.0, 1.0))
                  for i in range(100000)]

    simulator = ParticleSimulator(particles)
    simulator.evolve(0.001)

我们可以通过 IPython 壳中的魔法命令 %mprun 使用 memory_profiler：

In [1]: %load_ext memory_profiler
In [2]: from simul import benchmark_memory
In [3]: %mprun -f benchmark_memory benchmark_memory()

Line #    Mem usage    Increment   Line Contents
==============================================
   135     45.5 MiB      0.0 MiB   def benchmark_memory():
   136     45.5 MiB      0.0 MiB       particles = [Particle(uniform(-1.0, 1.0),
   137                                                       uniform(-1.0, 1.0),
   138                                                       uniform(-1.0, 1.0))
 139     71.2 MiB     25.7 MiB                     for i in range(100000)]
   140                                 
   141     71.2 MiB      0.0 MiB       simulator = ParticleSimulator(particles)
   142     71.3 MiB      0.1 MiB       simulator.evolve(0.001)

小贴士

在添加了 @profile 装饰器后，可以使用 mprof run 命令从壳中运行 memory_profiler。

从输出中我们可以看到，100000 个 Particle 对象占用了 25.7 MiB 的内存。

小贴士

1 MiB（梅比字节）相当于 1024² = 1,048,576 字节。它与 1 MB（兆字节）不同，1 MB 相当于 1000² = 1,000,000 字节。

我们可以在 Particle 类上使用 __slots__ 来减少其内存占用。这个特性通过避免在内部字典中存储实例变量来节省一些内存。这种优化的缺点是：它阻止了添加 __slots__ 中未指定的属性（要在 Python 2 中使用此功能，你应该确保你正在使用新式类）：

class Particle: 
# class Particle(object):  # New-style class for Python 2

    __slots__ = ('x', 'y', 'ang_speed')

    def __init__(self, x, y, ang_speed):
        self.x = x
        self.y = y
        self.ang_speed = ang_speedWe can now re-run our benchmark:
In [1]: %load_ext memory_profiler
In [2]: from simul import benchmark_memory
In [3]: %mprun -f benchmark_memory benchmark_memory()

Line #    Mem usage    Increment   Line Contents
==============================================
   138     45.5 MiB      0.0 MiB   def benchmark_memory():
   139     45.5 MiB      0.0 MiB       particles = [Particle(uniform(-1.0, 1.0),
   140                                                       uniform(-1.0, 1.0),
   141                                                       uniform(-1.0, 1.0))
 142     60.2 MiB     14.7 MiB                     for i in range(100000)]
   143                                 
   144     60.2 MiB      0.0 MiB       simulator = ParticleSimulator(particles)
   145     60.3 MiB      0.1 MiB       simulator.evolve(0.001)

通过使用 __slots__ 重新编写 Particle 类，我们可以节省 11 MiB 的内存。

纯 Python 代码的性能调优技巧

作为一条经验法则，在优化纯 Python 代码时，你应该查看标准库中可用的内容。标准库包含针对最常见的数据结构（如列表、字典和集合）的巧妙算法。此外，许多标准库模块是用 C 实现的，具有快速的执行时间。然而，始终计时不同的解决方案是很重要的——结果往往是不可预测的。

collections模块提供了额外的数据容器，可以高效地处理一些常见操作。例如，当你需要从开始处弹出项目并在末尾追加新项目时，你可以使用deque代替列表。collections模块还包括一个Counter类，可以用来计算可迭代对象中的重复元素。请注意，Counter可能比用标准循环在字典中编写的等效代码要慢：

def counter_1():
    items = [random.randint(0, 10) for i in range(10000)]
    return Counter(items)

def counter_2():
    items = [random.randint(0, 10) for i in range(10000)]
    counter = {}
    for item in items:
        if item not in counter:
            counter[item] = 0
        else:
            counter[item] += 1
    return counter

你可以将代码放入一个名为purepy.py的文件中，并通过 IPython 进行计时：

In [1]: import purepy
In [2]: %timeit purepy.counter_1()
100 loops, best of 3: 10.1 ms per loop
In [3]: %timeit purepy.counter_2()
100 loops, best of 3: 9.11 ms per loop

通常，应该优先使用列表推导和生成器，而不是显式循环。即使与标准循环相比速度提升不大，这也是一个好习惯，因为它提高了可读性。在以下示例中，我们可以看到，列表推导和生成器表达式与函数sum结合时比显式循环更快：

def loop():
    res = []
    for i in range(100000):
        res.append(i * i)
    return sum(res)

def comprehension():
    return sum([i * i for i in range(100000)])

def generator():
    return sum(i * i for i in range(100000))

我们可以将这些函数添加到purepy.py中，并用 IPython 进行测试：

In [1]: import purepy
In [2]: %timeit purepy.loop()
100 loops, best of 3: 8.26 ms per loop
In [3]: %timeit purepy.comprehension()
100 loops, best of 3: 5.39 ms per loop
In [4]: %timeit purepy.generator()
100 loops, best of 3: 5.07 ms per loop

bisect模块可以帮助快速插入和检索元素，同时保持列表排序。

纯 Python 代码的原始优化并不非常有效，除非有实质性的算法优势。加快你代码的第二种最佳方式是使用专门为此目的设计的库，如numpy，或者使用Cython在更接近底层语言（如 C）中编写扩展模块。

摘要

在本章中，我们介绍了优化的基本原理，并将这些原理应用于我们的测试应用。最重要的是在编辑代码之前识别应用中的瓶颈。我们看到了如何使用time Unix 命令和 Python 的timeit模块编写和计时基准测试。我们学习了如何使用cProfile、line_profiler和memory_profiler来分析我们的应用，以及如何使用 KCachegrind 图形化地导航分析数据。我们还概述了一些利用标准库中可用的工具来优化纯 Python 代码的策略。

在下一章中，我们将看到如何使用numpy以简单方便的方式显著加快计算速度。

第二章. 使用 NumPy 进行快速数组操作

NumPy 是 Python 中科学计算的事实标准。它通过提供灵活的多维数组扩展了 Python，允许快速进行数学计算。

NumPy 作为一个框架，允许使用简洁的语法进行复杂操作。多维数组(numpy.ndarray)在内部基于 C 数组：这样，开发者可以轻松地将 NumPy 与现有的 C 和 FORTRAN 代码接口。NumPy 构成了 Python 和这些语言编写的遗留代码之间的桥梁。

在本章中，我们将学习如何创建和操作 NumPy 数组。我们还将探索 NumPy 广播功能，以高效简洁的方式重写复杂的数学表达式。

在过去几年中，开发了许多包以进一步提高 NumPy 的速度。我们将探索这些包中的一个，numexpr，它优化数组表达式并利用多核架构。

开始使用 NumPy

NumPy 围绕其多维数组对象numpy.ndarray建立。NumPy 数组是相同数据类型元素的集合；这种基本限制允许 NumPy 以高效的方式打包数据。通过这种方式存储数据，NumPy 可以以高速处理算术和数学运算。

创建数组

你可以使用numpy.array函数创建 NumPy 数组。它接受一个类似列表的对象（或另一个数组）作为输入，并可选地接受一个表示其数据类型的字符串。你可以使用 IPython shell 交互式测试数组创建，如下所示：

In [1]: import numpy as np
In [2]: a = np.array([0, 1, 2])

每个 NumPy 数组都有一个数据类型，可以通过dtype属性访问，如下面的代码所示。在下面的代码示例中，dtype是一个 64 位整数：

In [3]: a.dtype
Out[3]: dtype('int64')

如果我们想将这些数字视为float类型的变量，我们可以在np.array函数中传递dtype参数，或者使用astype方法将数组转换为另一种数据类型，如下面的代码所示：

In [4]: a = np.array([1, 2, 3], dtype='float32')
In [5]: a.astype('float32')
Out[5]: array([ 0.,  1.,  2.], dtype=float32)

要创建一个二维数组（数组数组），我们可以使用嵌套序列初始化数组，如下所示：

In [6]: a = np.array([[0, 1, 2], [3, 4, 5]])
In [7]: print(a)
Out[7]: [[0 1 2][3 4 5]]

以这种方式创建的数组有两个维度——在 NumPy 术语中称为轴。这样的数组就像一个包含两行三列的表格。我们可以使用ndarray.shape属性访问轴结构：

In [7]: a.shape
Out[7]: (2, 3)

数组也可以被重塑，只要形状维度的乘积等于数组中元素的总数。例如，我们可以以下列方式重塑包含 16 个元素的数组：(2, 8)，(4, 4)，或(2, 2, 4)。要重塑数组，我们可以使用ndarray.reshape方法或直接更改ndarray.shape属性。以下代码展示了ndarray.reshape方法的使用：

In [7]: a = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8,9, 10, 11, 12, 13, 14, 15])
In [7]: a.shape
Out[7]: (16,)
In [8]: a.reshape(4, 4) # Equivalent: a.shape = (4, 4)
Out[8]:
array([[ 0,  1,  2,  3],[ 4,  5,  6,  7],[ 8,  9, 10, 11],[12, 13, 14, 15]])

由于这个特性，你也可以自由地添加大小为 1 的维度。你可以将包含 16 个元素的数组重塑为(16, 1)，(1, 16)，(16, 1, 1)等等。

NumPy 提供了便捷的函数，如下面的代码所示，用于创建填充零的数组、填充一的数组或没有初始化值（empty—它们的实际值没有意义，取决于内存状态）的数组。这些函数接受数组形状作为元组，并可选地接受其 dtype：

In [8]: np.zeros((3, 3))
In [9]: np.empty((3, 3))
In [10]: np.ones((3, 3), dtype='float32')

在我们的例子中，我们将使用 numpy.random 模块在 (0, 1) 区间内生成随机浮点数。在下面的代码中，我们使用 np.random.rand 函数生成一个形状为 (3, 3) 的随机数数组：

In [11]: np.random.rand(3, 3)

有时，初始化与其它数组形状相似的数组很方便。同样，NumPy 提供了一些方便的函数来实现这个目的，如 zeros_like，empty_like 和 ones_like。这些函数如下所示：

In [12]: np.zeros_like(a)
In [13]: np.empty_like(a)
In [14]: np.ones_like(a)

访问数组

NumPy 数组接口在浅层上与 Python 列表相似。它们可以使用整数进行索引，也可以使用 for 循环进行迭代。以下代码展示了如何索引和迭代一个数组：

In [15]: A = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8])
In [16]: A[0]
Out[16]: 0
In [17]: [a for a in A]
Out[17]: [0, 1, 2, 3, 4, 5, 6, 7, 8]

也可以在多维中对数组进行索引。如果我们取一个 (3, 3) 的数组（包含 3 个三元组的数组）并索引第一个元素，我们将获得如下所示的第一个三元组：

In [18]: A = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
In [19]: A[0]
Out[19]: array([0, 1, 2])

我们可以通过添加另一个索引并用逗号分隔来再次索引三元组。要获取第一个三元组的第二个元素，我们可以使用 [0, 1] 进行索引，如下面的代码所示：

In [20]: A[0, 1]
Out[20]: 1

NumPy 允许你在单维和多维中切片数组。如果我们对第一维进行索引，我们将得到如下所示的一组三元组：

In [21]: A[0:2]
Out[21]: array([[0, 1, 2],
               [3, 4, 5]])

如果我们用 [0:2] 切片数组，对于每个选定的三元组，我们提取前两个元素，结果是一个 (2, 2) 的数组，如下面的代码所示：

In [22]: A[0:2, 0:2]
Out[22]: array([[0, 1],
                [3, 4]])

直观地，你可以使用数值索引和切片来更新数组中的值。语法如下：

In [23]: A[0, 1] = 8
In [24]: A[0:2, 0:2] = [[1, 1], [1, 1]]

小贴士

使用切片语法进行索引速度快，因为它不会复制数组。在 NumPy 术语中，它返回对相同内存区域的 视图。如果我们从原始数组中取一个切片并更改其值之一；原始数组也将被更新。以下代码演示了同样的例子：

In [25]: a = np.array([1, 1, 1, 1])
In [26]: a_view = a[0:2]
In [27]: a_view[0] = 2
In [28]: print(a)
Out[28]: [2 1 1 1]

我们可以看看另一个例子，展示如何在现实场景中使用切片语法。我们定义一个数组 r_i，如下面的代码行所示，它包含一组 10 个坐标（x, y）；其形状将是 (10, 2)：

In [29]: r_i = np.random.rand(10, 2)

一个典型的操作是提取每个坐标的 x 分量。换句话说，你想要提取项目 [0, 0]，[1, 0]，[2, 0]，依此类推，结果是一个形状为 (10,) 的数组。有助于思考的是，第一个索引是 移动的，而第二个索引是 固定的（在 0）。有了这个想法，我们将切片第一轴（移动的那个）上的每个索引，并在第二轴（固定的那个）上取第一个元素，如下面的代码行所示：

In [30]: x_i = r_i[:, 0]

另一方面，以下代码表达式将第一个索引固定，第二个索引移动，给出第一个（x, y）坐标：

In [31]: r_0 = r_i[0, :]

在最后一个轴上对索引进行切片是可选的；使用r_i[0]与r_i[0, :]具有相同的效果。

NumPy 允许使用另一个由整数或布尔值组成的 NumPy 数组来索引数组——这被称为花式索引。

如果您使用整数数组进行索引，NumPy 将解释这些整数为索引，并返回一个包含它们对应值的数组。如果我们用[0, 2, 3]索引包含 10 个元素的数组，我们将得到一个大小为 3 的数组，包含位置 0、2 和 3 的元素。以下代码展示了这一概念：

In [32]: a = np.array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])
In [33]: idx = np.array([0, 2, 3])
In [34]: a[idx]
Out[34]: array([9, 7, 6])

您可以通过为每个维度传递一个数组来在多个维度上使用花式索引。如果我们想提取元素[0, 2]和[1, 3]，我们必须将作用于第一个轴的所有索引打包在一个数组中，将作用于第二个轴的索引打包在另一个数组中。这可以在以下代码中看到：

In [35]: a = np.array([[0, 1, 2], [3, 4, 5],
                       [6, 7, 8], [9, 10, 11]])
In [36]: idx1 = np.array([0, 1])
In [37]: idx2 = np.array([2, 3])
In [38]: a[idx1, idx2]

提示

您还可以使用普通列表作为索引数组，但不能使用元组。例如，以下两个语句是等价的：

>>> a[np.array([0, 1])] # is equivalent to
>>> a[[0, 1]]

然而，如果您使用一个元组，NumPy 将解释以下语句为在多个维度上的索引：

>>> a[(0, 1)] # is equivalent to>>> a[0, 1]

索引数组不必是一维的；我们可以以任何形状从原始数组中提取元素。例如，我们可以从原始数组中选择元素以形成一个（2, 2）的数组，如下所示：

In [39]: idx1 = [[0, 1], [3, 2]]
In [40]: idx2 = [[0, 2], [1, 1]]
In [41]: a[idx1, idx2]
Out[41]: array([[ 0,  5],[10,  7]])

可以组合数组切片和花式索引功能。例如，如果我们想在坐标数组中交换 x 和 y 列，这很有用。在下面的代码中，第一个索引将遍历所有元素（一个切片），然后对于这些元素，我们首先提取位置 1（y）的元素，然后是位置 0（x）的元素：

In [42]: r_i = np.random(10, 2)
In [43]: r_i[:, [0, 1]] = r_i[:, [1, 0]]

当索引数组是布尔值时，有一些不同的规则。布尔数组将像掩码一样作用；每个对应于True的元素将被提取并放入输出数组中。这个过程如下所示：

In [44]: a = np.array([0, 1, 2, 3, 4, 5])
In [45]: mask = np.array([True, False, True, False, False, False])
In [46]: a[mask]
Out[46]: array([0, 2])

在处理多个维度时，同样适用这些规则。此外，如果索引数组与原始数组具有相同的形状，则对应于True的元素将被选中并放入结果数组中。

NumPy 中的索引操作是一个相对快速的操作。无论如何，当速度至关重要时，您可以使用稍微快一点的numpy.take和numpy.compress函数来挤出更多速度。numpy.take的第一个参数是我们想要操作的数组，第二个参数是我们想要提取的索引列表。最后一个参数是axis；如果没有提供，索引将作用于展平后的数组，否则将沿着指定的轴进行。以下代码显示了np.take的使用及其与花式索引的执行时间比较：

In [47]: r_i = np.random(100, 2)
In [48]: idx = np.arange(50) # integers 0 to 50
In [49]: %timeit np.take(r_i, idx, axis=0)
1000000 loops, best of 3: 962 ns per loop
In [50]: %timeit r_i[idx]
100000 loops, best of 3: 3.09 us per loop

使用布尔数组索引的类似但更快的方法是 numpy.compress，它的工作方式与 numpy.take 相同。以下是如何使用 numpy.compress 的示例：

In [51]: idx = np.ones(100, dtype='bool') # all True values
In [52]: %timeit np.compress(idx, r_i, axis=0)
1000000 loops, best of 3: 1.65 us per loop
In [53]: %timeit r_i[idx]
100000 loops, best of 3: 5.47 us per loop

广播

NumPy 的真正力量在于其快速的数学运算。NumPy 使用的策略是通过在匹配的数组之间进行逐元素计算来避免进入 Python。

无论何时你在两个数组（如乘积）上进行算术运算，如果两个操作数具有相同的形状，该操作将以逐元素的方式应用。例如，在乘以两个 (2, 2) 的数组时，操作将在对应元素对之间进行，产生另一个 (2, 2) 的数组，如下面的代码所示：

In [54]: A = np.array([[1, 2], [3, 4]])
In [55]: B = np.array([[5, 6], [7, 8]])
In [56]: A * B
Out[56]: array([[ 5, 12],[21, 32]])

如果操作数的形状不匹配，NumPy 将尝试使用某些规则来匹配它们——这被称为 广播 功能。如果一个操作数是一个单独的值，它将被应用到数组的每个元素上，如下面的代码所示：

In [57]: A * 2
Out[58]: array([[2, 4],
                [6, 8]])

如果操作数是另一个数组，NumPy 将从最后一个轴开始尝试匹配形状。例如，如果我们想将一个形状为 (3, 2) 的数组与一个形状为 (2,) 的数组组合，第二个数组将被重复三次以生成一个 (3, 2) 的数组。数组被 广播 以匹配另一个操作数的形状，如下面的图所示：

如果形状不匹配，例如通过将一个形状为 (3, 2) 的数组与一个形状为 (2, 2) 的数组组合，NumPy 将抛出异常。

如果某个轴的大小是 1，数组将在这个轴上重复，直到形状匹配。为了说明这一点，如果我们有一个形状如下的数组：

5, 10, 2

并且我们想将它与一个形状为 (5, 1, 2) 的数组广播，数组将在第二个轴上重复 10 次，如下所示：

5, 10, 2
5,  1, 2 → repeated
- - - -
5, 10, 2

我们之前看到，我们可以自由地重塑数组以添加大小为 1 的轴。在索引时使用 numpy.newaxis 常量将引入一个额外的维度。例如，如果我们有一个 (5, 2) 的数组，并且我们想将它与一个形状为 (5, 10, 2) 的数组结合，我们可以在中间添加一个额外的轴，如下面的代码所示，以获得兼容的 (5, 1, 2) 数组：

In [59]: A = np.random.rand(5, 10, 2)
In [60]: B = np.random.rand(5, 2)
In [61]: A * B[:, np.newaxis, :]

这个功能可以用来操作两个数组的所有可能的组合。其中一个应用是 外积。如果我们有以下两个数组：

a = [a1, a2, a3]
b = [b1, b2, b3]

外积是一个矩阵，包含两个数组元素所有可能的组合（i, j）的乘积，如下面的代码所示：

a x b = a1*b1, a1*b2, a1*b3
        a2*b1, a2*b2, a2*b3
        a3*b1, a3*b2, a3*b3

要使用 NumPy 计算这个，我们将在一个维度上重复元素 [a1, a2, a3]，在另一个维度上重复元素 [b1, b2, b3]，然后取它们的乘积，如下面的图所示：

我们的战略是将数组 a 从形状 (3,) 转换为形状 (3, 1)，将数组 b 从形状 (3,) 转换为形状 (1, 3)。两个数组在两个维度上广播，并使用以下代码相乘：

AB = a[:, np.newaxis] * b[np.newaxis, :]

这个操作非常快且非常有效，因为它避免了 Python 循环并能处理大量元素。

数学运算

NumPy 默认包含了广播中最常见的数学运算，从简单的代数到三角学、舍入和逻辑。例如，为了取数组中每个元素的平方根，我们可以使用numpy.sqrt函数，如下所示：

In [59]: np.sqrt(np.array([4, 9, 16]))
Out[59]: array([2., 3., 4.])

比较运算符在尝试根据条件过滤某些元素时非常有用。想象一下，我们有一个范围在[0, 1]内的随机数数组，我们想要提取大于 0.5 的数字。我们可以在数组上使用>运算符；结果将是一个布尔数组，如下所示：

In [60]: a = np.random.rand(5, 3)
In [61]: a > 0.5
Out[61]: array([[ True, False,  True],[ True,  True,  True],[False,  True,  True],[ True,  True, False],[ True,  True, False]], dtype=bool)

结果布尔数组然后可以被用作索引来检索大于 0.5 的元素，如下所示：

In [62]: a[a > 0.5]
In [63]: print(a[a>0.5])
[ 0.9755  0.5977  0.8287  0.6214  0.5669  0.9553  0.58940.7196  0.9200  0.5781  0.8281 ]

NumPy 还实现了ndarray.sum等方法，该方法在轴上对所有元素求和。如果我们有一个形状为(5, 3)的数组，我们可以使用ndarray.sum方法，如下所示，来在第一个轴、第二个轴或整个数组上的元素进行求和：

In [64]: a = np.random.rand(5, 3)
In [65]: a.sum(axis=0)
Out[65]: array([ 2.7454,  2.5517,  2.0303])
In [66]: a.sum(axis=1)
Out[66]: array([ 1.7498,  1.2491,  1.8151,  1.9320,  0.5814])
In [67]: a.sum() # With no argument operates on flattened array
Out[67]: 7.3275

注意，通过在一个轴上求和元素，我们消除了该轴。从上一个例子中，轴 0 上的求和产生一个(3,)数组，而轴 1 上的求和产生一个(5,)数组。

计算范数

我们可以通过计算一组坐标的范数来回顾本节中阐述的基本概念。对于二维向量，范数被定义为：

norm = sqrt(x² + y²)

给定一个包含 10 个坐标(x, y)的数组，我们想要找到每个坐标的范数。我们可以通过以下步骤来计算范数：

1. 对坐标进行平方：得到一个包含(x**2, y**2)元素的数组。

2. 使用numpy.sum在最后一个轴上对这些值进行求和。

3. 使用numpy.sqrt逐元素取平方根。

最终表达式可以压缩成一行：

In [68]: r_i = np.random.rand(10, 2)
In [69]: norm = np.sqrt((r_i ** 2).sum(axis=1))
In [70]: print(norm)
[ 0.7314  0.9050  0.5063  0.2553  0.0778   0.91431.3245  0.9486  1.010   1.0212]

用 NumPy 重写粒子模拟器

在本节中，我们将通过用 NumPy 重写其部分来优化我们的粒子模拟器。从我们在第一章中进行的性能分析来看，我们程序中最慢的部分是包含在ParticleSimulator.evolve方法中的以下循环：

for i in range(nsteps):
  for p in self.particles:

    norm = (p.x**2 + p.y**2)**0.5
    v_x = (-p.y)/norm
    v_y = p.x/norm

    d_x = timestep * p.ang_speed * v_x
    d_y = timestep * p.ang_speed * v_y

    p.x += d_x
    p.y += d_y

我们可能会注意到循环的主体仅作用于当前粒子。如果我们有一个包含粒子位置和角速度的数组，我们可以使用广播操作重写循环。相比之下，时间步长上的循环依赖于前一步，不能以并行方式处理。

因此，将所有数组坐标存储在一个形状为(nparticles, 2)的数组中，并将角速度存储在一个形状为(nparticles,)的数组中是很自然的。我们将这些数组称为r_i和ang_speed_i，并使用以下代码进行初始化：

r_i = np.array([[p.x, p.y] for p in self.particles])
ang_speed_i = np.array([p.ang_speed for p in self.particles])

速度方向，垂直于向量（x, y），被定义为：

v_x = -y / norm
v_y = x / norm

可以使用在“NumPy 入门”标题下的“计算范数”部分中说明的策略来计算范数。最终的表达式在以下代码行中显示：

norm_i = ((r_i ** 2).sum(axis=1))**0.5

对于 (-y, x) 分量，我们首先需要在 r_i 中交换 x 和 y 列，然后将第一列乘以 -1，如下所示：

v_i = r_i[:, [1, 0]] / norm_i
v_i[:, 0] *= -1

要计算位移，我们需要计算 v_i、ang_speed_i 和 timestep 的乘积。由于 ang_speed_i 的形状为 (nparticles,)，它需要一个新轴才能与形状为 (nparticles, 2) 的 v_i 操作。我们将使用 numpy.newaxis 常量来完成此操作，如下所示：

d_i = timestep * ang_speed_i[:, np.newaxis] * v_i
r_i += d_i

循环外部，我们需要按照以下方式更新粒子实例的新坐标 x 和 y：

for i, p in enumerate(self.particles):
  p.x, p.y = r_i[i]

总结来说，我们将实现一个名为 ParticleSimulator.evolve_numpy 的方法，并将其与重命名为 ParticleSimulator.evolve_python 的纯 Python 版本进行基准测试。完整的 ParticleSimulator.evolve_numpy 方法如下所示：

def evolve_numpy(self, dt):
  timestep = 0.00001
  nsteps = int(dt/timestep)

  r_i = np.array([[p.x, p.y] for p in self.particles])
  ang_speed_i = np.array([p.ang_speed for p in self.particles])

  for i in range(nsteps):

    norm_i = np.sqrt((r_i ** 2).sum(axis=1))
    v_i = r_i[:, [1, 0]]
    v_i[:, 0] *= -1
    v_i /= norm_i[:, np.newaxis]
    d_i = timestep * ang_speed_i[:, np.newaxis] * v_i
    r_i += d_i

    for i, p in enumerate(self.particles):
      p.x, p.y = r_i[i]

我们还更新了基准测试，以便方便地更改粒子数量和模拟方法，如下所示：

def benchmark(npart=100, method='python'):
  particles = [Particle(uniform(-1.0, 1.0),uniform(-1.0, 1.0),uniform(-1.0, 1.0))for i in range(npart)]

  simulator = ParticleSimulator(particles)

  if method=='python':
    simulator.evolve_python(0.1)

  elif method == 'numpy':
    simulator.evolve_numpy(0.1)

我们可以在 IPython 会话中运行更新后的基准测试，如下所示：

In [1]: from simul import benchmark
In [2]: %timeit benchmark(100, 'python')
1 loops, best of 3: 614 ms per loop
In [3]: %timeit benchmark(100, 'numpy')
1 loops, best of 3: 415 ms per loop

我们有一些改进，但看起来速度提升并不大。NumPy 的强大之处在于处理大型数组。如果我们增加粒子数量，我们会注意到更明显的性能提升。我们可以使用以下代码重新运行具有更多粒子数量的基准测试：

In [4]: %timeit benchmark(1000, 'python')
1 loops, best of 3: 6.13 s per loop
In [5]: %timeit benchmark(1000, 'numpy')
1 loops, best of 3: 852 ms per loop

下一个图中的图表是通过运行具有不同粒子数量的基准测试生成的：

图表显示，两种实现都与粒子大小成线性关系，但纯 Python 版本（用菱形表示）的运行时间比 NumPy 版本（用圆圈表示）增长得更快；在更大的尺寸上，NumPy 的优势更加明显。一般来说，当使用 NumPy 时，应尝试将事物打包到大型数组中，并使用广播功能按组进行计算。

使用 numexpr 达到最佳性能

当处理复杂表达式时，NumPy 会将中间结果存储在内存中。David M. Cooke 编写了一个名为 numexpr 的包，该包可以即时优化和编译数组表达式。它通过优化 CPU 缓存的使用并利用多处理器来实现。

其用法通常很简单，基于一个单一的功能——numexpr.evaluate。该函数将包含数组表达式的字符串作为其第一个参数。语法基本上与 NumPy 相同。例如，我们可以以下这种方式计算一个简单的 a + b * c 表达式：

a = np.random.rand(10000)
b = np.random.rand(10000)
c = np.random.rand(10000)
d = ne.evaluate('a + b * c')

numexpr 包在几乎所有情况下都会提高性能，但要实现实质性的优势，你应该使用它来处理大型数组。一个涉及大型数组的应用是计算 距离矩阵。在粒子系统中，距离矩阵包含粒子之间所有可能距离。为了计算它，我们首先应该计算连接任意两个粒子（i，j）的所有向量，如下定义：

x_ij = x_j - x_i
y_ij = y_i - y_j

然后，我们通过取其范数来计算这个向量的长度，如下面的代码所示：

d_ij = sqrt(x_ij**2 + y_ij**2)

我们可以通过使用常规广播规则（操作类似于外积）来用 NumPy 实现这一点：

r = np.random.rand(10000, 2)
r_i = r[:, np.newaxis]
r_j = r[np.newaxis, :]
r_ij = r_j – r_i

最后，我们使用以下代码行在最后一个轴上计算范数：

d_ij = np.sqrt((r_ij ** 2).sum(axis=2))

使用 numexpr 语法重写相同的表达式非常简单。numexpr 包不支持数组表达式的切片，因此我们首先需要通过添加一个额外的维度来准备广播操作符，如下所示：

r = np.random(10000, 2)
r_i = r[:, np.newaxis]
r_j = r[np.newaxis, :]

在这一点上，我们应该尝试在一个单独的表达式中尽可能多地组合操作，以允许进行显著的优化。

大多数 NumPy 数学函数也存在于 numexpr 中；然而，有一个限制。减少操作——比如求和——必须放在最后执行。因此，我们必须首先计算总和，然后退出 numexpr，并在另一个表达式中计算平方根。这些操作的 numexpr 代码如下：

d_ij = ne.evaluate('sum((r_j – r_i)**2, 2)')
d_ij = ne.evaluate('sqrt(d_ij)')

numexpr 编译器将通过避免存储中间结果并利用多个处理器来优化内存使用。在 distance_matrix.py 文件中，你可以找到实现距离矩阵计算两种版本的函数：distance_matrix_numpy 和 distance_matrix_numexpr。我们可以这样导入和基准测试它们：

from distance_matrix import (distance_matrix_numpy,
                             distance_matrix_numexpr)
%timeit distance_matrix_numpy(10000)
1 loops, best of 3: 3.56 s per loop
%timeit distance_matrix_numexpr(10000)
1 loops, best of 3: 858 ms per loop

通过简单地使用 numexpr 复制表达式，我们能够在实际场景中获得比标准 NumPy 高出 4.5 倍的性能提升。numexpr 包可以在你需要优化涉及大型数组和复杂操作的 NumPy 表达式时使用，并且你可以通过最小的代码更改来实现这一点。

摘要

在本章中，我们学习了如何操作 NumPy 数组以及如何使用数组广播编写快速数学表达式。这些知识将帮助你设计更好的程序，同时获得巨大的性能提升。我们还介绍了 numexpr 库，以最小的努力进一步加快我们的计算速度。

NumPy 在处理独立输入集时工作得非常好，但当表达式变得复杂且无法分割为逐元素操作时，它就不适用了。在这种情况下，我们可以通过使用 Cython 包将 Python 与 C 进行接口，利用 Python 作为粘合语言的能力。

第三章。使用 Cython 的 C 性能

Cython 是一种通过向函数、变量和类添加静态类型来扩展 Python 的语言。Cython 结合了 Python 的简洁性和 C 语言的效率。在将脚本重写为 Cython 之后，您可以编译它们为 C 或 C++，以简单直接的方式生成高效代码。

Cython 还充当 Python 和 C 之间的桥梁，因为它可以用来创建对外部 C 代码的接口。通过创建绑定，您可以在脚本中重用快速的 C 例程，有效地使用 Python 作为粘合语言。

在本章中，我们将学习：

• Cython 语法基础

• 如何编译 Cython 程序

• 如何使用 静态类型 生成快速代码

• 如何通过使用类型化的 memoryviews 高效地操作数组。

最后，我们将应用我们新的 Cython 技能来分析和优化粒子模拟器。

虽然对 C 语言有最小了解有帮助，但本章仅关注在 Python 优化背景下使用 Cython。因此，它不需要任何 C 语言背景知识。

编译 Cython 扩展

按照设计，Cython 语法是 Python 的超集。Cython 通常可以编译 Python 模块而不需要任何更改。Cython 源文件具有 .pyx 扩展名，并且可以使用 cython 命令编译为 C 语言。

我们的第一篇 Cython 脚本将包含一个简单的函数，该函数将 Hello, World! 打印到输出。创建一个名为 hello.pyx 的新文件，包含以下代码：

def hello():
  print('Hello, World!')

cython 命令将读取 hello.pyx 并生成 hello.c 文件：

$ cython hello.pyx

要将 hello.c 编译为 Python 扩展模块，我们将使用 gcc 编译器。我们需要添加一些依赖于操作系统的特定于 Python 的编译选项。在 Ubuntu 13.10 上，使用默认的 Python 安装，您可以使用以下选项进行编译：

$ gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing -lm -I/usr/include/python3.3/ -o hello.so hello.c

这将生成一个名为 hello.so 的文件：一个可以从 Python 导入的 C 扩展模块。

>>> import hello
>>> hello.hello()
Hello, World!

注意

Cython 既可以接受 Python 2，也可以接受 Python 3 作为输入和输出语言。换句话说，您可以使用 -3 选项编译 Python 3 的 hello.pyx 文件：

$ cython -3 hello.pyx

生成的 hello.c 可以通过在 gcc 中包含相应的头文件（使用 -I 选项）而不做任何更改地编译为 Python 2 和 Python 3：

$ gcc -I/usr/include/python3.3 # ... other options
$ gcc -I/usr/include/python2.7 # ... other options

通过使用 distutils（标准的 Python 打包工具），Cython 程序可以以更直接的方式编译。通过编写 setup.py 脚本，我们可以直接将 .pyx 文件编译为扩展模块。为了编译我们的 hello.pyx 示例，我们需要编写一个包含以下代码的 setup.py：

from distutils.core import setup
from Cython.Build import cythonize

setup(
  name='Hello',ext_modules = cythonize('hello.pyx'),)

在前面代码的前两行中，我们导入了 setup 函数和 cythonize 辅助函数。setup 函数包含一些键值对，告诉 distutils 应用程序的名称以及需要构建哪些扩展。

cythonize 辅助函数接受一个字符串或包含我们想要编译的 Cython 模块的字符串列表。您还可以使用以下代码使用 glob 模式：

cythonize(['hello.pyx', 'world.pyx', '*.pyx'])

要使用 distutils 编译我们的扩展模块，您可以执行以下代码的 setup.py 脚本：

$ python setup.py build_ext --inplace

build_ext 选项指示脚本构建 ext_modules 中指定的扩展模块，而 --inplace 选项将输出文件 hello.so 放置在源文件相同的目录中（而不是构建目录）。

Cython 模块可以自动使用 pyximport 编译。通过在脚本开头添加 pyximport.install()（或在解释器中发出该命令），您可以直接导入 .pyx 文件；pyximport 将透明地编译相应的 Cython 模块。

>>> import pyximport
>>> pyximport.install()
>>> import hello # This will compile hello.pyx

不幸的是，pyximport 并不适用于所有类型的配置（例如涉及 C 和 Cython 文件组合的情况），但它对于测试简单的脚本来说很有用。

自 0.13 版本以来，IPython 包含了 cythonmagic 扩展，可以交互式地编写和测试一系列 Cython 语句。您可以使用 load_ext 在 IPython 壳中加载扩展：

In [1]: %load_ext cythonmagic

一旦加载了扩展，您就可以使用 %%cython 单元格魔法 来编写多行 Cython 片段。在以下示例中，我们定义了一个 hello_snippet 函数，该函数将被编译并添加到 session 命名空间中：

In [2]: %%cython
  ....:  def hello_snippet():
   ....:    print("Hello, Cython!")
  ....:
In [3]:  hello_snippet()
Hello,  Cython!

添加静态类型

在 Python 中，变量有一个关联的类型，该类型在执行过程中可以改变。虽然这个特性是可取的，因为它使语言更加灵活，但解释器需要进行类型检查和方法查找来正确处理变量之间的操作——这是一个引入了显著开销的额外步骤。Cython 通过静态类型声明扩展了 Python 语言；这样，它可以通过避免 Python 解释器来生成高效的 C 代码。

在 Cython 中声明数据类型的主要方式是通过使用 cdef 语句。cdef 关键字可以在多个上下文中使用：声明变量、函数和扩展类型（cdef 类）。

变量

在 Cython 中，您可以通过在变量前加上 cdef 和相应的类型来声明变量的类型。例如，我们可以以下这种方式声明变量 i 为 16 位整数：

cdef int i

cdef 语句支持在同一行上使用多个变量名，以及可选的初始化值，如下面的代码行所示：

cdef double a, b = 2.0, c = 3.0

与标准变量相比，类型化变量被处理得不同。在 Python 中，变量通常被视为 标签，指代内存中的对象。在任何程序点，我们都可以将一个字符串赋给变量，如下所示：

a = 'hello'

字符串 hello 将绑定到变量 a。在程序的不同位置，我们可以将另一个值赋给相同的变量，例如一个整数：

a = 1

Python 会将整数对象 1 赋给变量 a 而不会出现任何问题。

带有类型的变量可以被认为是更类似于 数据容器；我们在变量中 存储 值，并且只有相同类型的值被允许进入。例如，如果我们将变量 a 声明为 int 类型的变量，然后我们尝试将其赋值为 double，Cython 将触发错误，如下面的代码所示：

In [4]: %%cython
  ....: cdef int i
  ....: i = 3.0
  ....:
# Output has been cut
...cf4b.pyx:2:4 Cannot assign type 'double' to 'int'

静态类型允许有用的优化。如果我们声明在循环中使用的索引为整数，Cython 将将循环重写为纯 C，而无需进入 Python 解释器。在下面的示例中，我们迭代 100 次，每次都递增 int 变量 j：

In [5]: %%cython
  ....: def example():
  ....:   cdef int i, j=0
  ....:   for i in range(100):....:      j += 1
  ....:   return j
  ....:
In [6]: example()
Out[6]: 100

为了了解改进有多大，我们将将其速度与类似的纯 Python 循环进行比较：

In [7]: def example_python():
  ....:    j=0
  ....:    for i in range(100):....:      j += 1
  ....:    return j
  ....:
In [8]: %timeit example()
10000000 loops, best of 3: 25 ns per loop
In [9]: %timeit example_python()
100000 loops, best of 3: 2.74 us per loop

通过编写带有类型信息的循环获得的加速效果高达 100 倍！这是因为 Cython 循环首先被转换为纯 C，然后转换为高效的机器代码，而 Python 循环仍然依赖于慢速的解释器。

我们可以声明任何可用的 C 类型的变量，我们也可以通过使用 C 结构体、枚举和 typedef 来定义自定义类型。一个有趣的例子是，如果我们声明一个变量为 object 类型，该变量将接受任何类型的 Python 对象：

cdef object a_py
# both 'hello' and 1 are Python objects
a_py = 'hello'
a_py = 1

有时，某些类型的变量是兼容的（例如 float 和 int 数字），但并不完全相同。在 Cython 中，可以通过将目标类型用 < 和 > 括号包围来在类型之间进行转换（cast），如下面的代码片段所示：

cdef int a = 0
cdef double b
b = <double> a

函数

你可以通过在参数名称前指定类型来向 Python 函数的参数添加类型信息。这样的函数将像常规 Python 函数一样工作并执行，但其参数将进行类型检查。我们可以编写一个 max_python 函数，以下是这样返回两个整数之间较大值的示例：

def max_python(int a, int b):
  return a if a > b else b

那个函数除了类型检查外没有提供太多好处。为了利用 Cython 优化，我们必须使用 cdef 语句和可选的返回类型来声明该函数，如下面的代码所示：

cdef int max_cython(int a, int b):
  return a if a > b else b

以这种方式声明的函数会被转换为原生 C 函数，这些函数不能从 Python 中调用。与 Python 函数相比，它们的开销要小得多，使用它们会导致性能显著提升。它们的范围限制在相同的 Cython 文件中，除非它们在定义文件中公开（参考 共享声明 部分）。

Cython 允许你定义既可以从 Python 调用也可以转换为原生 C 函数的函数。如果你使用关键字 cpdef 声明一个函数，Cython 将生成该函数的两个版本——一个可供解释器使用的 Python 版本和一个可从 Cython 使用的快速 C 函数，从而实现便利性和速度。cpdef 语法与 cdef 相当，如下所示：

cpdef int max_hybrid(int a, int b):
  return a if a > b else b

有时，即使有 C 函数，调用开销也可能成为性能问题，尤其是在关键循环中多次调用同一函数时。当函数体较小时，在函数定义前添加 inline 关键字是很方便的；函数调用将被移除并替换为函数体。例如，我们下面的 max 函数是 内联 的一个很好的候选：

cdef inline int max_inline(int a, int b):
  return a if a > b else b

类

cdef 关键字也可以放在类定义的前面来创建一个 扩展类型。扩展类型类似于 Python 类，但其属性必须具有类型，并且存储在高效的 C 结构体 中。

我们可以使用 cdef class 语句定义扩展类型，并在类体中声明其属性。例如，我们可以创建一个扩展类型 Point，如下面的代码所示，它存储两个坐标（x，y）的类型为 double：

cdef class Point:
 cdef double x
 cdef double y

  def __init__(self, double x,double y):
    self.x = x
    self.y = y

在类方法中访问声明的属性允许 Cython 通过将其替换为直接访问 struct 字段来避免 Python 属性查找。这样，属性访问就变得非常快速。

为了利用 struct 访问，Cython 需要知道在编译时变量是一个扩展类型。你可以在任何可以使用标准类型（如 double、float、int）的地方使用扩展类型名称（如 Point）。例如，如果我们想编写一个计算 Point 的 norm 的 Cython 函数，我们必须将输入变量声明为 Point，如下面的代码所示：

cdef double norm(Point p):
  return p.x**2 + p.y**2

默认情况下，对属性的访问限制在 Cython 代码中。如果你尝试从 Python 访问扩展类型属性，你会得到一个 AttributeError，如下所示：

>>> a = Point(0.0, 0.0)
>>> a.x
AttributeError: 'Point' object has no attribute 'x'

为了从 Python 代码中访问属性，你必须使用属性声明中的 public（用于读写访问）或 readonly 说明符，如下面的代码所示：

cdef class Point:
  cdef public double x

扩展类型不支持添加额外的属性。对于这个问题的一个解决方案是子类化扩展类型，创建一个派生的 Python 类。

分享声明

当编写你的 Cython 模块时，你可能希望在一个单独的文件中封装通用的函数和类型。Cython 允许你通过编写 定义文件 来使用 cimport 语句重用这些组件。

假设我们有一个包含 max 和 min 函数的模块，并且我们想在多个 Cython 程序中重用这些函数。如果我们简单地编写一个 .pyx 文件——也称为 实现文件——声明的函数将局限于同一模块中。

注意

定义文件也用于将 Cython 与外部 C 代码接口。其思路是将定义文件中的类型和函数原型复制过来，并将实现留给外部 C 代码。

要共享这些函数，我们需要编写一个定义文件，具有.pxd扩展名。这样的文件只包含我们想要与其他模块共享的类型和函数原型——一个公共接口。我们可以在名为mathlib.pxd的文件中编写我们的max和min函数的原型，如下所示：

cdef int max(int a, int b)
cdef int min(int a, int b)

如您所见，我们只编写了函数名和参数，而没有实现函数体。

函数实现将放入具有相同基本名称但.pyx扩展名的实现文件中——mathlib.pyx：

cdef int max(int a, int b):
  return a if a > b else b

cdef int min(int a, int b):
  return a if a < b else b

mathlib模块现在可以从另一个 Cython 模块导入。

要测试我们的 Cython 模块，我们将创建一个名为distance.pyx的文件，其中包含一个名为chebyshev的函数。该函数将计算两点之间的 Chebyshev 距离，如下面的代码所示。两点坐标（x1, y1）和（x2, y2）之间的 Chebyshev 距离定义为每个坐标之间差异的最大值。

max(abs(x1 – x2), abs(y1 – y2))

要实现chebyshev函数，我们将使用max函数，它在mathlib.pxd中声明，通过使用cimport语句导入，如下面的代码片段所示：

from mathlib cimport max

def chebyshev(int x1,int y1,int x2,int y2):
  return max(abs(x1 - x2), abs(y1 - y2))

cimport语句将读取hello.pxd，并将max定义用于生成distance.c文件。

使用数组

数值和高性能计算通常使用数组。Cython 提供了一种简单的方法来与之交互，从 C 数组的底层方法到更通用的类型化内存视图。

C 数组和指针

C 数组是一系列大小相同的项，在内存中连续存储。在深入了解细节之前，了解（或复习）C 语言中内存的管理方式是有帮助的。

C 语言中的变量就像容器。当创建一个变量时，会在内存中预留一个空间来存储其值。例如，如果我们创建一个包含 64 位浮点数（double）的变量，程序将分配 64 位（16 字节）的内存。这部分内存可以通过该内存位置的地址来访问。

要获取变量的地址，我们可以使用地址运算符，用&符号表示。我们还可以使用printf函数，如下所示，在libc.stdio Cython 模块中可用，以打印该变量的地址：

In [1]: %%cython
  ...: cdef double a
  ...: from libc.stdio cimport printf
  ...: printf("%p", &a)
  ...:
0x7fc8bb611210

内存地址可以存储在特殊的变量中——指针，通过在变量名前放置*前缀来声明，如下所示：

from libc.stdio cimport printf
cdef double a
cdef double *a_pointer
a_pointer = &a # They are of the same data type

如果我们有一个指针，并且想要获取它所指向的地址中的值，我们可以使用解引用运算符，用*符号表示，如下面的代码所示。请注意，在此上下文中使用的*与在变量声明中使用的*有不同的含义。

cdef double a
cdef double *a_pointer
a_pointer = &a
a = 3.0
print(*a_pointer) # prints 3.0

当声明一个 C 数组时，程序会分配足够的空间来包含指定大小的多个元素。例如，要创建一个包含 10 个double值（每个 8 字节）的数组，程序将在内存中保留8 * 10 = 80字节的连续空间。在 Cython 中，我们可以使用以下语法声明这样的数组：

cdef double arr[10]

我们也可以使用以下语法声明一个多维数组，例如一个有 5 行 2 列的数组：

cdef double arr[5][2]

内存将在一个单独的内存块中分配，一行接一行。这种顺序通常被称为行主序，如下图中所示。数组也可以按列主序排序，正如在 FORTRAN 编程语言中发生的那样。

小贴士

数组排序有重要的影响。当我们对 C 数组的最后一个维度进行迭代时，我们访问连续的内存块（在我们的例子中是 0, 1, 2, 3 …），而当我们对第一个维度进行迭代时，我们会跳过一些位置（0, 2, 4, 6, 8, 1 …）。你应该始终尝试连续访问内存，因为这优化了缓存的使用。

我们可以通过使用标准索引来存储和检索数组中的元素，C 数组不支持花哨的索引或切片：

arr[0] = 1.0

C 数组也可以用作指针。实际上，arr变量是一个指向数组第一个元素的指针。我们可以验证数组第一个元素的地址与变量arr中包含的地址相同：

In [1]: %%cython
  ...: from libc.stdio cimport printf
  ...: cdef double arr[10]
  ...: printf("%p\n", arr)
  ...: printf("%p\n", &arr[0])
  ...:
0x7ff6de204220
0x7ff6de204220

当与现有的 C 库接口或需要精细控制内存时，应使用 C 数组和指针。对于更常见的用例，可以使用 NumPy 数组或类型化内存视图。

NumPy 数组

NumPy 数组可以作为常规 Python 对象在 Cython 中使用，通过使用它们的已优化的广播操作。

当我们想要高效地遍历数组时，问题就出现了。当我们对 NumPy 数组进行索引操作时，解释器级别会发生一些其他操作，导致大量开销。Cython 可以通过直接作用于 NumPy 数组使用的底层内存区域来优化这些索引操作，使我们能够像 C 数组一样处理它们。

NumPy 数组支持以ndarray数据类型的形式提供。我们首先必须cimport numpy模块。我们将其分配给名称c_np以区分它和常规的numpy Python 模块，如下所示：

cimport numpy as c_np

我们现在可以通过指定数组元素的类型和维度数，使用一种称为缓冲区语法的特殊语法来声明 NumPy 数组。要声明一个类型为double的两维数组，我们可以使用以下代码：

cdef c_np.ndarray[double, ndim=2] arr

以这种方式定义的数组将通过直接作用于底层内存区域来进行索引；这种操作将避免 Python 解释器给我们带来巨大的速度提升。

在下一个示例中，我们将展示缓冲区语法的用法，并将其与常规 Python 版本进行比较。

我们首先编写了 numpy_bench_py 函数，该函数将 py_arr 的每个元素增加 1000。我们将索引 i 声明为整数，以避免 for 循环的开销：

In [1]:  %%cython
  ...:  import numpy as np
  ...:  def numpy_bench_py():...:    py_arr = np.random.rand(1000)
  ...:    cdef int i
  ...:    for i in range(1000):
  ...:       py_arr[i] += 1

然后，我们使用缓冲区语法编写相同的函数。注意，在定义 c_arr 变量使用 c_np.ndarray 之后，我们可以从 numpy Python 模块给它分配一个数组：

In [2]:  %%cython
  ...:  import numpy as np
  ...:  cimport numpy as c_np
  ...:  def numpy_bench_c():
  ...:    cdef c_np.ndarray[double, ndim=1] c_arr
  ...:    c_arr = np.random.rand(1000)
  ...:    cdef int i
  ...:
  ...:    for i in range(1000):
  ...:       c_arr[i] += 1

我们可以使用 timeit 来计时结果，获得令人印象深刻的 50 倍速度提升：

In [10]: %timeit numpy_bench_c()
100000 loops, best of 3: 11.5 us per loop
In [11]: %timeit numpy_bench_py()
1000 loops, best of
 3: 603 us per loop

类型化内存视图

C 和 NumPy 数组都是作用于内存区域的对象。Cython 提供了一个通用的对象——类型化内存视图——以访问数组和其他暴露所谓 缓冲区接口 的数据结构，例如内置的 bytes、bytearray 和 array.array。

内存视图是一个维护对某个内存区域的引用的对象。它实际上并不拥有内存，但它可以读取和更改其内容（它是一个 视图）。通过使用类型化内存视图，我们可以以相同的方式与 C 和 NumPy 数组进行交互。

内存视图可以使用特殊语法定义。我们可以以下这种方式定义一个 int 类型的内存视图和一个 2D double 类型的内存视图：

cdef int[:] a
cdef double[:, :] b

相同的语法适用于函数定义、类属性等。任何暴露缓冲区接口的对象都将自动绑定到内存视图。我们可以通过以下简单的赋值来绑定内存视图到数组：

import numpy as np

cdef int[:] arr
arr_np = np.zeros(10, dtype='int32')
arr = arr_np # We bind the array to the memoryview

新的内存视图将与 NumPy 数组共享数据。数组元素的变化将在两个数据结构之间共享：

arr[2] = 1 # Changing memoryview
print(arr_np)
# [0 0 1 0 0 0 0 0 0 0]

在某种意义上，内存视图是 NumPy 数组的一种推广。正如我们在第二章中看到的，“使用 NumPy 的快速数组操作”，切片 NumPy 数组不会复制数据，而是返回对同一内存区域的视图。

内存视图还支持以下标准 NumPy 语法进行数组切片：

cdef int[:, :, :] a
arr[0, :, :] # Is a 2-dimensional memoryview
arr[0, 0, :] # Is a 1-dimensional memoryview
arr[0, 0, 0] # Is an int

要在内存视图和另一个对象之间复制数据，你可以使用类似于切片赋值的语法，如下面的代码所示：

import numpy as np

cdef double[:, :] b
cdef double[:] r
b = np.random.rand(10, 3)
r = np.zeros(3, dtype='float64')

b[0, :] = r # Copy the value of r in the first row of b

在下一节中，我们将使用类型化内存视图来处理我们的粒子模拟器应用程序中的数组。

Cython 中的粒子模拟器

现在我们对 Cython 的工作原理有了基本的了解，我们可以重写 ParticleSimulator.evolve 方法。多亏了 Cython，我们可以将我们的循环转换为 C，从而消除由 Python 解释器引入的开销。

在第二章中，我们编写了一个相当高效的 evolve 方法版本，使用 NumPy。我们可以将旧版本重命名为 evolve_numpy 以区分新旧版本：

  def evolve_numpy(self, dt):
    timestep = 0.00001
    nsteps = int(dt/timestep)

    r_i = np.array([[p.x, p.y] for p in self.particles])    
    ang_speed_i = np.array([p.ang_speed for pin self.particles])
    v_i = np.empty_like(r_i)

    for i in range(nsteps):
      norm_i = np.sqrt((r_i ** 2).sum(axis=1))

      v_i = r_i[:, [1, 0]]
      v_i[:, 0] *= -1
      v_i /= norm_i[:, np.newaxis]        

      d_i = timestep * ang_speed_i[:, np.newaxis] * v_i

      r_i += d_i

    for i, p in enumerate(self.particles):
      p.x, p.y = r_i[i]

我们希望将此代码转换为 Cython。我们的策略将是利用快速索引操作，通过删除 NumPy 数组广播，从而回到基于索引的算法。由于 Cython 生成高效的 C 代码，我们可以自由地使用尽可能多的循环，而不会产生任何性能惩罚。

作为设计选择，我们可以决定将循环封装在一个函数中，我们将用名为 cevolve.pyx 的 Cython 模块重写这个函数。该模块将包含一个单一的 Python 函数 c_evolve，它将接受粒子位置、角速度、时间步长和步数作为输入。

起初，我们并没有添加类型信息；我们只是想隔离函数并确保我们可以无错误地编译我们的模块。

# file: simul.py
# ... other code
  def evolve_cython(self, dt):
    timestep = 0.00001
    nsteps = int(dt/timestep)

    r_i = np.array([[p.x, p.y] for p in self.particles])    
    ang_speed_i = np.array([p.ang_speed forp in self.particles])

    c_evolve(r_i, ang_speed_i, timestep, nsteps)

    for i, p in enumerate(self.particles):
      p.x, p.y = r_i[i]

# file: cevolve.pyx
import numpy as np

def c_evolve(r_i, ang_speed_i, timestep, nsteps):
  v_i = np.empty_like(r_i)

  for i in range(nsteps):
    norm_i = np.sqrt((r_i ** 2).sum(axis=1))

    v_i = r_i[:, [1, 0]]
    v_i[:, 0] *= -1
    v_i /= norm_i[:, np.newaxis]        

    d_i = timestep * ang_speed_i[:, np.newaxis] * v_i

    r_i += d_i

注意，我们不需要为 c_evolve 返回值，因为值是在 r_i 数组中就地更新的。我们可以通过稍微改变我们的基准函数，将无类型的 Cython 版本与旧的 NumPy 版本进行基准测试，如下所示：

def benchmark(npart=100, method='python'):
  particles = [Particle(uniform(-1.0, 1.0),uniform(-1.0, 1.0),uniform(-1.0, 1.0))for i in range(npart)]

  simulator = ParticleSimulator(particles)
  if method=='python':
    simulator.evolve_python(0.1)

 if method == 'cython':
 simulator.evolve_cython(0.1)

  elif method == 'numpy':
 simulator.evolve_numpy(0.1)

我们可以在 IPython 壳中计时不同的版本：

In [4]: %timeit benchmark(100, 'cython')
1 loops, best of 3: 401 ms per loop
In [5]: %timeit benchmark(100, 'numpy')
1 loops, best of 3: 413 ms per loop

这两个版本的速度相同。编译不带静态类型的 Cython 模块并不比纯 Python 有任何优势。下一步，是声明所有重要变量的类型，以便 Cython 可以执行其优化。

我们可以先从添加函数参数的类型开始。我们将声明数组为包含 double 值的已类型化内存视图。值得一提的是，如果我们传递一个 int 或 float32 类型的数组，转换不会自动发生，我们会得到一个错误。

def c_evolve(double[:, :] r_i, double[:] ang_speed_i,double timestep, int nsteps):

在这一点上，我们想要重写粒子和时间步的循环。我们可以将迭代变量 i、j 和粒子数 nparticles 声明为 int：

  cdef int i, j
  cdef int nparticles = r_i.shape[0]

在这一点上，算法与纯 Python 版本非常相似；我们遍历粒子和时间步，并计算每个粒子坐标的速度和位移向量，如下所示：

  for i in range(nsteps):
    for j in range(nparticles):
      x = r_i[j, 0]
      y = r_i[j, 1]
      ang_speed = ang_speed_i[j]

      norm = sqrt(x ** 2 + y ** 2)

      vx = (-y)/norm
      vy = x/norm

      dx = timestep * ang_speed * vx
      dy = timestep * ang_speed * vy

      r_i[j, 0] += dx
      r_i[j, 1] += dy

在之前的代码中，我们添加了 x、y、ang_speed、norm、vx、vy、dx 和 dy 变量。为了避免 Python 解释器的开销，我们必须在函数开始时声明它们对应的类型，如下所示：

cdef double norm, x, y, vx, vy, dx, dy, ang_speed

我们还使用了一个名为 sqrt 的函数来计算范数。如果我们使用 math 模块或 numpy 中的 sqrt，我们又会将一个慢速的 Python 函数包含在我们的关键循环中，从而降低我们的性能。标准 C 库中有一个快速的 sqrt，已经包含在 libc.math Cython 模块中：

from libc.math cimport sqrt

重新编译后，我们可以重新运行我们的基准测试来评估我们的改进，如下所示：

In [4]: %timeit benchmark(100, 'cython')
100 loops, best of 3: 13.4 ms per loop
In [5]: %timeit benchmark(100, 'numpy')
1 loops, best of 3: 429 ms per loop

对于小的粒子数，速度提升非常巨大，我们获得了比上一个版本 40 倍的性能提升。然而，我们也应该尝试使用更多的粒子来测试性能缩放，如下面的代码所示：

In [2]: %timeit benchmark(1000, 'cython')
10 loops, best of 3: 134 ms per loop
In [3]: %timeit benchmark(1000, 'numpy')
1 loops, best of 3: 877 ms per loop

随着粒子数的增加，两个版本的速度越来越接近。通过将粒子大小增加到 1000，我们已经将速度提升降低到更适度的 6 倍。这很可能是由于随着粒子数的增加，Python for 循环的开销与其它操作的速度相比变得越来越不显著。

分析 Cython

Cython 给我们一个很棒的工具，可以快速找到由于 Python 解释器引起的缓慢区域——一个称为 注释视图 的功能。我们可以通过使用 -a 选项编译 Cython 文件来打开此功能，使用以下命令行。Cython 将生成一个包含我们代码的 HTML 文件，并附有有用的注释：

$ cython -a cevolve.pyx
$ google-chrome cevolve.html

下面的截图显示的 HTML 文件显示了我们的 Cython 文件逐行内容：

每一行都有不同深浅的黄色背景色；深色表示代码有很多与解释器相关的调用，而白色行则被转换为纯 C 代码。由于解释器调用通常较慢，目标是使函数体尽可能为白色。通过点击任何一行，我们可以看到 Cython 编译器生成的 C 代码。例如，行 v_y = x/norm 检查 norm 是否为 0，否则引发 ZeroDivisionError。行 x = r_i[j, 0] 显示 Cython 检查索引是否在数组的边界内。您可能会注意到最后一行颜色非常强烈，通过检查代码我们可以看到这实际上是一个错误；代码引用了与函数结束相关的样板代码。

Cython 可以通过其编译器指令关闭这些检查以提高速度。有三种不同的方式来添加编译器指令：

• 使用装饰器或上下文管理器

• 在文件开头使用注释

• 使用 Cython 命令行选项

  注意

  您可以在官方文档中找到 Cython 编译器指令的完整列表，链接为 docs.cython.org/src/reference/compilation.html#compiler-directives

例如，要禁用数组的 "bounds" 检查，只需以下方式装饰一个函数 cython.boundscheck：

cimport cython

@cython.boundscheck(False)
def myfunction:
  # Code here

我们可以使用 cython.boundscheck 将代码块包装到上下文管理器中，如下所示：

with cython.boundscheck(False):
  # Code here

如果我们想要禁用整个模块的边界检查，我们可以在文件开头添加以下代码行：

# cython: boundscheck=False

要使用命令行选项更改指令，您可以使用 -X 如下所示：

$ cython -X boundscheck=True

我们现在可以尝试通过禁用 boundscheck 指令并启用 cdivision（这将禁用 ZeroDivisionError 的检查）来避免函数中的额外检查，如下面的代码所示：

cimport cython

@cython.boundscheck(False)
@cython.cdivision(True)
def c_evolve(double[:, :] r_i,double[:] ang_speed_i,double timestep,int nsteps):

如果我们再次查看注释视图，循环体完全为白色；我们从循环中移除了所有解释器的痕迹。然而，在以下情况下，我们没有获得性能提升：

In [3]: %timeit benchmark(100, 'cython')
100 loops, best of 3: 13.4 ms per loop

我们可以通过在文件中包含 profile=True 指令来使用 cProfile 对 Cython 代码进行性能分析。为了展示其用法，我们可以编写一个函数来计算两个坐标数组之间的 Chebyshev 距离。创建一个名为 cheb.py 的文件：

import numpy as np
from distance import chebyshev

def benchmark():
  a = np.random.rand(100, 2)
  b = np.random.rand(100, 2)
  for x1, y1 in a:
    for x2, y2 in b:
      chebyshev(x1, x2, y1, y2)

如果我们尝试以当前状态分析此脚本，我们将不会得到关于我们在 Cython 中实现的函数的任何统计信息。如果我们想了解max和min函数的配置文件指标，我们必须将profile=True选项添加到mathlib.pyx文件中，如下所示：

# cython: profile=True

cdef int max(int a, int b):
  # Code here

现在，我们可以使用 IPython 中的%prun来分析我们的脚本，如下所示：

In [2]: import cheb
In [3]: %prun cheb.benchmark()
     2000005 function calls in 2.066 seconds

  Ordered by: internal time

  ncalls tottime percall cumtime percall filename:lineno(function)
    1  1.664  1.664  2.066  2.066 cheb.py:4(benchmark)
 1000000  0.351  0.000  0.401  0.000 {distance.chebyshev}
 1000000  0.050  0.000  0.050  0.000 mathlib.pyx:2(max)
    2  0.000  0.000  0.000  0.000 {method 'rand' of 'mtrand.RandomState' objects}
    1  0.000  0.000  2.066  2.066 <string>:1(<module>)
    1  0.000  0.000  0.000  0.000 {method 'disable' of '_lsprof.Profiler' objects}

从输出中，我们可以看到max函数存在并且不是瓶颈。问题似乎出在benchmark函数上；问题很可能是 Python 循环的开销。在这种情况下，最佳策略是将循环重写为 NumPy 或移植代码到 Cython。

摘要

Cython 将把您程序的运行速度提升到另一个层次。与 C 语言相比，Cython 程序更容易维护，这得益于与 Python 的紧密集成以及可用性分析工具。

在本章中，我们介绍了 Cython 语言的基础知识以及如何通过添加静态类型来使我们的程序更快。我们还学习了如何处理 C 数组、NumPy 数组和内存视图。

我们通过重写关键的evolve函数优化了我们的粒子模拟器，获得了巨大的速度提升。最后，我们学会了如何使用注释视图快速定位与解释器相关的调用，以及如何为 Cython 脚本启用cProfile。

在下一章中，我们将学习并行处理的基础知识，并了解如何编写利用多个处理器的 Python 程序，以便您可以编写更快的程序并解决更大的问题。

第四章 并行处理

使用并行处理，你可以在不使用更快的处理器的情况下，在给定时间内增加程序可以完成的计算量。主要思想是将任务划分为许多子单元，并使用多个处理器独立解决它们。

包含多个核心（2、4、6、8、...）的 CPU 已成为技术界的普遍趋势。提高单个处理器的速度成本高昂且问题重重；而利用价格更低的多个核心处理器的并行能力是提高性能的可行途径。

并行处理让你能够处理大规模问题。科学家和工程师通常在超级计算机上运行并行代码——由大量标准处理器组成的庞大网络——以模拟庞大的系统。并行技术还可以利用图形芯片（一种针对并行化优化的硬件）。

Python 可以用于所有这些领域，使我们能够以简单和优雅的方式将并行处理应用于各种问题，开启无限可能的大门。

在本章中，我们将：

• 简要介绍并行处理的基本原理

• 使用 multiprocessing Python 库说明如何并行化简单问题

• 学习如何使用IPython parallel框架编写程序

• 使用 Cython 和 OpenMP 进一步优化我们的程序

并行编程简介

为了并行化一个程序，我们需要将问题划分为可以独立（或几乎独立）运行的子单元。

当子单元之间完全独立时，这样的问题被称为令人尴尬的并行。数组上的元素级操作是一个典型例子——操作只需要知道它当前处理的元素。另一个例子是我们的粒子模拟器——由于没有相互作用，每个粒子可以独立于其他粒子在时间上进化。令人尴尬的并行问题很容易实现，并且在并行架构上表现最优。

其他问题可能可以划分为子单元，但必须共享一些数据以执行它们的计算。在这些情况下，实现方式不太直接，并且由于通信成本可能导致性能问题。

我们将通过一个例子来说明这个概念。想象你有一个粒子模拟器，但这次粒子在特定距离内会吸引其他粒子（如图所示）。为了并行化这个问题，我们将模拟区域划分为区域，并将每个区域分配给不同的处理器。如果我们对系统进行一步进化，一些粒子将与相邻区域的粒子相互作用。为了执行下一次迭代，需要相邻区域的新粒子位置：

进程间的通信成本高昂，可能会严重阻碍并行程序的性能。在并行程序中处理数据通信存在两种主要方法：

• 共享内存

• 分布式内存

在共享内存中，子单元可以访问相同的内存空间。这种方法的优点是，你不需要显式地处理通信，因为从共享内存中写入或读取就足够了。然而，当多个进程同时尝试访问和更改相同的内存位置时，就会出现问题。应小心使用同步技术来避免此类冲突。

在分布式内存模型中，每个进程与其他进程完全分离，并拥有自己的内存空间。在这种情况下，进程间的通信是显式处理的。与共享内存相比，通信开销通常更昂贵，因为数据可能需要通过网络接口传输。

在共享内存模型中实现并行性的一个常见方法是 线程。线程是从进程派生出来的独立子任务，并共享资源，如内存。

Python 可以创建和处理线程，但由于 Python 解释器的设计，它们不能用来提高性能——一次只能允许一个 Python 指令运行。这种机制被称为 全局解释器锁（GIL）。其工作原理是，每次线程执行 Python 语句时，都会获取一个锁，这阻止了其他线程在锁释放之前运行。GIL 避免了线程之间的冲突，简化了 CPython 解释器的实现。尽管存在这种限制，但在可以释放锁的情况，例如耗时的 I/O 操作或 C 扩展中，线程仍然可以用来提供并发性。

可以通过使用进程而不是线程来完全避免 GIL。进程不共享相同的内存区域，彼此独立——每个进程都有自己的解释器。通过使用进程，我们将几乎没有缺点：进程间通信比共享内存效率低，但它更灵活且更明确。

多进程模块

标准的 multiprocessing 模块可以通过创建多个进程来快速并行化简单任务。它的接口易于使用，并包括一些用于处理任务提交和同步的实用工具。

进程和池类

你可以通过继承 multiprocessing.Process 来创建一个独立运行的进程。你可以扩展 __init__ 方法来初始化资源，并且可以通过实现 Process.run 方法来编写分配给子进程的代码部分。在以下代码中，我们定义了一个将等待一秒并打印其分配的 id 的进程：

import multiprocessing
import time

class Process(multiprocessing.Process):
    def __init__(self, id):
        super(Process, self).__init__()
        self.id = id

    def run(self):
        time.sleep(1)
        print("I'm the process with id: {}".format(self.id))

要创建进程，我们必须初始化我们的 Process 对象并调用 Process.start 方法。请注意，您不会直接调用 Process.run：调用 Process.start 将创建一个新的进程，并依次调用 Process.run 方法。我们可以在脚本末尾添加以下行来初始化并启动新进程：

if __name__ == '__main__':
    p = Process(0)
    p.start()

在 Process.start 之后的指令将立即执行，而不需要等待进程 p 完成。要等待任务完成，您可以使用 Process.join 方法，如下所示：

 if __name__ == '__main__':
    p = Process(0)
    p.start()
    p.join()

我们可以用相同的方式启动四个不同的进程，它们将并行运行。在一个串行程序中，总共需要四秒钟。由于我们并行运行，每个进程将同时运行，从而实现 1 秒的墙钟时间。在以下代码中，我们创建了四个进程并并行启动它们：

if __name__ == '__main__':
    processes = Process(1), Process(2), Process(3), Process(4)
    [p.start() for p in processes]

注意，并行进程的执行顺序是不可预测的，它最终取决于操作系统如何调度进程执行。您可以通过多次运行程序来验证此行为——每次运行的顺序都会不同。

multiprocessing 模块提供了一个方便的接口，使得将任务分配和分发到一组进程变得容易，multiprocessing.Pool 类。

multiprocessing.Pool 类创建了一组进程——称为 工作者——并允许通过 apply/apply_async 和 map/map_async 方法提交任务。

Pool.map 方法将函数应用于列表中的每个元素，并返回结果列表。它的用法与内置（串行）map 相当。

要使用并行映射，您首先需要初始化一个 multiprocessing.Pool 对象。它将工作进程的数量作为其第一个参数；如果没有提供，则该数量将与系统中的核心数相等。您可以通过以下方式初始化一个 multiprocessing.Pool 对象：

pool = multiprocessing.Pool()
pool = multiprocessing.Pool(processes=4)

让我们看看 Pool.map 的实际应用。如果您有一个计算数字平方的函数，您可以通过调用 Pool.map 并传递函数及其输入列表作为参数，将该函数映射到列表上，如下所示：

def square(x):
    return x * x

inputs = [0, 1, 2, 3, 4]
outputs = pool.map(square, inputs)

Pool.map_async 方法与 Pool.map 类似，但返回一个 AsyncResult 对象而不是实际的结果。当我们调用正常的 map 时，主程序的执行将停止，直到所有工作进程完成处理结果。使用 map_async，AsyncResult 对象将立即返回，而不会阻塞主程序，计算将在后台进行。然后我们可以通过使用 AsyncResult.get 方法在任何时候检索结果，如下所示：

outputs_async = pool.map_async(square, inputs)
outputs = outputs_async.get()

Pool.apply_async 将由单个函数组成的任务分配给一个工作进程。它接受函数及其参数，并返回一个 AsyncResult 对象。我们可以通过使用 apply_async 来获得类似于 map 的效果，如下所示：

results_async = [pool.apply_async(square, i) for i in range(100))]
results = [r.get() for r in results_async]

作为例子，我们将实现一个典型的、令人尴尬的并行程序：蒙特卡洛法估算π。

蒙特卡洛法估算π

想象我们有一个边长为 2 个单位的正方形；其面积将是 4 个单位。现在，我们在正方形内画一个半径为 1 个单位的圆，圆的面积将是 π * r²。通过将 r 的值代入前面的方程，我们得到圆的面积数值为 π * (1)² = π。你可以参考以下图示进行图形表示。

如果我们在该图形上随机射击很多点，一些点将落在圆内——我们将它们称为命中——而剩余的点——未命中——将位于圆外。蒙特卡洛方法的思想是圆的面积将与命中的数量成正比，而正方形的面积将与射击的总数成正比。为了得到π的值，只需将圆的面积（等于π）除以正方形的面积（等于 4）并解出π：

hits/total = area_circle/area_square = pi/4
pi = 4 * hits/total

我们程序中将采用的战略是：

• 在范围（-1，1）内生成大量的样本（x，y）数字

• 通过检查 x**2 + y**2 == 1 来测试这些数字是否位于圆内

我们首先编写一个串行版本并检查它是否工作。然后，我们可以编写并行版本。串行程序的实现如下：

import random

samples = 1000000
hits = 0

for i in range(samples):
    x = random.uniform(-1.0, 1.0)
    y = random.uniform(-1.0, 1.0)

    if x**2 + y**2 <= 1:
        hits += 1

pi = 4.0 * hits/samples

随着样本数量的增加，我们的近似精度将提高。你可以注意到每个循环迭代都是独立的——这个问题是令人尴尬的并行。

为了并行化此代码，我们可以编写一个名为 sample 的函数，它对应于单个命中-未命中检查。如果样本击中圆，则函数返回 1；否则返回 0。通过多次运行 sample 并汇总结果，我们将得到总命中数。我们可以通过 apply_async 在多个处理器上运行 sample 并以下列方式获取结果：

def sample():
    x = random.uniform(-1.0, 1.0)
    y = random.uniform(-1.0, 1.0)

    if x**2 + y**2 <= 1:
        return 1
    else:
        return 0

pool = multiprocessing.Pool()
results_async = [pool.apply_async(sample) for i in range(samples)]
hits = sum(r.get() for r in results_async)

我们可以将这两个版本封装在函数 pi_serial 和 pi_apply_async 中（你可以在 pi.py 文件中找到它们的实现）并按如下方式基准测试执行速度：

$ time python -c 'import pi; pi.pi_serial()'
real    0m0.734s
user    0m0.731s
sys    0m0.004s
$ time python -c 'import pi; pi.pi_apply_async()'
real    1m36.989s
user    1m55.984s
sys    0m50.386

如前基准测试所示，我们的第一个并行版本实际上削弱了我们的代码。原因是实际计算所需的时间与发送和分配任务到工作者的开销相比非常小。

为了解决这个问题，我们必须使开销与计算时间相比可以忽略不计。例如，我们可以要求每个工作者一次处理多个样本，从而减少任务通信开销。我们可以编写一个名为 sample_multiple 的函数，通过将问题分解为 10 个更密集的任务来修改我们的并行版本，如下面的代码所示：

def sample_multiple(samples_partial):
 return sum(sample() for i in range(samples_partial))

ntasks = 10
chunk_size = int(samples/ntasks)
pool = multiprocessing.Pool()
results_async = [pool.apply_async(sample_multiple, chunk_size)
 for i in range(ntasks)]
hits = sum(r.get() for r in results_async)

我们可以将这个功能封装在一个名为 pi_apply_async_chunked 的函数中，并按如下方式运行：

$ time python -c 'import pi; pi.pi_apply_async_chunked()'
real    0m0.325s
user    0m0.816s
sys    0m0.008s

结果要好得多；我们使程序的速度提高了不止一倍。你还可以注意到user指标大于real：总 CPU 时间大于总时间，因为同时使用了多个 CPU。如果你增加样本数量，你会注意到通信与计算的比率降低，从而提供更好的加速。

在处理令人尴尬的并行问题时，一切都很简单。但有时，你必须在进程之间共享数据。

同步和锁

即使multiprocessing使用进程（它们有自己独立的内存），它也允许你定义某些变量和数组作为共享内存。你可以通过使用multiprocessing.Value并将数据类型作为字符串传递（i表示整数，d表示双精度，f表示浮点数等）来定义一个共享变量。你可以通过value属性更新变量的内容，如下面的代码片段所示：

shared_variable = multiprocessing.Value('f')
shared_variable.value = 0

当使用共享内存时，你应该意识到并发访问。想象你有一个共享的整数变量，每个进程多次增加它的值。你可以定义一个进程类如下：

class Process(multiprocessing.Process):

    def __init__(self, counter):
        super(Process, self).__init__()
        self.counter = counter

    def run(self):
        for i in range(1000):
            self.counter.value += 1

你可以在主程序中初始化共享变量，并将其传递给4个进程，如下面的代码所示：

def main():
    counter = multiprocessing.Value('i', lock=True)
    counter.value = 0

    processes = [Process(counter) for i in range(4)]
    [p.start() for p in processes]
    [p.join() for p in processes] # processes are done
    print(counter.value)
main()

如果你运行这个程序（代码目录中的shared.py），你会注意到counter的最终值不是 4000，而是有随机值（在我的机器上它们在 2000 到 2500 之间）。如果我们假设算术是正确的，我们可以得出结论，并行化存在问题。

发生的情况是多个进程同时尝试访问相同的共享变量。这种情况最好通过以下图示来解释。在串行执行中，第一个进程读取（数字0），增加它，并写入新值（1）；第二个进程读取新值（1），增加它，并再次写入（2）。在并行执行中，两个进程同时读取值（0），增加它，并写入它（1），导致错误的结果。

为了解决这个问题，我们需要同步对这个变量的访问，以确保一次只有一个进程可以访问、增加并在共享变量上写入值。这个功能由multiprocessing.Lock类提供。锁可以通过acquire和release方法获取和释放，或者通过将锁用作上下文管理器来使用。当一个进程获取锁时，其他进程将无法获取它，直到锁被释放。

我们可以定义一个全局锁，并使用它作为上下文管理器来限制对计数器的访问，如下面的代码片段所示：

lock = multiprocessing.Lock()

class Process(multiprocessing.Process):

    def __init__(self, counter):
        super(Process, self).__init__()
        self.counter = counter

    def run(self):
        for i in range(1000):
            with lock: # acquire the lock
 self.counter.value += 1
 # release the lock

同步原语，如锁，对于解决许多问题是必不可少的，但你应该避免过度使用它们，因为它们可能会降低你程序的性能。

注意

multiprocessing 包含其他通信和同步工具，你可以参考官方文档以获取完整的参考信息：

docs.python.org/3/library/multiprocessing.html

IPython parallel

IPython 的强大之处不仅限于其高级壳。它的 parallel 包包括一个框架，用于在单核和多核机器上以及连接到网络的多个节点上设置和运行计算。IPython 很棒，因为它为并行计算增添了交互式特性，并为不同的通信协议提供了一个统一的接口。

要使用 IPython.parallel，你必须启动一组由 Controller（一个在客户端和引擎之间进行通信调度的实体）管理的 Engines（引擎）。这种方法与多进程完全不同；你将单独启动工作进程，它们将无限期地等待，监听来自客户端的命令。

要启动控制器和一组引擎（默认情况下，每个处理单元一个引擎），你可以使用 ipcluster 壳命令，如下所示：

$ ipcluster start

使用 ipcluster，你还可以通过编写自定义配置文件来设置多个节点，以在网络中分配你的计算。你可以参考以下网站上的官方文档以获取具体说明：

ipython.org/ipython-doc/dev/parallel/parallel_process.html

在启动控制器和引擎后，我们可以使用 IPython 壳来并行执行计算。IPython 提供了两个基本接口（或视图）：直接和基于任务的。

直接接口

直接接口允许你明确地向每个计算单元发送命令。该接口直观、灵活且易于使用，尤其是在交互式会话中使用时。

在启动引擎后，你必须在单独的壳中启动 IPython 会话来与之交互。通过创建一个客户端，你可以建立与控制器的连接。在以下代码中，我们导入 Client 类并创建一个实例：

In [1]: from IPython.parallel import Client
In [2]: rc = Client()

属性 Client.ids 将会给你一个表示可用引擎的整数列表，如下代码片段所示：

In [3]: rc.ids
Out[4]: [0, 1, 2, 3]

我们可以通过获取一个 DirectView 实例来向引擎发送命令。你可以通过索引 Client 实例或调用 DirectView.direct_view 方法来获取一个 DirectView 实例。以下代码展示了从先前创建的 Client 中获取 DirectView 实例的不同方法：

In [5]: dview = rc[0] # Select the first engine
In [6]: dview = rc[::2] # Select every other engine
In [7]: dview = rc[:] # Selects all the engines
In [8]: dview = rc.direct_view('all') # Alternative

你可以将引擎视为新的 IPython 会话。在最细粒度上，你可以通过使用 DirectView.execute 方法远程执行命令：

In [9]: dview.execute('a = 1')

命令将由每个引擎单独发送和执行。返回值将是一个 AsyncResult 对象，实际返回值可以通过 get 方法检索。

如以下代码所示，你可以使用DirectView.pull方法检索远程变量中的数据，并使用DirectView.push方法将数据发送到远程变量。DirectView类还支持方便的类似字典的接口：

In [10]: dview.pull('a').get() # equivalent: dview['a']
Out[10]: [1, 1, 1, 1]
In [11]: dview.push({'a': 2}) # equivalent: dview['a'] = 2

可以使用pickle模块发送和检索所有可序列化的对象。除此之外，还专门处理了如NumPy数组等数据结构以提高效率。

如果你发出一个会导致异常的语句，你将收到每个引擎中异常的摘要：

In [12]: res = dview.execute('a = *__*') # Invalid
In [13]: res.get()
[0:execute]:
  File "<ipython-input-3-945a473d5cbb>", line 1
    a = *__*
        ^
SyntaxError: invalid syntax

[1:execute]:
  File "<ipython-input-3-945a473d5cbb>", line 1
    a = *__*
        ^
SyntaxError: invalid syntax
[2: execute]:
...

应将引擎视为独立的 IPython 会话，并且必须通过网络同步导入和自定义函数。要导入一些库，包括本地和引擎中，可以使用DirectView.sync_imports上下文管理器：

with dview.sync_imports():
    import numpy
# The syntax import _ as _ is not supported

要将计算提交给引擎，DirectView提供了一些用于常见用例的实用工具，例如 map 和 apply。DirectView.map方法的工作方式类似于Pool.map_async，如下面的代码片段所示。你将函数映射到序列，返回一个AsyncResult对象。

In [14]: a = range(100)
In [15]: def square(x): return x * x
In [16]: result_async = dview.map(square, a)
In [17]: result = result_async.get()

IPython 通过DirectView.parallel装饰器提供了一个更方便的 map 实现。如果你在函数上应用装饰器，该函数现在将具有一个可以应用于序列的map方法。在以下代码中，我们将并行装饰器应用于square函数并将其映射到一系列数字上：

In [18]: @dview.parallel()
    ...: def square(x):
    ...:     return x * x
In [19]: square.map(range(100))

小贴士

要获取map的非阻塞版本，你可以使用DirectView.map_sync方法，或者将block=True选项传递给DirectView.parallel装饰器。

DirectView.apply方法的行为与Pool.apply_async不同。函数将在每个引擎上执行。例如，如果我们选择了四个引擎并应用square函数，该函数在每个引擎上执行一次，并返回四个结果，如下面的代码片段所示：

In [20]: def square(x):
            return x * x
In [21]: result_async = dview.apply(square, 2)
In [22]: result_async.get()
Out[22]: [4, 4, 4, 4]

DirectiView.remote装饰器允许你创建一个将在每个引擎上直接运行的函数。其用法如下：

In [23]: @dview.remote()
    ...: def square(x):
    ...:     return x * x
    ...:
In [24]: square(2)
Out[24]: [4, 4, 4, 4]

DirectView还提供了两种其他类型的通信方案：scatter和gather。

Scatter 将输入列表分配给引擎。想象一下你有四个输入和四个引擎；你可以使用DirectView.scatter在远程变量中分配这些输入，如下所示：

In [25]: dview.scatter('a', [0, 1, 2, 3])
In [26]: dview['a']
Out[26]: [[0], [1], [2], [3]]

Scatter 会尽可能均匀地分配输入，即使输入的数量不是引擎数量的倍数。以下代码展示了如何将 11 个计算处理成三个每批三个项目的批次和一个每批两个项目的批次：

In [13]: dview.scatter('a', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
In [14]: dview['a']
Out[14]: [[0, 1, 2], [3, 4, 5], [6, 7, 8], [9, 10]]

gather函数简单地检索分散的值并将它们合并。在以下代码片段中，我们将分散的结果合并回来：

In [17]: dview.gather('a').get()
Out[17]: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

我们可以使用 scatter 和 gather 函数来并行化我们的模拟之一。在我们的系统中，每个粒子与其他粒子是独立的，因此我们可以使用 scatter 和 gather 将粒子平均分配到可用的引擎之间，演化它们，并从引擎中获取粒子。

首先，我们必须设置引擎。ParticleSimulator 类应该对所有引擎可用。

请记住，引擎是在一个单独的进程中启动的，simul 模块应该可以被它们导入。你可以通过两种方式实现这一点：

• 通过在 simul.py 所在目录中启动 ipcluster

• 通过将那个目录添加到 PYTHONPATH

如果你正在使用代码示例，别忘了使用 setup.py 编译 Cython 扩展。

在以下代码中，我们创建粒子并获取一个 DirectView 实例：

from random import uniform
from simul import Particle
from IPython.parallel import Client

particles = [Particle(uniform(-1.0, 1.0),
                      uniform(-1.0, 1.0),
                      uniform(-1.0, 1.0)) for i in range(10000)]
rc = Client()
dview = rc[:]

现在，我们可以将粒子散射到远程变量 particle_chunk，使用 DirectView.execute 执行粒子演化并检索粒子。我们使用 scatter、execute 和 gather 来完成这项工作，如下面的代码所示：

dview.scatter('particle_chunk', particles, block=True)

dview.execute('from simul import ParticleSimulator')
dview.execute('simulator = ParticleSimulator(particle_chunk)')
dview.execute('simulator.evolve_cython(0.1)')

particles = dview.gather('particle_chunk', block=True)

我们现在可以封装并行版本，并将其与串行版本（参考文件 simul_parallel.py）进行基准测试，如下所示：

In [1]: from simul import benchmark
In [2]: from simul_parallel import scatter_gather
In [5]: %timeit benchmark(10000, 'cython')
1 loops, best of 3: 1.34 s per loop
In [6]: %timeit scatter_gather(10000)
1 loops, best of 3: 720 ms per loop

代码非常简单，给我们带来了 2 倍的速度提升，并且可以在任何数量的引擎上进行扩展。

基于任务的接口

IPython 有一个可以智能处理计算任务的接口。虽然这从用户的角度来看意味着一个不太灵活的接口，但它可以通过平衡引擎的负载和重新提交失败的作业来提高性能。在本节中，我们将介绍基于任务的接口中的 map 和 apply 函数。

任务接口由 LoadBalancedView 类提供，可以通过客户端使用 load_balanced_view 方法获得，如下所示：

In [1]: from IPython.parallel import Client
In [2]: rc = Client()
In [3]: tview = rc.load_balanced_view()

在这一点上，我们可以使用 map 和 apply 运行一些任务。LoadBalancedView 类与 multiprocessing.Pool 类似，任务由调度器提交和处理；在 LoadBalancedView 的情况下，任务分配基于特定时间引擎上的负载量，确保所有引擎都在工作，没有停机时间。

解释 DirectView 中的 apply 和 LoadBalancedView 之间的一个重要区别是有帮助的。对 DirectView.apply 的调用将在 每个 选定的引擎上运行，而对 LoadBalancedView.apply 的调用将调度一个 单个 任务到某个引擎。在前一种情况下，结果将是一个列表，在后一种情况下，它将是一个单个值，如下面的代码片段所示：

In [4]: dview = rc[:]
In [5]: tview = rc.load_balanced_view()
In [6]: def square(x):
   ...:     return x * x
   ...:
In [7]: dview.apply(square, 2).get()
Out[7]: [4, 4, 4, 4]
In [8]: tview.apply(square, 2).get()
Out[8]: 4

LoadBalancedView 也能够处理失败，并在满足某些条件时在引擎上运行任务。这个功能是通过一个依赖系统提供的。我们不会在本书中涵盖这个方面，但感兴趣的读者可以参考以下链接中的官方文档：

ipython.org/ipython-doc/rel-1.1.0/parallel/parallel_task.html

基于 OpenMP 的并行 Cython

Cython 提供了一个方便的接口，通过OpenMP执行共享内存并行处理。这使得你可以在 Cython 中直接编写非常高效的并行代码，而无需创建 C 包装器。

OpenMP 是一个用于编写多线程程序的规范，包括一系列 C 预处理器指令来管理线程；这些包括通信模式、负载均衡和同步功能。几个 C/C++和 Fortran 编译器（包括 GCC）实现了 OpenMP API。

让我们用一个小的例子来介绍 Cython 的并行功能。Cython 在cython.parallel模块中提供了一个基于 OpenMP 的简单 API。最简单的结构是prange：一个自动在多个线程中分配循环操作的构造。

首先，我们可以在hello_parallel.pyx文件中编写一个计算 NumPy 数组每个元素平方的程序。我们得到一个缓冲区作为输入，并通过填充输入数组元素的平方来创建一个输出数组。

以下代码片段显示了串行版本square_serial：

import numpy as np

def square_serial(double[:] inp):
    cdef int i, size
    cdef double[:] out
    size = inp.shape[0]
    out_np = np.empty(size, 'double')
    out = out_np

    for i in range(size):
        out[i] = inp[i]*inp[i]

    return out_np  

现在，我们可以通过将范围调用替换为prange来更改并行版本中的循环。有一个注意事项，你需要确保循环体中没有解释器。正如已经解释的，为了使用线程，我们需要释放 GIL，因为解释器调用会获取和释放 GIL，所以我们应避免它们。这样做失败会导致编译错误。

在 Cython 中，你可以使用nogil来释放 GIL，如下所示：

with nogil:
    for i in prange(size):
        out[i] = inp[i]*inp[i]

或者，你可以使用prange的方便选项nogil=True，这将自动将循环包装在nogil块中：

for i in prange(size, nogil=True):
    out[i] = inp[i]*inp[i]

注意

在prange块中尝试调用 Python 代码会导致错误。这包括赋值操作、函数调用、对象初始化等。要在prange块中包含此类操作（你可能想这样做以进行调试），你必须使用with gil语句重新启用 GIL：

for i in prange(size, nogil=True):
    out[i] = inp[i]*inp[i]
    with gil: 
 x = 0 # Python assignment

在这一点上，我们需要重新编译我们的扩展。我们需要更改setup.py以启用 OpenMP 支持。你必须使用distutils中的Extension类指定 GCC 选项-fopenmp，并将其传递给cythonize函数。以下代码显示了完整的setup.py文件：

from distutils.core import setup
from distutils.extension import Extension
from Cython.Build import cythonize

hello_parallel = Extension('hello_parallel',
 ['hello_parallel.pyx'],
 extra_compile_args=['-fopenmp'],
 extra_link_args=['-fopenmp'])

setup(
   name='Hello',
   ext_modules = cythonize(['cevolve.pyx', hello_parallel]),
)

现在我们知道了如何使用prange，我们可以快速并行化我们的ParticleSimulator的 Cython 版本。

在下面的代码中，我们可以查看 Cython 模块cevolve.pyx中包含的c_evolve函数，这是我们第二章“使用 NumPy 进行快速数组操作”中编写的：

def c_evolve(double[:, :] r_i,double[:] ang_speed_i,
             double timestep,int nsteps):

    # cdef declarations

    for i in range(nsteps):
        for j in range(nparticles):
            # loop body

我们首先要做的是反转循环的顺序；我们希望最外层的循环是并行循环，其中每个迭代都是独立的。由于粒子之间没有相互作用，我们可以安全地改变迭代的顺序，如下面的代码片段所示：

  for j in range(nparticles):
        for i in range(nsteps):

            # loop body

在这一点上，我们可以使用 prange 并行化循环，因为我们添加静态类型时已经移除了与解释器相关的调用，所以可以安全地应用 nogil 块，如下所示：

for j in prange(nparticles, nogil=True)

现在，我们可以将两个不同的版本封装成单独的函数，并对其进行计时，如下所示：

In [3]: %timeit benchmark(10000, 'openmp')
1 loops, best of 3: 599 ms per loop
In [4]: %timeit benchmark(10000, 'cython')
1 loops, best of 3: 1.35 s per loop

使用 OpenMP，我们能够通过更改一行代码与串行 Cython 版本相比获得显著的加速。

摘要

并行处理是提高程序速度或处理大量数据的有效方法。令人尴尬的并行问题是非常好的并行化候选者，并且导致实现简单和最佳扩展。

在本章中，我们介绍了 Python 并行编程的基础。我们学习了如何使用 Python 中已包含的工具轻松并行化程序。另一个更强大的并行处理工具是 IPython parallel。这个包允许你交互式地原型设计并行程序并有效地管理计算节点网络。最后，我们探讨了 Cython 和 OpenMP 的易于使用的多线程功能。

在本书的整个过程中，我们学习了设计、基准测试、分析和优化 Python 应用程序的最有效技术。NumPy 可以优雅地重写 Python 循环，如果还不够，你可以使用 Cython 生成高效的 C 代码。在最后阶段，你可以使用本章中介绍的工具轻松地并行化你的程序。