IPython-交互式计算和可视化秘籍第二版-二-IPython 交互式计算和可视化秘籍第二版（二）第三章：掌握笔记本

IPython 交互式计算和可视化秘籍第二版（二）

原文：annas-archive.org/md5/62516af4e05f6a3b0523d8aa07fd5acb

译者：飞龙

协议：CC BY-NC-SA 4.0

第三章：掌握笔记本

本章将涵盖以下主题：

在笔记本中使用 IPython 块教授编程
使用 nbconvert 将 IPython 笔记本转换为其他格式
在笔记本工具栏中添加自定义控件
自定义笔记本中的 CSS 样式
使用交互式小部件——笔记本中的钢琴
在笔记本中创建自定义 JavaScript 小部件——一个用于 pandas 的电子表格编辑器
从笔记本实时处理网络摄像头图像

介绍

在本章中，我们将了解笔记本的许多功能，包括 IPython 2.0 带来的交互式小部件。由于在前几章中我们仅了解了基本功能，这里将深入探讨笔记本的架构。

什么是笔记本？

笔记本于 2011 年发布，比 IPython 创建晚了十年。它的发展有着悠久且复杂的历史，费尔南多·佩雷斯在他的博客上对其进行了很好的总结，blog.fperez.org/2012/01/ipython-notebook-historical.html。受到数学软件如 Maple、Mathematica 或 Sage 的启发，笔记本极大地推动了 IPython 的流行。

通过将代码、文本、图像、图表、超链接和数学方程式混合在一个文档中，笔记本为交互式计算带来了可复现性。正确使用笔记本时，它能够彻底改变科学计算中的工作流程。在笔记本出现之前，人们不得不在文本编辑器和交互式提示之间来回切换；而现在，人们可以在一个统一的环境中保持专注。

笔记本不仅是一个工具，还是一个强大且健壮的架构。此外，这个架构大多是语言独立的，因此它不再仅限于 Python。笔记本定义了一套消息传递协议、API 和 JavaScript 代码，其他语言也可以使用这些协议和代码。实际上，我们现在已经看到非 Python 内核能够与笔记本进行交互，如 IJulia、IHaskell、IRuby 等。

在 2014 年 7 月的 SciPy 大会上，IPython 开发者甚至宣布他们决定将项目拆分为以下两部分：

新的Jupyter 项目将实现所有语言独立部分：笔记本、消息传递协议和整体架构。欲了解更多详情，请访问jupyter.org。
IPython 将作为 Python 内核的名称。

在本书中，我们并不做语义上的区分，而是使用 IPython 一词来指代整个项目（包括语言独立部分和 Python 内核）。

笔记本生态系统

笔记本被表示为JavaScript 对象 表示法（JSON）文档。JSON 是一种语言独立的、基于文本的文件格式，用于表示结构化文档。因此，笔记本可以被任何编程语言处理，并且可以转换为其他格式，如 Markdown、HTML、LaTeX/PDF 等。

一个围绕笔记本的生态系统正在建立，我们可以期待在不久的将来看到越来越多的使用案例。例如，Google 正在努力将 IPython 笔记本引入 Google Drive，以实现协作数据分析。此外，笔记本还被用来创建幻灯片、教学材料、博客文章、研究论文，甚至是书籍。实际上，本书完全是用笔记本写成的。

IPython 2.0 引入了交互式小部件（widgets）。这些小部件使 Python 和浏览器之间的联系更加紧密。我们现在可以创建实现 IPython 内核与浏览器之间双向通信的应用程序。此外，任何 JavaScript 交互式库原则上都可以集成到笔记本中。例如，D3.js JavaScript 可视化库现在被多个 Python 项目使用，以便为笔记本提供交互式可视化功能。我们可能会在不久的将来看到这些交互式功能的许多有趣应用。

IPython 笔记本的架构

IPython 实现了一个双进程模型，包含一个内核和一个客户端。客户端是提供给用户发送 Python 代码到内核的接口。内核执行代码并将结果返回给客户端以供显示。在读-评估-打印循环（REPL）术语中，内核实现了评估，而客户端则实现了读取和打印过程。

客户端可以是 Qt 控件（如果我们运行 Qt 控制台），也可以是浏览器（如果我们运行笔记本）。在笔记本中，内核一次接收完整的单元，因此并不认识笔记本的存在。包含笔记本的线性文档与底层内核之间有很强的解耦。这是一个很强的限制，可能会限制某些可能性，但它仍然带来了极大的简洁性和灵活性。

整个架构中的另一个基本假设是，每个笔记本至多只能连接一个内核。然而，IPython 3.0 提供了选择任何笔记本的语言内核的可能性。

在考虑笔记本的新使用场景时，牢记这些要点非常重要。

在笔记本中，除了 Python 内核和浏览器客户端外，还有一个基于Tornado的 Python 服务器（www.tornadoweb.org）。该进程为 HTML 基础的笔记本界面提供服务。

不同进程之间的所有通信过程都建立在ZeroMQ（或ZMQ）消息协议之上（zeromq.org）。笔记本通过WebSocket（一种基于 TCP 的协议，现代浏览器中已实现）与底层内核进行通信。

官方支持 IPython 2.x 笔记本的浏览器如下：

Chrome ≥ 13
Safari ≥ 5
Firefox ≥ 6

该笔记本还应在 Internet Explorer ≥ 10 上运行。这些要求本质上是 WebSocket 的要求。

将多个客户端连接到一个内核

在笔记本中，在一个单元格输入 %connect_info 会显示我们连接新客户端（例如 Qt 控制台）到底层内核所需的信息：

In [1]: %connect_info
{
  "stdin_port": 53978,
  "ip": "127.0.0.1", 
  "control_port": 53979, 
  "hb_port": 53980, 
  "signature_scheme": "hmac-sha256", 
  "key": "053...349", 
  "shell_port": 53976, 
  "transport": "tcp", 
  "iopub_port": 53977
}
Paste the above JSON code into a file, and connect with:
    $> ipython <app> --existing <file>
or, if you are local, you can connect with just:
    $> ipython <app> --existing kernel-6e0...b92.json
or even just:
    $> ipython <app> --existing
if this is the most recent IPython session you have started.

这里，<app> 是 console、qtconsole 或 notebook。

甚至可以让内核和客户端运行在不同的机器上。你可以在 IPython 文档中找到运行公共笔记本服务器的说明，详见 ipython.org/ipython-doc/dev/notebook/public_server.html#running-a-public-notebook-server。

笔记本中的安全性

有可能有人在 IPython 笔记本中插入恶意代码。由于笔记本中可能包含单元格输出中的隐藏 JavaScript 代码，因此理论上，当用户打开笔记本时，恶意代码有可能在不被察觉的情况下执行。

因此，IPython 2.0 引入了一个安全模型，在该模型中，笔记本中的 HTML 和 JavaScript 代码可以被标记为受信任或不受信任。用户生成的输出始终是受信任的。然而，当用户首次打开一个现有笔记本时，已存在的输出被视为不受信任。

安全模型基于每个笔记本中的加密签名。这个签名是通过每个用户拥有的私钥生成的。

你可以在接下来的部分中找到有关安全模型的更多参考资料。

参考文献

以下是一些关于笔记本架构的参考资料：

IPython 双进程模型，详见 ipython.org/ipython-doc/stable/overview.html#decoupled-two-process-model
该笔记本的文档，请参阅 ipython.org/ipython-doc/stable/interactive/notebook.html
笔记本中的安全性，详见 ipython.org/ipython-doc/dev/notebook/security.html
笔记本服务器，详见 ipython.org/ipython-doc/dev/interactive/public_server.html
IPython 消息协议，详见 ipython.org/ipython-doc/dev/development/messaging.html
关于如何为笔记本编写自定义内核的教程，详见 andrew.gibiansky.com/blog/ipython/ipython-kernels/

以下是一些（主要是实验性的）非 Python 语言内核，用于笔记本：

IJulia，请参阅 github.com/JuliaLang/IJulia.jl
IRuby，请参阅 github.com/isotope11/iruby
IHaskell，请参阅 github.com/gibiansky/IHaskell
IGo，可在github.com/takluyver/igo找到
IScala，可在github.com/mattpap/IScala找到

在 IPython 块中进行编程教学

IPython 笔记本不仅是科学研究和数据分析的工具，还是一款很棒的教学工具。在这个示例中，我们展示了一个简单有趣的 Python 库，用于教授编程概念：IPython Blocks（可在ipythonblocks.org找到）。这个库允许你或你的学生创建五颜六色的方块网格。你可以改变每个方块的颜色和大小，甚至可以让网格动起来。通过这个工具，你可以演示许多基本的技术概念。这个工具的可视化特点使得学习过程更加高效且富有吸引力。

在这个示例中，我们将主要执行以下任务：

用动画演示矩阵乘法
将图像显示为方块网格

这个示例部分灵感来源于nbviewer.ipython.org/gist/picken19/b0034ba7ec690e89ea79中的一个例子。

准备工作

你需要为这个示例安装 IPython Blocks。你只需在终端中输入pip install ipythonblocks即可。请注意，你也可以通过在 IPython 笔记本中在命令前加上!来执行这个命令。

In [1]: !pip install ipythonblocks

在这个示例的最后部分，你还需要安装 Pillow，相关文档可在pillow.readthedocs.org/en/latest/找到；更多说明请参见第十一章，图像和音频处理。如果你使用 Anaconda，可以在终端执行conda install pillow。

最后，你需要从本书网站下载Portrait数据集（github.com/ipython-books/cookbook-data）并在当前目录下解压。你也可以使用你自己的图像进行实验！

如何操作...

首先，我们导入一些模块，代码如下：

In [1]: import time
        from IPython.display import clear_output
        from ipythonblocks import BlockGrid, colors

现在，我们创建一个五列五行的方块网格，并将每个方块填充为紫色：

In [2]: grid = BlockGrid(width=5, height=5,
                         fill=colors['Purple'])
        grid.show()

如何操作...

我们可以通过二维索引来访问单独的方块。这演示了 Python 中的索引语法。我们还可以用:（冒号）来访问整个行或列。每个方块都是由 RGB 颜色表示的。这个库附带了一个便捷的颜色字典，将 RGB 元组分配给标准颜色名称，如下所示：
```
In [3]: grid[0,0] = colors['Lime']
        grid[-1,0] = colors['Lime']
        grid[:,-1] = colors['Lime']
        grid.show()
```

现在，我们将演示矩阵乘法，这是线性代数中的一个基本概念。我们将表示两个（n,n）矩阵，A（青色）和B（青柠色），并与C = A B（黄色）对齐。为了实现这一点，我们采用一个小技巧，创建一个大小为（2n+1,2n+1）的大白网格。矩阵A、B和C只是网格部分的视图。

In [4]: n = 5
        grid = BlockGrid(width=2*n+1, 
                         height=2*n+1, 
                         fill=colors['White'])
        A = grid[n+1:,:n]
        B = grid[:n,n+1:]
        C = grid[n+1:,n+1:]
        A[:,:] = colors['Cyan']
        B[:,:] = colors['Lime']
        C[:,:] = colors['Yellow']
        grid.show()

如何操作...

让我们转向矩阵乘法本身。我们对所有行和列进行循环，并突出显示在矩阵乘法过程中相乘的 A 和 B 中对应的行和列。我们结合使用 IPython 的 clear_output() 方法与 grid.show() 和 time.sleep()（暂停）来实现动画，如下所示：

In [5]: for i in range(n):
            for j in range(n):
                # We reset the matrix colors.
                A[:,:] = colors['Cyan']
                B[:,:] = colors['Lime']
                C[:,:] = colors['Yellow']
                # We highlight the adequate rows
                # and columns in red.
                A[i,:] = colors['Red']
                B[:,j] = colors['Red']
                C[i,j] = colors['Red']
                # We animate the grid in the loop.
                clear_output()
                grid.show()
                time.sleep(0.25)

如何操作...

最后，我们将使用 IPython Blocks 显示一张图片。我们使用 Image.open() 导入 JPG 图片，并通过 getdata() 获取数据，具体如下：
```
In [6]: from PIL import Image
        imdata = Image.open('data/photo.jpg').getdata()
```

现在，我们创建一个 BlockGrid 实例，设置适当的行数和列数，并将每个块的颜色设置为图片中相应像素的颜色。我们使用较小的块大小，并去除块之间的线条，如下所示：

In [7]: rows, cols = imdata.size
        grid = BlockGrid(width=rows, height=cols,
                         block_size=4, lines_on=False)
        for block, rgb in zip(grid, imdata):
            block.rgb = rgb
        grid.show()

如何操作...

还有更多...

正如本食谱中所演示的，笔记本是一个理想的教育平台，适用于所有层次的教育活动。

该库是在 IPython 2.0 引入交互式笔记本功能之前开发的。现在，我们可以期待更多交互式的发展。

使用 nbconvert 将 IPython 笔记本转换为其他格式

一个 IPython 笔记本会以 JSON 文本文件的形式保存。该文件包含笔记本的所有内容：文本、代码和输出。matplotlib 图形会以 base64 字符串的形式编码在笔记本内，导致笔记本文件既独立又可能较大。

注意

JSON 是一种人类可读、基于文本的开放标准格式，可以表示结构化数据。尽管源自 JavaScript，它是语言独立的。其语法与 Python 字典有一些相似之处。JSON 可以在多种语言中解析，包括 JavaScript 和 Python（Python 标准库中的 json 模块）。

IPython 提供了一个名为 nbconvert 的工具，可以将笔记本转换为其他格式：纯文本、Markdown、HTML、LaTeX/PDF，甚至是使用 reveal.js 库的幻灯片。你可以在 nbconvert 文档中找到有关不同支持格式的更多信息。

在本食谱中，我们将了解如何操作笔记本的内容，并将其转换为其他格式。

准备工作

你需要安装 pandoc，可以在 johnmacfarlane.net/pandoc/ 获取，这是一个用于将文件从一种标记语言转换为另一种标记语言的工具。

要将笔记本转换为 PDF，你需要一个 LaTeX 发行版，可以在 latex-project.org/ftp.html 找到。你还需要从本书网站下载 Notebook 数据集（github.com/ipython-books/cookbook-data），并将其解压到当前目录中。

在 Windows 上，你可能需要安装 pywin32 包。如果你使用 Anaconda，可以通过 conda install pywin32 来安装它。

如何操作...

让我们打开 data 文件夹中的 test 笔记本。笔记本只是一个普通的文本文件（JSON），所以我们以文本模式（r模式）打开它，如下所示：

In [1]: with open('data/test.ipynb', 'r') as f:
            contents = f.read()
        print(len(contents))
3787

这是test.ipynb文件的摘录：

{
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "name": "",
  "signature": "sha256:50db..."
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
...
     "source": [
      "# First chapter"
     ]
    },
  ...
   ],
   "metadata": {}
  }
 ]
}

现在我们已经将笔记本加载为字符串，让我们使用 json 模块进行解析，如下所示：
```
In [3]: import json
        nb = json.loads(contents)
```
让我们来看看笔记本字典中的键：
```
In [4]: print(nb.keys())
        print('nbformat ' + str(nb['nbformat']) + 
              '.' + str(nb['nbformat_minor']))
[u'nbformat', u'nbformat_minor', u'worksheets', u'metadata']
nbformat 3.0
```
注意

笔记本格式的版本在nbformat和nbformat_minor中指示。笔记本格式的向后不兼容更改预计将在未来的 IPython 版本中出现。此食谱已在 IPython 2.x 分支和笔记本格式 v3 中进行测试。
主要字段是worksheets，默认情况下只有一个。一个工作表包含单元格列表和一些元数据。worksheets字段在未来版本的笔记本格式中可能会消失。让我们来看看一个工作表的内容：
```
In [5]: nb['worksheets'][0].keys()
Out[5]: [u'cells', u'metadata']
```

每个单元格都有类型、可选的元数据、一些内容（文本或代码）、可能有一个或多个输出以及其他信息。让我们看看一个 Markdown 单元格和一个代码单元格：

In [6]: nb['worksheets'][0]['cells'][1]
Out[6]: {u'cell_type': u'markdown',
 u'metadata': {u'my_field': [u'value1', u'2405']},
               u'source': [u"Let's write ...:\n", ...]}
In [7]: nb['worksheets'][0]['cells'][2]
Out[7]: {u'cell_type': u'code',
         u'collapsed': False,
         u'input': [u'import numpy as np\n', ...],
         u'language': u'python',
         u'metadata': {},
         u'outputs': [
                      {u'metadata': {},
                       u'output_type': u'display_data',
                       u'png': u'iVB...mCC\n',
                       u'prompt_number': 1}]}

一旦被解析，笔记本将表示为一个 Python 字典。因此，在 Python 中操作它是非常方便的。在这里，我们按照如下方式计算 Markdown 和代码单元格的数量：

In [8]: cells = nb['worksheets'][0]['cells']
        nm = len([cell for cell in cells
                  if cell['cell_type'] == 'markdown'])
        nc = len([cell for cell in cells
                  if cell['cell_type'] == 'code'])
        print(("There are {nm} Markdown cells and "
               "{nc} code cells.").format(nm=nm, nc=nc))
There are 2 Markdown cells and 1 code cells.

让我们仔细看看包含 matplotlib 图形的单元格的图像输出：

In [9]: png = cells[2]['outputs'][0]['png']
        cells[2]['outputs'][0]
Out[9]: {u'metadata': {},
         u'output_type': u'display_data',
         u'png': u'iVBORwoAAAANSUhE...ErAAAElTkQmCC\n'}

通常，可以有零个、一个或多个输出。此外，每个输出可以有多个表示。在这里，matplotlib 图形有 PNG 表示（base64 编码的图像）和文本表示（图形的内部表示）。
现在，我们将使用 nbconvert 将我们的文本笔记本转换为其他格式。此工具可以从命令行使用。请注意，nbconvert 的 API 在未来版本中可能会发生变化。在这里，我们将笔记本转换为 HTML 文档，如下所示：
```
In [10]: !ipython nbconvert --to html data/test.ipynb
[NbConvertApp] Writing 187617 bytes to test.html
```
让我们在 <iframe> 中显示此文档（这是一个在笔记本中显示外部 HTML 文档的小窗口）：
```
In [11]: from IPython.display import IFrame
         IFrame('test.html', 600, 200)
```
我们还可以将笔记本转换为 LaTeX 和 PDF。为了指定文档的标题和作者，我们需要扩展默认的 LaTeX 模板。首先，我们创建一个名为mytemplate.tplx的文件，扩展 nbconvert 提供的默认article.tplx模板。我们指定作者和标题块的内容如下：
```
In [12]: %%writefile mytemplate.tplx
         ((*- extends 'article.tplx' -*))

         ((* block author *))
         \author{Cyrille Rossant}
         ((* endblock author *))

         ((* block title *))
         \title{My document}
         ((* endblock title *))
Writing mytemplate.tplx
```
然后，我们可以通过指定自定义模板来运行 nbconvert，如下所示：
```
In [13]: !ipython nbconvert --to latex --template mytemplate data/test.ipynb
         !pdflatex test.tex
[NbConvertApp] PDF successfully created
```
我们使用 nbconvert 将笔记本转换为 LaTeX，然后使用pdflatex（随 LaTeX 发行版一起提供）将 LaTeX 文档编译为 PDF。以下截图展示了笔记本的 PDF 版本：

它是如何工作的...

正如我们在这个食谱中看到的，.ipynb文件包含了笔记本的结构化表示。这种 JSON 文件可以在 Python 中轻松解析和操作。

nbconvert 是一个用于将笔记本转换为其他格式的工具。转换可以通过多种方式进行自定义。在这里，我们使用 jinja2 模板包扩展了现有模板。你可以在 nbconvert 的文档中找到更多信息。

还有更多...

有一个免费的在线服务，nbviewer，它允许我们在云中动态地将 IPython 笔记本渲染为 HTML。其理念是，我们提供一个原始笔记本（JSON 格式）的 URL 给 nbviewer，便可获得一个渲染后的 HTML 输出。nbviewer 的主页（nbviewer.ipython.org）包含了一些示例。

该服务由 IPython 开发者维护，并托管在 Rackspace（www.rackspace.com）上。

这里有一些更多的参考资料：

nbconvert 的文档，参见 ipython.org/ipython-doc/dev/interactive/nbconvert.html
一个 nbconvert 转换示例的列表，参见 github.com/ipython/nbconvert-examples
Wikipedia 上的 JSON 文章，网址为 en.wikipedia.org/wiki/JSON

在笔记本工具栏中添加自定义控件

HTML 笔记本的 CSS 和 JavaScript 可以通过~/.ipython/profile_default/static/custom目录中的文件进行自定义，其中~代表你的主目录，default是你的 IPython 配置文件。

在本食谱中，我们将使用这些自定义选项在笔记本工具栏中添加一个新的按钮，能够线性地重新编号所有代码单元。

如何操作...

首先，我们将直接在笔记本中注入 JavaScript 代码。这对于测试目的很有用，或者如果我们不希望更改是永久性的，JavaScript 代码将仅在该笔记本中加载。为了做到这一点，我们可以使用如下的 %%javascript 单元魔法：

In [1]: %%javascript
        // This function allows us to add buttons 
        // to the notebook toolbar.
        IPython.toolbar.add_buttons_group([
        {
            // The button's label.
            'label': 'renumber all code cells',

            // The button's icon.
            // See a list of Font-Awesome icons here:
            // http://fortawesome.github.io/Font-
            //                              Awesome/icons/
            'icon': 'icon-list-ol',

            // The callback function.
            'callback': function () {

                // We retrieve the lists of all cells.
                var cells = IPython.notebook.get_cells();

                // We only keep the code cells.
                cells = cells.filter(function(c)
                  {
                      return c instanceof IPython.CodeCell; 
                  })

                // We set the input prompt of all code 
                // cells.
                for (var i = 0; i < cells.length; i++) {
                    cells[i].set_input_prompt(i + 1);
                }
            }
        }]);

运行上述代码单元将在工具栏中添加一个按钮，如下图所示。点击此按钮会自动更新所有代码单元的提示编号。

添加重新编号工具栏按钮
为了使这些更改成为永久性的，也就是说，要在当前配置文件下的每个笔记本中添加此按钮，我们可以打开 ~/.ipython/profile_default/static/custom/custom.js 文件，并添加以下代码：
```
$([IPython.events]).on('app_initialized.NotebookApp',
    function(){
        // Copy the JavaScript code (in step 1) here.
    });
```
上述代码将自动加载到笔记本中，并且重新编号按钮将出现在当前配置文件下每个笔记本的顶部。

还有更多...

允许我们将按钮添加到笔记本工具栏的 IPython 笔记本 JavaScript API 在撰写本文时仍不稳定。它可能随时发生变化，且文档不完整。本食谱仅在 IPython 2.0 中测试过。你仍然可以在此页面找到一个非正式且部分的 API 文档：ipjsdoc.herokuapp.com。

我们应当期待未来会有更稳定的 API。

另请参见

在笔记本中自定义 CSS 样式配方

在笔记本中自定义 CSS 样式

在本配方中，我们展示了如何在笔记本界面和导出的 HTML 笔记本中自定义 CSS。

准备工作

你需要对 CSS3 有一定了解才能进行此配方。你可以在网上找到许多教程（参见本配方末尾的参考资料）。

你还需要从书籍网站下载Notebook数据集（链接：ipython-books.github.io），并将其解压到当前目录。

如何操作...

首先，我们创建一个新的 IPython 配置文件，以避免使默认配置文件杂乱无章，具体步骤如下：
```
In [1]: !ipython profile create custom_css
```

在 Python 中，我们检索此配置文件（~/.ipython）的路径以及custom.css文件的路径（默认为空）。

In [2]: dir = !ipython locate profile custom_css
        dir = dir[0]
In [3]: import os
        csspath = os.path.realpath(os.path.join(
                    dir, 'static/custom/custom.css'))
In [4]: csspath
Out[4]: '~\.ipython\profile_custom_css\static\
                                       custom\custom.css'

我们现在可以在这里编辑此文件。我们更改背景颜色、代码单元的字体大小、某些单元的边框，并在编辑模式中突出显示选定的单元。

In [5]: %%writefile {csspath}

        body {
            /* Background color for the whole notebook. */
            background-color: #f0f0f0;
        }

        /* Level 1 headers. */
        h1 {
            text-align: right;
            color: red;
        }

        /* Code cells. */
        div.input_area > div.highlight > pre {
            font-size: 10px;
        }

        /* Output images. */
        div.output_area img {
            border: 3px #ababab solid;
            border-radius: 8px;
        }

        /* Selected cells. */
        div.cell.selected {
            border: 3px #ababab solid;
            background-color: #ddd;
        }

        /* Code cells in edit mode. */
        div.cell.edit_mode {
            border: 3px red solid;
            background-color: #faa;
        }
Overwriting C:\Users\Cyrille\.ipython\profile_custom_css\static\custom\custom.css

使用custom_css配置文件打开笔记本（使用ipython notebook --profile=custom_css命令）会得到如下的自定义样式：

如何操作...

交互式笔记本界面中的自定义 CSS

我们还可以将此样式表与 nbconvert 一起使用。只需将笔记本转换为静态 HTML 文档，并将custom.css文件复制到当前目录。在这里，我们使用的是从书籍网站下载的测试笔记本（参见准备工作）：
```
In [6]: !cp {csspath} custom.css
        !ipython nbconvert --to html data/test.ipynb
[NbConvertApp] Writing 187617 bytes to test.html
```

这是该 HTML 文档的显示效果：

In [7]: from IPython.display import IFrame
        IFrame('test.html', 600, 650)

如何操作...

还有更多...

这里有一些关于 CSS 的教程和参考资料：

w3schools 上的 CSS 教程，链接：www.w3schools.com/css/DEFAULT…
Mozilla 开发者网络上的 CSS 教程，链接：developer.mozilla.org/en-US/docs/Web/Guide/CSS/Getting_started
Matthias Bussonnier 撰写的关于如何自定义笔记本 CSS 的博客文章，链接：nbviewer.ipython.org/github/Carreau/posts/blob/master/Blog1.ipynb

另请参见

在笔记本工具栏中添加自定义控件配方

使用交互式小部件——笔记本中的钢琴

从 IPython 2.0 开始，我们可以将交互式小部件放入笔记本中，以创建与 Python 内核交互的丰富 GUI 应用程序。IPython 提供了一组丰富的图形控件，如按钮、滑块和下拉菜单。我们可以完全控制它们的位置和外观。我们可以将不同的小部件组合成复杂的布局。我们甚至可以像在下一个配方中所展示的那样，从头开始创建我们自己的交互式小部件——在笔记本中创建自定义 Javascript 小部件——一个用于 pandas 的电子表格编辑器。

在本配方中，我们将展示 IPython 2.0+交互式小部件 API 提供的多种可能性。我们将在笔记本中创建一个非常基础的钢琴。

准备工作

你需要从本书网站下载 Piano 数据集（http://ipython-books.github.io）。该数据集包含在 archive.org 上获得的合成钢琴音符（CC0 1.0 公共领域授权）。它可以在 archive.org/details/SynthesizedPianoNotes 下载。

如何实现...

让我们导入一些模块，如下所示：

In [1]: import numpy as np
        import os
        from IPython.display import (Audio, display,
                                     clear_output)
        from IPython.html import widgets
        from functools import partial

要创建一个钢琴，我们将为每个音符绘制一个按钮。用户点击按钮时，播放相应的音符。这是通过如下方式显示 <audio> 元素来实现的：

In [2]: dir = 'data/synth'
In [3]: # This is the list of notes.
        notes = 'C,C#,D,D#,E,F,F#,G,G#,A,A#,B,C'.split(',')
In [4]: def play(note, octave=0):
            """This function displays an HTML Audio element
            that plays automatically when it appears."""
            f = os.path.join(dir, 
                 "piano_{i}.mp3".format(i=note+12*octave))
            clear_output()
            display(Audio(filename=f, autoplay=True))

我们将把所有按钮放置在一个容器小部件内。在 IPython 2.0 中，小部件可以按层次结构组织。一个常见的使用场景是将多个小部件组织在一个特定的布局中。在这里，piano 将包含 12 个按钮，代表 12 个音符：
```
In [5]: piano = widgets.ContainerWidget()
```
注意

用于创建容器小部件（如水平或垂直框）的 API 在 IPython 3.0 中发生了变化。有关更多详细信息，请参考 IPython 的文档。

我们创建第一个小部件：一个滑块控件，用于指定音阶（此处为 0 或 1）：

In [6]: octave_slider = widgets.IntSliderWidget()
        octave_slider.max = 1
        octave_slider

如何实现...

现在，我们创建按钮。有几个步骤。首先，我们为每个音符实例化一个 ButtonWidget 对象。然后，我们指定一个 callback() 函数，用于在给定的音阶（由当前的音阶滑块值决定）上播放相应的音符（由索引给出）。最后，我们设置每个按钮的 CSS，特别是白色或黑色的颜色。

In [7]: buttons = []
        for i, note in enumerate(notes):
            button = widgets.ButtonWidget(description=note)

            def on_button_clicked(i, _):
                play(i+1, octave_slider.value)

            button.on_click(partial(on_button_clicked, i))

            button.set_css({
                      'width': '30px', 
                      'height': '60px',
                      'padding': '0',
                      'color': 
                          ('black', 'white')['#' in note],
                      'background':
                          ('white', 'black')['#' in note],
                           'border': '1px solid black',
                           'float': 'left'})

            buttons.append(button)

最后，我们将所有小部件安排在容器内。piano 容器包含按钮，而主容器（container）包含滑块和钢琴。这可以通过以下方式实现：
```
In [8]: piano.children = buttons
In [9]: container = widgets.ContainerWidget()
        container.children = [octave_slider,
                              piano]
```
默认情况下，小部件在容器内垂直组织。在这里，音阶滑块将位于钢琴上方。在钢琴中，我们希望所有音符水平排列。我们通过将默认的 vbox CSS 类替换为 hbox 类来实现这一点。下图显示了 IPython 笔记本中的钢琴：
```
In [10]: display(container)
         piano.remove_class('vbox')
         piano.add_class('hbox')
```

它是如何工作的...

IPython 小部件由丰富的对象表示，这些对象在 Python 内核和浏览器之间共享。一个小部件包含特殊的属性，称为 trait 属性。例如，SliderWidget 的 value trait 属性动态且自动地与用户在笔记本滑块中选择的值相关联。

这个链接是双向的。在 Python 中更改这个属性会更新笔记本中的滑块。

小部件的位置由容器小部件和 CSS 类控制。你可以在文档中找到更多信息。

这种架构使得在笔记本中创建丰富的图形应用程序成为可能，且这些应用程序由 Python 代码支持。

还有更多内容...

小部件示例在 nbviewer.ipython.org/github/ipython/ipython/blob/master/examples/Interactive%20Widgets/Index.ipynb

另见

在笔记本中创建自定义的 JavaScript 小部件——一个用于 pandas 的电子表格编辑器 食谱

在笔记本中创建自定义的 JavaScript 小部件——一个用于 pandas 的电子表格编辑器

我们之前介绍了 IPython 笔记本 2.0 的新交互式功能。在这个食谱中，我们通过展示如何超越 IPython 2.0 提供的现有小部件，深入探讨了这个主题。具体来说，我们将创建一个自定义的基于 JavaScript 的小部件，它与 Python 内核进行通信。

具体来说，我们将在 IPython 笔记本中创建一个基本的交互式类 Excel 数据网格编辑器，兼容 pandas 的 DataFrame。从一个 DataFrame 对象开始，我们将能够在笔记本中的 GUI 内进行编辑。该编辑器基于 Handsontable JavaScript 库（handsontable.com）。也可以使用其他 JavaScript 数据网格编辑器。

准备工作

你需要 IPython 2.0+ 和 Handsontable JavaScript 库才能执行此食谱。以下是将此 JavaScript 库加载到 IPython 笔记本中的说明：

首先，访问 github.com/handsontable/jquery-handsontable/tree/master/dist。
然后，下载 jquery.handsontable.full.css 和 jquery.handsontable.full.js，并将这两个文件放入 ~\.ipython\profile_default\static\custom\。

在此文件夹中，在 custom.js 中添加以下行：

require(['/static/custom/jquery.handsontable.full.js']);

在此文件夹中，在 custom.css 中添加以下行：

@import "/static/custom/jquery.handsontable.full.css"

现在，刷新笔记本！

如何操作...

让我们导入以下几个函数和类：

In [1]: from IPython.html import widgets
        from IPython.display import display
        from IPython.utils.traitlets import Unicode

我们创建一个新小部件。value 特性将包含整个表格的 JSON 表示。由于 IPython 2.0 的小部件机制，这个特性将在 Python 和 JavaScript 之间进行同步。

In [2]: class HandsonTableWidget(widgets.DOMWidget):
            _view_name = Unicode('HandsonTableView',
                                 sync=True)
            value = Unicode(sync=True)

现在，我们为小部件编写 JavaScript 代码。负责同步的三个重要函数如下：

render 用于小部件初始化
update 用于 Python 到 JavaScript 的更新

handle_table_change 用于 JavaScript 到 Python 的更新

In [3]: %%javascript
var table_id = 0;
require(["widgets/js/widget"], function(WidgetManager){    
    // Define the HandsonTableView
    var HandsonTableView = IPython.DOMWidgetView.extend({

        render: function(){
            // Initialization: creation of the HTML elements
            // for our widget.

            // Add a <div> in the widget area.
            this.$table = $('<div />')
                .attr('id', 'table_' + (table_id++))
                .appendTo(this.$el);

            // Create the Handsontable table.
            this.$table.handsontable({
            });

        },

        update: function() {
            // Python --> Javascript update.

            // Get the model's JSON string, and parse it.
            var data = $.parseJSON(this.model.get('value'));

            // Give it to the Handsontable widget.
            this.$table.handsontable({data: data});

            return HandsonTableView.__super__.
                                         update.apply(this);
        },

        // Tell Backbone to listen to the change event 
        // of input controls.
        events: {"change": "handle_table_change"},

        handle_table_change: function(event) {
            // Javascript --> Python update.

            // Get the table instance.
            var ht = this.$table.handsontable('getInstance');

            // Get the data, and serialize it in JSON.
            var json = JSON.stringify(ht.getData());

            // Update the model with the JSON string.
            this.model.set('value', json);

            this.touch();
        },
    });

    // Register the HandsonTableView with the widget manager.
    WidgetManager.register_widget_view(
        'HandsonTableView', HandsonTableView);
});

现在，我们有了一个已经可以使用的同步表格小部件。然而，我们希望将其与 pandas 集成。为此，我们在 DataFrame 实例周围创建一个轻量级的包装器。我们创建了两个回调函数，用于将 pandas 对象与 IPython 小部件同步。GUI 中的更改将自动触发 DataFrame 的更改，但反之则不行。如果我们在 Python 中更改了 DataFrame 实例，我们需要重新显示小部件：

In [4]: from io import StringIO
        import numpy as np
        import pandas as pd
In [5]: class HandsonDataFrame(object):
            def __init__(self, df):
                self._df = df
                self._widget = HandsonTableWidget()
                self._widget.on_trait_change(
                           self._on_data_changed, 'value')
                self._widget.on_displayed(self._on_displayed)

            def _on_displayed(self, e):
                # DataFrame ==> Widget (upon initialization)
                json = self._df.to_json(orient='values')
                self._widget.value = json

            def _on_data_changed(self, e, val):
                # Widget ==> DataFrame (called every time the
                # user changes a value in the widget)
                buf = StringIO(val)
                self._df = pd.read_json(buf, orient='values')

            def to_dataframe(self):
                return self._df

            def show(self):
                display(self._widget)

现在，让我们测试一下所有这些！我们首先创建一个随机的 DataFrame 实例：

In [6]: data = np.random.randint(size=(3, 5),
                                 low=100, high=900)
        df = pd.DataFrame(data)
        df
Out[6]:      
352  201  859  322  352
326  519  848  802  642
171  480  213  619  192

我们将其包装在 HandsonDataFrame 中，并显示如下：
```
In [7]: ht = HandsonDataFrame(df)
        ht.show()
```

现在我们可以交互式地更改值，它们将在 Python 中相应地改变：

In [8]: ht.to_dataframe()
Out[8]:
352  201  859   322  352
326  519  848  1024  642
171  480  213   619  192

它是如何工作的……

让我们简要解释一下 IPython 2.0+ 中交互式 Python-JavaScript 通信的架构。

该实现遵循模型-视图-控制器（MVC）设计模式，这是图形界面应用程序中常用的模式。在后端（Python 内核）有一个模型，保存一些数据。在前端（浏览器），有一个或多个该模型的视图。视图与模型动态同步。当 Python 端的模型属性发生变化时，JavaScript 端的视图也会发生变化，反之亦然。我们可以实现 Python 和 JavaScript 函数来响应模型的变化。这些变化通常是由用户操作触发的。

在 Python 中，动态属性实现为特性（traits）。这些特殊的类属性在更新时会自动触发回调函数。在 JavaScript 中，使用的是 Backbone.js MVC 库。Python 和浏览器之间的通信是通过Comms完成的，这是一种在 IPython 中的特殊通信协议。

要创建一个新的小部件，我们需要创建一个继承自 DOMWidget 的类。然后，我们定义可以在 Python 和 JavaScript 之间同步的特性属性，如果传递 sync=True 给特性构造函数。我们可以注册回调函数来响应特性变化（无论是来自 Python 还是 JavaScript），使用 widget.on_trait_change(callback, trait_name)。callback() 函数可以具有以下签名之一：

callback()
callback(trait_name)
callback(trait_name, new_value)
callback(trait_name, old_value, new_value)

在 JavaScript 中，render() 函数会在初始化时创建单元格小部件区域的 HTML 元素。update() 方法允许我们响应后端（Python）中模型的变化。此外，我们还可以使用 Backbone.js 来响应前端（浏览器）中的变化。通过用 {"change": "callback"} 事件扩展小部件，我们告诉 Backbone.js 在 HTML 输入控件变化时，调用 callback() JavaScript 函数。这就是我们在这里响应用户触发的操作的方式。

还有更多……

以下是此概念验证的改进方式：

只同步变化，而不是每次都同步整个数组（这里使用的方法在大型表格上会比较慢）
避免在每次更改时重新创建一个新的 DataFrame 实例，而是在原地更新相同的 DataFrame 实例
支持命名列
隐藏包装器，即使得 DataFrame 在 notebook 中的默认丰富表示为 HandsonDataFrame
将一切实现为一个易于使用的扩展

以下是关于 IPython notebook 2.0+ 中小部件架构的一些参考资料：

关于自定义小部件的官方示例，网址：nbviewer.ipython.org/github/ipython/ipython/tree/master/examples/Interactive%20Widgets
Wikipedia 上的 MVC 模式，网址：en.wikipedia.org/wiki/Model%E2%80%93view%E2%80%93controller
Backbone.js，网址：backbonejs.org/
Backbone.js课程，网址：www.codeschool.com/courses/ana…
IPEP 21：小部件消息（comms），网址：github.com/ipython/ipython/wiki/IPEP-21%3A-Widget-Messages
IPEP 23：IPython 小部件，网址：github.com/ipython/ipython/wiki/IPEP-23%3A-Backbone.js-Widgets

另见

实时处理笔记本中的摄像头图像教程

从笔记本实时处理摄像头图像

在这个教程中，我们展示了如何让笔记本与 Python 内核进行双向通信。

具体来说，我们将使用 HTML5 的<video>元素从浏览器获取摄像头视频流，并通过 IPython notebook 2.0+的交互功能实时传递给 Python。然后，我们将在 Python 中使用边缘检测器（在 scikit-image 中实现）处理图像，并实时在笔记本中显示它。

这个教程的大部分代码来自 Jason Grout 的示例，网址：github.com/jasongrout/ipywidgets。

准备工作

你需要安装 Pillow 和 scikit-image 库来实现这个教程。（更多信息，请参见第十一章，图像和音频处理。）

你还需要一个支持 HTML5 捕获 API 的现代浏览器。你可以在dev.w3.org/2011/webrtc/editor/getusermedia.html找到规范。

如何实现...

我们需要导入几个模块，如下所示：

In [1]: from IPython.html.widgets import DOMWidget
        from IPython.utils.traitlets import (Unicode, Bytes,
                                             Instance)
        from IPython.display import display

        from skimage import io, filter, color
        import urllib
        import base64
        from PIL import Image
        from io import BytesIO # to change in Python 2
        import numpy as np
        from numpy import array, ndarray
        import matplotlib.pyplot as plt

我们定义了两个函数，用于将图像转换为和从 base64 字符串转换。这种转换是进程之间传递二进制数据的常见方法（在我们的案例中是浏览器和 Python 之间）：

In [2]: def to_b64(img):
            imgdata = BytesIO()
            pil = Image.fromarray(img)
            pil.save(imgdata, format='PNG')
            imgdata.seek(0)
            return urllib.parse.quote(
                        base64.b64encode(
                            imgdata.getvalue()))
In [3]: def from_b64(b64):
            im = Image.open(BytesIO(
                                    base64.b64decode(b64)))
            return array(im)

我们定义了一个 Python 函数，它将在实时处理网络摄像头图像时工作。它接受并返回一个 NumPy 数组。这个函数应用了一个边缘检测器，使用的是 scikit-image 中的roberts()函数，如下所示：
```
In [4]: def process_image(image):
            img = filter.roberts(image[:,:,0]/255.)
            return (255-img*255).astype(np.uint8)
```

现在，我们创建一个自定义小部件来处理浏览器和 Python 之间视频流的双向通信：

In [5]: 
class Camera(DOMWidget):
    _view_name = Unicode('CameraView', sync=True)

    # This string contains the base64-encoded raw
    # webcam image (browser -> Python).
    imageurl = Unicode('', sync=True)

    # This string contains the base64-encoded processed 
    # webcam image(Python -> browser).
    imageurl2 = Unicode('', sync=True)

    # This function is called whenever the raw webcam
    # image is changed.
    def _imageurl_changed(self, name, new):
        head, data = new.split(',', 1)
        if not data:
            return

        # We convert the base64-encoded string
        # to a NumPy array.
        image = from_b64(data)

        # We process the image.
        image = process_image(image)

        # We convert the processed image
        # to a base64-encoded string.
        b64 = to_b64(image)

        self.imageurl2 = 'data:image/png;base64,' + b64

下一步是为小部件编写 JavaScript 代码。由于代码较长，这里仅突出重要部分。完整代码可以在本书网站上找到：

In [6]: %%javascript 

var video        = $('<video>')[0];
var canvas       = $('<canvas>')[0];
var canvas2       = $('<img>')[0];
[...]

require(["widgets/js/widget"], function(WidgetManager){
    var CameraView = IPython.DOMWidgetView.extend({
        render: function(){
            var that = this;

            // We append the HTML elements.
            setTimeout(function() {
                that.$el.append(video).
                         append(canvas).
                         append(canvas2);}, 200);

            // We initialize the webcam.
            [...]

            // We initialize the size of the canvas.
            video.addEventListener('canplay', function(ev){
                if (!streaming) {
                  height = video.videoHeight / (
                      video.videoWidth/width);
                  video.setAttribute('width', width);
                  video.setAttribute('height', height);
                  [...]
                  streaming = true;
                }
            }, false);

            // Play/Pause functionality.
            var interval;
            video.addEventListener('play', function(ev){
                // We get the picture every 100ms.    
                interval = setInterval(takepicture, 100);
            })
            video.addEventListener('pause', function(ev){
                clearInterval(interval);
            })
            // This function is called at each time step.
            // It takes a picture and sends it to the model.
            function takepicture() {
                canvas.width = width; canvas.height = height;
                canvas2.width = width; canvas2.height = height;

                video.style.display = 'none';
                canvas.style.display = 'none';

                // We take a screenshot from the webcam feed and 
                // we put the image in the first canvas.
                canvas.getContext('2d').drawImage(video, 
                    0, 0, width, height);

                // We export the canvas image to the model.
                that.model.set('imageurl',
                               canvas.toDataURL('image/png'));
                that.touch();
            }
        },

        update: function(){
            // This function is called whenever Python modifies
            // the second (processed) image. We retrieve it and
            // we display it in the second canvas.
            var img = this.model.get('imageurl2');
            canvas2.src = img;
            return CameraView.__super__.update.apply(this);
        }
    });

    // Register the view with the widget manager.
    WidgetManager.register_widget_view('CameraView', 
                                       CameraView);
});

最后，我们创建并显示小部件，如下所示：
```
In [7]: c = Camera()
        display(c)
```

它是如何工作的…

让我们解释一下这个实现的原理。该模型有两个属性：从浏览器传入的（原始）图像和从 Python 传出的（处理后的）图像。每 100 毫秒，JavaScript 都会捕获网络摄像头画面（在<video>HTML 元素中），并将其复制到第一个画布上。画布图像被序列化为 base64 格式，并赋值给第一个模型属性。然后，Python 函数_imageurl_changed()被调用。图像被反序列化，通过 scikit-image 进行处理，并重新序列化。接着，第二个属性由 Python 修改，并设置为序列化后的处理图像。最后，JavaScript 中的update()函数会反序列化处理后的图像，并将其显示在第二个画布中。

还有更多…

通过从 Python 而非浏览器捕获网络摄像头图像，可以大大提高这个例子的速度。在这里，瓶颈可能来自于每次时间步长中从浏览器到 Python 以及反向传输的两次操作。

使用诸如OpenCV或SimpleCV等库从 Python 捕获网络摄像头图像会更高效。然而，由于这些库可能难以安装，因此让浏览器直接访问网络摄像头设备要简单得多。

另见

在笔记本中创建自定义 JavaScript 小部件——用于 pandas 的电子表格编辑器配方

第四章：性能分析与优化

本章将涵盖以下主题：

在 IPython 中评估语句所花费的时间
使用 cProfile 和 IPython 轻松分析代码
使用 line_profiler 逐行分析代码的性能
使用 memory_profiler 分析代码的内存使用情况
理解 NumPy 的内部机制，以避免不必要的数组复制
使用 NumPy 的跨步技巧
使用跨步技巧实现高效的滚动平均算法
在 NumPy 中进行高效的数组选择
使用内存映射处理超大的 NumPy 数组
使用 HDF5 和 PyTables 操作大数组
使用 HDF5 和 PyTables 操作大规模异构数据表

引言

尽管 Python 通常被认为是（有点不公平地）较慢的语言，但通过使用正确的方法，实际上可以实现非常好的性能。这就是本章和下一章的目标。本章将介绍如何评估（分析）程序变慢的原因，以及如何利用这些信息来优化代码，使其更加高效。下一章将讨论一些更高级的高性能计算方法，只有在本章中描述的方法不足以解决问题时才应采用。

本章的内容分为三个部分：

时间和内存性能分析：评估代码的性能
NumPy 优化：更高效地使用 NumPy，特别是在处理大数组时
内存映射与数组：为超大数组的外存计算实现内存映射技术，特别是使用 HDF5 文件格式

在 IPython 中评估语句所花费的时间

%timeit 魔法命令和 %%timeit 单元格魔法命令（适用于整个代码单元）允许你快速评估一个或多个 Python 语句所花费的时间。对于更全面的性能分析，你可能需要使用本章后续介绍的更高级方法。

如何实现...

我们将估算计算所有正整数的倒数平方和，直到给定的 n 所需的时间：

让我们定义 n：
```
In [1]: n = 100000
```

让我们在纯 Python 中计时这段计算：

In [2]: %timeit sum([1\. / i**2 for i in range(1, n)])
10 loops, best of 3: 131 ms per loop

现在，我们使用 %%timeit 单元格魔法命令来计时将相同的计算分成两行代码：

In [3]: %%timeit s = 0.
        for i in range(1, n):
            s += 1\. / i**2
10 loops, best of 3: 137 ms per loop

最后，让我们计时使用 NumPy 版本的计算：

In [4]: import numpy as np
In [5]: %timeit np.sum(1\. / np.arange(1., n) ** 2)
1000 loops, best of 3: 1.71 ms per loop

它是如何工作的...

%timeit 命令接受多个可选参数。其中一个参数是语句评估的次数。默认情况下，这个次数会自动选择，以确保 %timeit 命令在几秒钟内返回。然而，你也可以通过 -r 和 -n 参数直接指定这个次数。在 IPython 中输入 %timeit? 以获取更多信息。

%%timeit 单元格魔法命令还接受一个可选的设置语句（位于 %%timeit 的同一行），该语句会被执行，但不计时。所有在此语句中创建的变量都可以在单元格内部使用。

还有更多内容...

如果你不在 IPython 交互式会话中，可以使用 timeit.timeit()。这个函数定义在 Python 的 timeit 模块中，用于基准测试存储在字符串中的 Python 语句。IPython 的 %timeit 魔法命令是 timeit() 的一个方便封装，适用于交互式会话。有关 timeit 模块的更多信息，请参阅 docs.python.org/3/library/timeit.html。

另请参见

使用 cProfile 和 IPython 轻松分析代码 配方
逐行分析代码性能的 line_profiler 配方

使用 cProfile 和 IPython 轻松分析你的代码

%timeit 魔法命令通常很有用，但当你需要详细了解代码中哪些部分占用了最多执行时间时，它的功能略显有限。这个魔法命令更适用于基准测试（比较不同版本函数的执行时间），而不是性能分析（获取按函数细分的执行时间报告）。

Python 包含一个名为 cProfile 的性能分析器，可以将执行时间分解为所有调用函数的贡献。IPython 提供了在交互式会话中方便使用此工具的方法。

如何实现...

IPython 提供了 %prun 行魔法命令和 %%prun 单元格魔法命令，可以轻松地分析一行或多行代码的性能。%run 魔法命令也接受 -p 标志，用于在性能分析器的控制下运行 Python 脚本。这些命令有许多选项，你可能希望查看它们的文档，可以通过 %prun? 和 %run? 进行查询。

在这个示例中，我们将分析一个从原点开始的随机漫步数值模拟。我们将在第十三章中更详细地介绍这些类型的模拟，随机动力学系统。

让我们导入 NumPy 和 matplotlib：

In [1]: import numpy as np
        import matplotlib.pyplot as plt
In [2]: %matplotlib inline

让我们创建一个生成随机 +1 和 -1 值的函数，并将其存储在数组中：

In [3]: def step(*shape):
            # Create a random n-vector with +1 or -1
            # values.
            return 2 * (np.random.random_sample(shape) 
                        < .5) - 1

现在，让我们在一个以 %%prun 开头的单元格中编写模拟代码，以便分析整个模拟过程的性能。各种选项允许我们将报告保存到文件中，并按累计时间对前 10 个结果进行排序。我们将在 原理介绍 部分更详细地解释这些选项。

In [4]: %%prun -s cumulative -q -l 10 -T prun0
        n = 10000
        iterations = 50
        x = np.cumsum(step(iterations, n), axis=0)
        bins = np.arange(-30, 30, 1)
        y = np.vstack([np.histogram(x[i,:], bins)[0]
                       for i in range(iterations)])

性能分析报告已保存为名为 prun0 的文本文件。让我们展示一下它（以下输出是经过简化的版本，以适应本页面）：

In [5]: print(open('prun0', 'r').read())
             2960 function calls in 0.075 seconds
Ordered by: cumulative time
ncalls  cumtime  percall function
    50    0.037    0.001 histogram
     1    0.031    0.031 step
    50    0.024    0.000 sort
     1    0.019    0.019 rand
     1    0.005    0.005 cumsum

在这里，我们观察了在代码中直接或间接涉及的不同函数的执行时间。

如果我们将模拟的迭代次数从 50 增加到 500，那么运行相同的模拟，我们将得到以下结果：

   29510 function calls in 1.359 seconds
   ncalls  cumtime  percall function
      500    0.566    0.001 histogram
        1    0.388    0.388 cumsum
        1    0.383    0.383 step
      500    0.339    0.001 sort
        1    0.241    0.241 rand

我们可以观察到，迭代次数对涉及的函数（特别是 cumsum 函数）的相对性能开销有很大影响。

原理介绍...

Python 的性能分析器会生成关于我们代码执行时间的详细报告，按函数进行分类。在这里，我们可以观察到 histogram、cumsum、step、sort 和 rand 函数的调用次数，以及在代码执行过程中这些函数的总时间。内部函数也会被分析。对于每个函数，我们会得到总调用次数、总时间和累积时间，以及每次调用的对应值（通过 ncalls 除以总值）。总时间表示解释器在某个函数中停留的时间，不包括在调用子函数时所花费的时间。累积时间类似，但包括在调用子函数时所花费的时间。文件名、函数名和行号会显示在最后一列。

%prun 和 %%prun 魔法命令接受多个可选参数（输入 %prun? 查看详细信息）。在示例中，-s 允许我们按特定列排序报告，-q 用于抑制（抑制）分页器输出（当我们想将输出整合到笔记本中时很有用），-l 用于限制显示的行数或按函数名筛选结果（当我们关注某个特定函数时非常有用），-T 用于将报告保存为文本文件。此外，我们还可以选择使用 -D 保存（转储）二进制报告到文件中，或者使用 -r 在 IPython 中返回报告。这个类似数据库的对象包含所有分析信息，可以通过 Python 的 pstats 模块进行分析。

注意

每个性能分析器都有其自身的开销，可能会影响分析结果（探测效应）。换句话说，一个被分析过的程序可能比未分析的程序运行得慢得多。这一点需要记住。

“过早的优化是万恶之源”

正如 Donald Knuth 的名言所示，过早地优化代码通常被认为是一种不良实践。代码优化应仅在真正需要时进行，也就是说，当代码在正常情况下真的运行得很慢时。此外，我们应当准确知道需要优化代码的地方；通常，执行时间的大部分来自于代码的相对小部分。了解这一点的唯一方法是对代码进行性能分析；优化永远不应在没有初步分析的情况下进行。

提示

我曾经处理一些相当复杂的代码，速度比预期慢。我以为我对问题的原因和如何解决有相当好的想法。解决方案将涉及对代码的重大更改。幸运的是，我首先对我的代码进行了性能分析，只是为了确保。我的诊断似乎完全错误；我在某处错误地写成了 max(x) 而不是 np.max(x)，其中 x 是一个非常大的向量。调用的是 Python 的内置函数，而不是 NumPy 为数组高度优化的例程。如果我没有对代码进行性能分析，我可能会永远错过这个错误。程序运行得非常好，只是慢了 150 倍！

有关编程优化的更一般建议，请参阅 en.wikipedia.org/wiki/Program_optimization。

还有更多...

在 IPython 中对代码进行性能分析特别简单（尤其在笔记本中），正如我们在本方法中所见。但是，从 IPython 执行我们需要分析的代码可能是不可取或困难的（例如 GUI）。在这种情况下，我们可以直接使用 cProfile。这比在 IPython 中稍微复杂一些。

首先，我们调用以下命令：
```
$ python -m cProfile -o profresults myscript.py
```
文件 profresults 将包含 myscript.py 的性能分析结果的转储。
然后，我们从 Python 或 IPython 执行以下代码，以以人类可读的形式显示性能分析结果：
```
import pstats
p = pstats.Stats('profresults')
p.strip_dirs().sort_stats("cumulative").print_stats()
```

探索 cProfile 和 pstats 模块的文档，以了解您可以对性能分析报告执行的所有分析。

提示

位于 github.com/rossant/easy_profiler 的存储库包含一个简单的命令行工具，可帮助分析 Python 脚本的性能。

有一些 GUI 工具可用于探索和可视化性能分析会话的输出。例如，RunSnakeRun 允许您在 GUI 程序中查看性能分析结果。

以下是一些参考资料：

cProfile 和 pstats 的文档，可在 docs.python.org/3/library/profile.html 获取
RunSnakeRun，在 www.vrplumber.com/programming… 上
Python 的性能分析工具，可在 blog.ionelmc.ro/2013/06/08/python-profiling-tools/ 获取

另请参阅

使用 line_profiler 逐行分析您的代码性能 方法

使用 `line_profiler` 逐行分析您的代码性能

Python 的原生 cProfile 模块和相应的 %prun 魔术将代码的执行时间逐个函数地分解。有时，我们可能需要更细粒度的代码性能分析，以逐行报告。这样的报告可能比 cProfile 的报告更易读。

要按行分析代码，我们需要一个名为line_profiler的外部 Python 模块，由 Robert Kern 创建，模块可从pythonhosted.org/line_profiler/获得。在本教程中，我们将演示如何在 IPython 中使用该模块。

准备好

要安装line_profiler，在终端中输入pip install line_profiler，或在 IPython 中输入!pip install line_profiler（你需要一个 C 编译器）。

在 Windows 上，你可以使用 Chris Gohlke 提供的非官方包，下载地址为www.lfd.uci.edu/~gohlke/pyt…。

怎么做...

我们将逐行分析与上一教程相同的模拟代码：

首先，让我们导入 NumPy 和随包一起提供的line_profiler IPython 扩展模块：
```
In [1]: import numpy as np
In [2]: %load_ext line_profiler
```

这个 IPython 扩展模块提供了一个%lprun魔法命令，用于逐行分析 Python 函数。它在函数定义在文件中而不是在交互式命名空间或笔记本中时效果最好。因此，在这里，我们将代码写入 Python 脚本，并使用%%writefile单元魔法：

In [3]: %%writefile simulation.py
        import numpy as np 
        def step(*shape):
            # Create a random n-vector with +1 or -1 
            # values.
            return (2 * (np.random.random_sample(shape) 
                         < .5) - 1)
        def simulate(iterations, n=10000): 
            s = step(iterations, n)
            x = np.cumsum(s, axis=0)
            bins = np.arange(-30, 30, 1)
            y = np.vstack([np.histogram(x[i,:], bins)[0] 
                           for i in range(iterations)])
            return y

现在，让我们将这个脚本导入到交互式命名空间中，以便执行和分析我们的代码：
```
In [4]: import simulation
```
我们在行级分析器的控制下执行函数。需要分析的函数必须在%lprun魔法命令中明确指定。我们还将报告保存在一个文件中，命名为lprof0：
```
In [5]: %lprun -T lprof0 -f simulation.simulate simulation.simulate(50)
```

让我们展示报告（以下输出是经过精简的版本，以适应页面）：

In [6]: print(open('lprof0', 'r').read())
File: simulation.py
Function: simulate at line 7
Total time: 0.114508 s
Line #   % Time  Line Contents
     7           def simulate(iterations, n=10000):
     8    36.3       s = step(iterations, n)
     9     5.6       x = np.cumsum(s, axis=0)
    10     0.1       bins = np.arange(-30, 30, 1)
    11    58.1       y = np.vstack([np.histogram(...)])
    12     0.0       return y

如果我们用比之前多 10 倍的迭代次数（simulation.simulate(500)）执行相同的分析，我们会得到如下报告：

Total time: 1.28704 s
     7           def simulate(iterations, n=10000):
     8    29.2       s = step(iterations, n)
     9    30.9       x = np.cumsum(s, axis=0)
    10     0.0       bins = np.arange(-30, 30, 1)
    11    39.9       y = np.vstack([np.histogram(...)])
    12     0.0       return y

它是如何工作的...

%lprun命令接受一个 Python 语句作为其主要参数。需要分析的函数必须通过-f明确指定。其他可选参数包括-D、-T和-r，它们的工作方式类似于%prun魔法命令的对应参数。

line_profiler模块显示每一行分析函数所花费的时间，可以以计时单位或总执行时间的分数形式显示。当我们在查找代码热点时，这些详细信息至关重要。

还有更多内容...

与上一教程一样，可能需要对一个独立的 Python 程序运行逐行分析器，该程序无法从 IPython 轻松启动。这个过程稍显复杂。

我们从github.com/rkern/line_profiler/blob/master/kernprof.py下载kernprof文件，并将其保存在代码的目录中。
在代码中，我们用@profile装饰器来装饰我们希望分析的函数。我们需要在分析会话结束后删除这些装饰器，因为如果代码正常执行（即不在行级分析器的控制下），它们会引发NameError异常：
```
@profile
def thisfunctionneedstobeprofiled():
    pass
```
提示

参见stackoverflow.com/questions/18229628/python-profiling-using-line-profiler-clever-way-to-remove-profile-statements链接，了解一种巧妙的方法来移除配置文件语句。
我们在终端中执行以下命令：
```
python -m kernprof -l -v myscript.py > lprof.txt
```
将执行myscript.py脚本，并将报告保存到lprof.txt中。

提示

github.com/rossant/easy_profiler上的代码库提供了一个稍微简化的逐行分析工具使用方法。

跟踪 Python 程序逐步执行过程

我们还将讨论跟踪工具，它们对于性能分析、调试程序或用于教育目的非常有用。

Python 的trace模块允许我们跟踪 Python 代码的程序执行。这在深入调试和性能分析过程中非常有用。我们可以跟踪 Python 解释器执行的所有指令序列。有关 trace 模块的更多信息，请访问docs.python.org/3/library/trace.html。

此外，在线 Python Tutor 是一个在线交互式教育工具，帮助我们逐步理解 Python 解释器在执行程序源代码时的操作。在线 Python Tutor 可通过pythontutor.com/访问。

另见

使用 cProfile 和 IPython 轻松进行代码性能分析的技巧
使用 memory_profiler 分析代码的内存使用情况的技巧

使用 memory_profiler 分析代码的内存使用情况

前一篇配方中描述的方法是关于 CPU 时间的性能分析。这可能是代码性能分析中最显著的因素。然而，内存也是一个关键因素。例如，运行np.zeros(500000000)很可能会立即崩溃你的计算机！这个命令可能会分配超过系统可用内存的内存；你的计算机会在几秒钟内进入无响应状态。

编写内存优化代码并不简单，但能显著提升程序的运行速度。尤其在处理大型 NumPy 数组时，这一点尤为重要，正如我们将在本章后面看到的那样。

在这个配方中，我们将介绍一个简单的内存分析工具。这个库，毫不奇怪地叫做memory_profiler，由 Fabian Pedregosa 创建。它的使用方式与line_profiler非常相似，且可以方便地从 IPython 中使用。你可以从pypi.python.org/pypi/memory_profiler下载它。

准备工作

你可以通过pip install memory_profiler来安装memory_profiler。

在 Windows 上，您还需要 psutil，它可以在 pypi.python.org/pypi/psutil 上找到。您可以使用 pip install psutil 安装，或者从 code.google.com/p/psutil/ 下载该包。您也可以从 Chris Gohlke 的网页 www.lfd.uci.edu/~gohlke/pyt… 下载安装程序。

本方法中的示例是前一个方法的延续。

如何操作...

假设仿真代码已经如前一个方法中所示加载，我们加载内存分析器的 IPython 扩展：
```
In [9]: %load_ext memory_profiler
```

现在，让我们在内存分析器的控制下运行代码：

In [10]: %mprun -T mprof0 -f simulation.simulate simulation.simulate(50)

让我们展示结果：

In [11]: print(open('mprof0', 'r').read())
Filename: simulation.py
Line #    Mem usage  Increment   Line Contents
     7    39.672 MB  0.000 MB   def simulate(...):
     8    41.977 MB  2.305 MB       s = step(iterations, n)
     9    43.887 MB  1.910 MB       x = np.cumsum(...)
    10    43.887 MB  0.000 MB       bins = np.arange(...)
    11    43.887 MB  0.000 MB       y = np.vstack(...)
    12    43.887 MB  0.000 MB       return y

最后，这是进行 500 次迭代的报告：

Line #    Mem usage Increment   Line Contents
     7    40.078 MB  0.000 MB   def simulate(...):
     8    59.191 MB 19.113 MB       s = step(iterations, n)
     9    78.301 MB 19.109 MB       x = np.cumsum(...)
    10    78.301 MB  0.000 MB       bins = np.arange(...)
    11    78.301 MB  0.000 MB       y = np.vstack(...)
    12    78.301 MB  0.000 MB       return y

它是如何工作的...

memory_profiler 包检查每行代码的内存使用情况。增量列帮助我们发现代码中分配大量内存的地方。当处理数组时，这一点尤其重要。不必要的数组创建和复制会显著减慢程序速度。我们将在接下来的几个方法中解决这个问题。

还有更多...

我们可以在没有 IPython 的情况下使用 memory_profiler，也可以在 IPython 中对单个命令进行快速内存基准测试。

在独立的 Python 程序中使用 memory_profiler

在独立的 Python 程序中使用内存分析器与使用 line_profiler 类似，但稍微简单一些。

首先，在我们的 Python 脚本中，我们通过 @profile 装饰器标记我们希望分析的函数。
然后，我们运行：
```
$ python -m memory_profiler myscript.py > mprof.txt
```
分析报告将保存在 myprof.txt 中。

在 IPython 中使用 %memit 魔法命令

memory_profiler 的 IPython 扩展还附带了一个 %memit 魔法命令，让我们可以基准测试单个 Python 语句所使用的内存。这里是一个简单的例子：

In [14]: %memit np.random.randn(1000, 1000)
maximum of 1: 46.199219 MB per loop

其他工具

还有其他工具可以监控 Python 程序的内存使用情况，特别是 Guppy-PE (guppy-pe.sourceforge.net/)、PySizer (pysizer.8325.org/) 和 Pympler (code.google.com/p/pympler/)。与 IPython 及 Python 的自省功能结合使用时，这些工具允许您分析命名空间或特定对象的内存使用情况。

另见

逐行分析代码并使用 line_profiler 方法
理解 NumPy 的内部机制以避免不必要的数组复制 方法

理解 NumPy 的内部机制以避免不必要的数组复制

使用 NumPy，我们可以显著提高性能，特别是当我们的计算遵循 单指令多数据 (SIMD) 模式时。然而，也有可能无意中写出非优化的 NumPy 代码。

在接下来的几个案例中，我们将看到一些可以帮助我们编写优化 NumPy 代码的技巧。在这个案例中，我们将看到如何避免不必要的数组复制，从而节省内存。在这方面，我们需要深入了解 NumPy 的内部实现。

准备工作

首先，我们需要一种方法来检查两个数组是否共享相同的底层数据缓冲区。我们可以定义一个返回底层数据缓冲区内存位置的函数id()：

def id(x):
    # This function returns the memory
    # block address of an array.
    return x.__array_interface__['data'][0]

两个具有相同数据位置（由id返回的）数组共享相同的底层数据缓冲区。然而，只有在数组具有相同偏移量（意味着它们的第一个元素相同）时，才会发生这种情况。具有不同偏移量的共享数组会有稍微不同的内存位置，下面的示例说明了这一点：

In [1]: id(a), id(a[1:])
Out[1]: (71211328, 71211336)

在接下来的几个案例中，我们将确保使用相同偏移量的数组。以下是一个更通用且可靠的解决方案，用于判断两个数组是否共享相同的数据：

In [2]: def get_data_base(arr):
            """For a given Numpy array, finds the
            base array that "owns" the actual data."""
            base = arr
            while isinstance(base.base, np.ndarray):
                base = base.base
            return base

        def arrays_share_data(x, y):
            return get_data_base(x) is get_data_base(y)

In [3]: print(arrays_share_data(a,a.copy()),
              arrays_share_data(a,a[1:]))
False True

感谢 Michael Droettboom 指出这一点并提出这种替代解决方案。

如何做到...

使用 NumPy 数组进行计算可能涉及内存块之间的内部复制。这些复制并非总是必要的，如果没有必要，应避免它们，正如我们将在以下提示中看到的：

我们有时需要对数组进行复制；例如，如果我们需要在保留原始副本的情况下操作数组：
```
In [3]: a = np.zeros(10); aid = id(a); aid
Out[3]: 65527008L
In [4]: b = a.copy(); id(b) == aid
Out[4]: False
```

数组计算可以涉及就地操作（以下代码中的第一个示例：数组被修改）或隐式复制操作（第二个示例：创建了一个新数组）：

In [5]: a *= 2; id(a) == aid
Out[5]: True
In [6]: a = a*2; id(a) == aid
Out[6]: False

确保选择你实际需要的操作类型。隐式复制操作显著较慢，如下所示：

In [7]: %%timeit a = np.zeros(10000000)
        a *= 2
10 loops, best of 3: 23.9 ms per loop
In [8]: %%timeit a = np.zeros(10000000)
        a = a*2
10 loops, best of 3: 77.9 ms per loop

重塑数组可能会或可能不会涉及复制。原因将在*它是如何工作的...*部分解释。例如，重塑一个二维矩阵不会涉及复制，除非它被转置（或者更一般地说，非连续）：
```
In [9]: a = np.zeros((10, 10)); aid = id(a)
In [10]: b = a.reshape((1, -1)); id(b) == aid
Out[10]: True
In [11]: c = a.T.reshape((1, -1)); id(c) == aid
Out[11]: False
```
因此，后者的操作将显著慢于前者。

数组的flatten和ravel方法都会将其重塑为一维向量（一个扁平化的数组）。然而，flatten方法总是返回一个副本，而ravel方法仅在必要时返回副本（因此它更快，尤其是在处理大数组时）。

In [12]: d = a.flatten(); id(d) == aid
Out[12]: False
In [13]: e = a.ravel(); id(e) == aid
Out[13]: True
In [14]: %timeit a.flatten()
1000000 loops, best of 3: 1.65 µs per loop
In [15]: %timeit a.ravel()
1000000 loops, best of 3: 566 ns per loop

广播规则允许我们对形状不同但兼容的数组进行计算。换句话说，我们不总是需要重塑或拼接数组来使它们的形状匹配。以下示例演示了在两个向量之间进行外积的两种方法：第一种方法涉及数组拼接，第二种方法（更快）涉及广播：

In [16]: n = 1000
In [17]: a = np.arange(n)
         ac = a[:, np.newaxis]  # Column vector.
         ar = a[np.newaxis, :]  # Row vector.
In [18]: %timeit np.tile(ac, (1, n)) * np.tile(ar, (n, 1))
10 loops, best of 3: 25 ms per loop
In [19]: %timeit ar * ac
100 loops, best of 3: 4.63 ms per loop

它是如何工作的...

在本节中，我们将看到在使用 NumPy 时，内部发生了什么，以及这些知识如何帮助我们理解本食谱中的技巧。

为什么 NumPy 数组高效？

NumPy 数组基本上由元数据描述（特别是维数、形状和数据类型）和实际数据组成。数据存储在一个同质且连续的内存块中，位于系统内存中的特定地址（随机存取存储器，或者RAM）。这个内存块称为数据缓冲区。与纯 Python 结构（如列表）相比，其中项目分散在系统内存中，这是主要区别。这个方面是使 NumPy 数组如此高效的关键特性。

为什么这么重要呢？以下是主要原因：

在低级语言（如 C）中，可以非常高效地编写数组上的计算（NumPy 的大部分实际上是用 C 编写的）。只要知道内存块的地址和数据类型，就可以简单地进行循环遍历所有项目，例如。在 Python 中使用列表进行这样的操作会有很大的开销。
空间局部性在内存访问模式中导致性能提升，这主要是由于 CPU 缓存。事实上，缓存会将字节以块的形式从 RAM 加载到 CPU 寄存器中。然后，相邻的项目会被非常高效地加载（顺序局部性，或者引用局部性）。
最后，项目在内存中连续存储的事实使得 NumPy 能够利用现代 CPU 的矢量化指令，例如英特尔的SSE和AVX，AMD 的 XOP 等。例如，多个连续的浮点数可以加载到 128、256 或 512 位寄存器中，用作实现为 CPU 指令的矢量化算术计算。

注意

另外，NumPy 可以通过ATLAS或英特尔数学核心库（MKL）与高度优化的线性代数库（如BLAS和LAPACK）链接。一些特定的矩阵计算也可以是多线程的，利用现代多核处理器的强大性能。

总之，将数据存储在连续的内存块中确保了现代 CPU 架构在内存访问模式、CPU 缓存和矢量化指令方面的最佳利用。

原地操作和隐式复制操作之间有什么区别？

让我们解释第 2 步中的示例。诸如 a *= 2 这样的表达式对应于原地操作，其中数组的所有值都乘以了二。相比之下，a = a*2 意味着创建了一个包含 a*2 值的新数组，并且变量 a 现在指向这个新数组。旧数组变得没有引用，并将被垃圾回收器删除。与第二种情况相反，第一种情况中不会发生内存分配。

更一般地，诸如 a[i:j] 这样的表达式是数组的视图；它们指向包含数据的内存缓冲区。使用原地操作修改它们会改变原始数组。因此，a[:] = a*2 是一个原地操作，不同于 a = a*2。

了解 NumPy 的这种微妙之处可以帮助你修复一些错误（因为一个数组由于对视图的操作而被隐式和无意中修改），并通过减少不必要的复制次数来优化代码的速度和内存消耗。

为什么有些数组不能在没有复制的情况下重塑？

我们在这里解释第 3 步的示例，其中一个转置的 2D 矩阵不能在没有复制的情况下被展平。一个 2D 矩阵包含由两个数字（行和列）索引的项目，但它在内部存储为一个 1D 连续的内存块，可以用一个数字访问。有多种将矩阵项目存储在 1D 内存块中的方法：我们可以先放第一行的元素，然后是第二行，依此类推，或者先放第一列的元素，然后是第二列，依此类推。第一种方法称为行主序，而后者称为列主序。在这两种方法之间的选择只是内部约定的问题：NumPy 使用行主序，类似于 C，但不同于 FORTRAN。

为什么有些数组不能在没有复制的情况下重塑？

内部数组布局：行主序和列主序

更一般地说，NumPy 使用strides的概念来在多维索引和底层（1D）元素序列的内存位置之间进行转换。array[i1, i2]与内部数据的相关字节地址之间的具体映射如下：

offset = array.strides[0] * i1 + array.strides[1] * i2

在重塑数组时，NumPy 通过修改strides属性来尽可能避免复制。例如，当转置矩阵时，步幅的顺序被颠倒，但底层数据保持不变。然而，简单通过修改步幅来展平一个转置数组是不可能的（试试看！），因此需要进行复制（感谢 Chris Beaumont 澄清了本段早期版本）。

内部数组布局还可以解释一些非常相似的 NumPy 操作之间的意外性能差异。作为一个小练习，你能解释以下基准测试吗？

In [20]: a = np.random.rand(5000, 5000)
In [21]: %timeit a[0,:].sum()
100000 loops, best of 3: 17.9 µs per loop
In [22]: %timeit a[:,0].sum()
10000 loops, best of 3: 60.6 µs per loop

NumPy 广播规则是什么？

广播规则描述了具有不同维度和/或形状的数组如何用于计算。一般规则是当两个维度相等或其中一个为 1 时，它们是兼容的。NumPy 使用此规则逐个元素地比较两个数组的形状，从尾部维度开始，逐步向前工作。最小的维度在内部被拉伸以匹配另一个维度，但此操作不涉及任何内存复制。

还有更多...

以下是一些参考资料：

广播规则和示例，可在docs.scipy.org/doc/numpy/user/basics.broadcasting.html找到。
NumPy 中的数组接口，位于docs.scipy.org/doc/numpy/reference/arrays.interface.html
引用自en.wikipedia.org/wiki/Locality_of_reference的参考资料
SciPy 讲座笔记中的 NumPy 内部，可在scipy-lectures.github.io/advanced/advanced_numpy/找到
Nicolas Rougier 编写的 100 个 NumPy 练习，可在www.loria.fr/~rougier/te…找到

另请参阅

使用 NumPy 的步幅技巧示例

使用 NumPy 的步幅技巧

在这个示例中，我们将深入研究 NumPy 数组的内部，通过将行优先和列优先顺序的概念推广到多维数组。一般概念是步幅，描述多维数组的项在一维数据缓冲区内的组织方式。步幅大多是一个实现细节，但在特定情况下也可以用于优化一些算法。

准备工作

我们假设已经导入了 NumPy，并且id函数已经定义（参见前一个示例，了解 NumPy 的内部以避免不必要的数组复制）。

如何做...

步幅是描述每个维度在连续内存块中的字节步长的整数。
```
In [3]: x = np.zeros(10); x.strides
Out[3]: (8,)
```
这个向量x包含双精度浮点数（float64，8 字节）；从一个项目到下一个项目需要向前移动 8 字节。
现在，让我们看一下一个二维数组的步幅：
```
In [4]: y = np.zeros((10, 10)); y.strides
Out[4]: (80, 8)
```
在第一维（垂直）中，从一个项目到下一个项目需要向前移动 80 字节（10 个 float64 项目），因为项目在内部以行优先顺序存储。在第二维（水平）中，从一个项目到下一个项目需要向前移动 8 字节。
让我们展示如何使用步幅重新审视前一个示例中的广播规则：
```
In [5]: n = 1000; a = np.arange(n)
```
我们将创建一个新数组b，指向与a相同的内存块，但形状和步幅不同。这个新数组看起来像是a的垂直平铺版本。我们使用 NumPy 中的一个特殊函数来改变数组的步幅：
```
In [6]: b = np.lib.stride_tricks.as_strided(a, (n, n), 
                                           (0, 4))
Out[7]: array([[ 0, 1, 2, ..., 997, 998, 999],
             ..., 
             [ 0, 1, 2, ..., 997, 998, 999]])
In [8]: b.size, b.shape, b.nbytes
Out[8]: (1000000, (1000, 1000), 4000000)
```
NumPy 认为这个数组包含一百万个不同的元素，而实际上数据缓冲区只包含与a相同的 1000 个元素。

现在我们可以使用与广播规则相同的原则执行高效的外积：

In [9]: %timeit b * b.T
100 loops, best of 3: 5.03 ms per loop
In [10]: %timeit np.tile(a, (n, 1)) \
               * np.tile(a[:, np.newaxis], (1, n))
10 loops, best of 3: 28 ms per loop

工作原理...

每个数组都有一定数量的维度、形状、数据类型和步幅。步幅描述了多维数组中的项是如何在数据缓冲区中组织的。有许多不同的方案来安排多维数组的项在一维块中。NumPy 实现了一个跨步索引方案，其中任何元素的位置是维度的线性组合，系数是步幅。换句话说，步幅描述了在任何维度中，我们需要跳过多少字节在数据缓冲区中从一个项到下一个项。

多维数组中任何元素的位置由其索引的线性组合给出，如下所示：

工作原理...

人为改变步幅允许我们使一些数组操作比标准方法更有效，后者可能涉及数组复制。在内部，这就是 NumPy 中广播的工作原理。

as_strided方法接受一个数组、一个形状和步幅作为参数。它创建一个新数组，但使用与原始数组相同的数据缓冲区。唯一改变的是元数据。这个技巧让我们像平常一样操作 NumPy 数组，只是它们可能比 NumPy 认为的占用更少的内存。在这里，使用步幅中的 0 意味着任何数组项都可以由许多多维索引来寻址，从而节省内存。

提示

对于跨步数组要小心！as_strided函数不会检查你是否保持在内存块边界内。这意味着你需要手动处理边缘效应；否则，你的数组中可能会出现垃圾值。

我们将在下一个食谱中看到跨步技巧的更有用的应用。

参见

使用跨步技巧实现高效的滚动平均算法食谱

使用跨步技巧实现高效的滚动平均算法

当在数组上进行局部计算时，跨步技巧可以很有用，当给定位置的计算值取决于相邻值时。示例包括动力系统、数字滤波器和细胞自动机。

在这个食谱中，我们将使用 NumPy 跨步技巧实现一个高效的滚动平均算法（一种特殊类型的基于卷积的线性滤波器）。1D 向量的滚动平均在每个位置包含原始向量中该位置周围元素的平均值。粗略地说，这个过程滤除信号的嘈杂成分，以保留只有较慢的成分。

准备工作

确保重用了解 NumPy 内部以避免不必要的数组复制食谱中的id()函数。这个函数返回 NumPy 数组的内部数据缓冲区的内存地址。

如何做...

思路是从一个 1D 向量开始，并创建一个虚拟的 2D 数组，其中每一行都是前一行的移位版本。使用跨步技巧时，这个过程非常高效，因为它不涉及任何复制。

让我们生成一个一维向量：

In [1]: import numpy as np
        from numpy.lib.stride_tricks import as_strided
In [2]: n = 5; k = 2
In [3]: a = np.linspace(1, n, n); aid = id(a)

让我们改变a的步幅以添加移位行：

In [4]: as_strided(a, (k, n), (8, 8))
Out[4]: array([[ 1e+00,  2e+00,  3e+00,  4e+00,  5e+00],
               [ 2e+00,  3e+00,  4e+00,  5e+00, -1e-23]])

最后一个值表示越界问题：步幅技巧可能存在危险，因为内存访问没有经过检查。在这里，我们应该通过限制数组的形状来考虑边缘效应。

In [5]: as_strided(a, (k, n-k+1), (8, 8))
Out[5]: array([[ 1.,  2.,  3.,  4.],
               [ 2.,  3.,  4.,  5.]])

现在，让我们实现滚动平均值的计算。第一个版本（标准方法）涉及显式数组复制，而第二个版本使用步幅技巧：

In [6]: def shift1(x, k):
            return np.vstack([x[i:n-k+i+1] for i in 
                              range(k)])
In [7]: def shift2(x, k):
            return as_strided(x, (k, n-k+1), 
                              (x.itemsize,)*2)

这两个函数返回相同的结果，只是第二个函数返回的数组指向原始数据缓冲区：

In [8]: b = shift1(a, k); b, id(b) == aid
Out[8]: (array([[ 1.,  2.,  3.,  4.],
                [ 2.,  3.,  4.,  5.]]), False) 
In [9]: c = shift2(a, k); c, id(c) == aid
Out[9]: (array([[ 1.,  2.,  3.,  4.],
                [ 2.,  3.,  4.,  5.]]), True)

让我们生成一个信号：

In [10]: n, k = 100, 10
         t = np.linspace(0., 1., n)
         x = t + .1 * np.random.randn(n)

我们通过创建信号的移位版本，并沿着垂直维度进行平均，来计算信号的滚动平均值。结果显示在下一个图中：
```
In [11]: y = shift2(x, k)
         x_avg = y.mean(axis=0)
```
一个信号及其滚动平均值

让我们评估第一种方法所花费的时间：

In [15]: %timeit shift1(x, k)
10000 loops, best of 3: 163 µs per loop
In [16]: %%timeit y = shift1(x, k)
         z = y.mean(axis=0)
10000 loops, best of 3: 63.8 µs per loop

以及第二种方法：

In [17]: %timeit shift2(x, k)
10000 loops, best of 3: 23.3 µs per loop
In [18]: %%timeit y = shift2(x, k)
         z = y.mean(axis=0)
10000 loops, best of 3: 35.8 µs per loop

在第一个版本中，大部分时间都花在数组复制上，而在步幅技巧版本中，大部分时间都花在平均值的计算上。

参见

使用 NumPy 的步幅技巧方法

在 NumPy 中进行高效的数组选择

NumPy 提供了几种选择数组切片的方法。数组视图指的是数组的原始数据缓冲区，但具有不同的偏移、形状和步幅。它们只允许步进选择（即，具有线性间隔的索引）。NumPy 还提供了特定函数来沿一个轴进行任意选择。最后，花式索引是最通用的选择方法，但也是最慢的，正如我们在这个方法中所看到的。在可能的情况下，应选择更快的替代方法。

准备工作

我们假设 NumPy 已经被导入，并且id函数已经被定义（参见了解 NumPy 的内部以避免不必要的数组复制方法）。

如何实现...

让我们创建一个具有大量行的数组。我们将沿着第一个维度选择这个数组的切片：
```
In [3]: n, d = 100000, 100
In [4]: a = np.random.random_sample((n, d)); aid = id(a)
```

让我们选择每 10 行中的一行，使用两种不同的方法（数组视图和花式索引）：

In [5]: b1 = a[::10]
        b2 = a[np.arange(0, n, 10)]
In [6]: np.array_equal(b1, b2)
Out[6]: True

视图指向原始数据缓冲区，而花式索引产生一个副本：
```
In [7]: id(b1) == aid, id(b2) == aid
Out[7]: (True, False)
```

让我们比较两种方法的性能：

In [8]: %timeit a[::10]
100000 loops, best of 3: 2.03 µs per loop
In [9]: %timeit a[np.arange(0, n, 10)]
10 loops, best of 3: 46.3 ms per loop

花式索引比涉及复制大数组的速度慢几个数量级。

当需要沿着一个维度进行非步进选择时，数组视图不是一个选择。然而，在这种情况下，仍然存在替代花式索引的方法。给定一个索引列表，NumPy 的take()函数可以沿着一个轴进行选择：
```
In [10]: i = np.arange(0, n, 10)
In [11]: b1 = a[i]
         b2 = np.take(a, i, axis=0)
In [12]: np.array_equal(b1, b2)
Out[12]: True
```
第二种方法更快：
```
In [13]: %timeit a[i]
10 loops, best of 3: 50.2 ms per loop
In [14]: %timeit np.take(a, i, axis=0)
100 loops, best of 3: 11.1 ms per loop
```
提示

最近的 NumPy 版本中改进了花式索引的性能；这个技巧在旧版本的 NumPy 中特别有用。

当沿着一个轴选择的索引由布尔掩码向量指定时，compress()函数是花式索引的一种替代方法：

In [15]: i = np.random.random_sample(n) < .5
In [16]: b1 = a[i]
         b2 = np.compress(i, a, axis=0) 
In [17]: np.array_equal(b1, b2)
Out[17]: True
In [18]: %timeit a[i]
1 loops, best of 3: 286 ms per loop
In [19]: %timeit np.compress(i, a, axis=0)
10 loops, best of 3: 41.3 ms per loop

第二种方法也比花式索引更快。

工作原理是怎样的...

花式索引是进行完全任意选择数组的最通用方式。然而，通常存在更具体且更快的方法，在可能的情况下应该优先选择这些方法。

当需要进行步进选择时，应使用数组视图，但我们需要小心，因为视图是引用原始数据缓冲区的。

还有更多内容…

这里有一些参考资料：

完整的 NumPy 例程列表可以在 NumPy 参考指南中找到，网址是docs.scipy.org/doc/numpy/reference/
索引例程的列表可以在docs.scipy.org/doc/numpy/reference/routines.indexing.html查看

使用内存映射处理超大 NumPy 数组

有时，我们需要处理无法完全放入系统内存的超大 NumPy 数组。一种常见的解决方案是使用内存映射并实现核心外计算。数组被存储在硬盘上的文件中，我们创建一个内存映射对象来访问该文件，该对象可以像普通的 NumPy 数组一样使用。访问数组的一部分会自动从硬盘中获取相应的数据。因此，我们只消耗我们使用的部分。

如何操作…

让我们创建一个内存映射数组：

In [1]: import numpy as np
In [2]: nrows, ncols = 1000000, 100
In [3]: f = np.memmap('memmapped.dat', dtype=np.float32, 
                      mode='w+', shape=(nrows, ncols))

让我们一次性喂入随机值，一次处理一列，因为我们的系统内存有限！
```
In [4]: for i in range(ncols):
            f[:,i] = np.random.rand(nrows)
```
我们保存数组的最后一列：
```
In [5]: x = f[:,-1]
```
现在，我们通过删除对象将内存中的更改刷新到磁盘：
```
In [6]: del f
```

从磁盘读取内存映射数组使用相同的memmap函数。数据类型和形状需要重新指定，因为这些信息不会保存在文件中：

In [7]: f = np.memmap('memmapped.dat', dtype=np.float32,
                      shape=(nrows, ncols))
In [8]: np.array_equal(f[:,-1], x)
Out[8]: True
In [9]: del f

提示

这种方法并不是最适合长期存储数据和数据共享的。接下来的章节将展示一种基于 HDF5 文件格式的更好的方法。

工作原理…

内存映射使你几乎可以像使用常规数组一样处理超大数组。接受 NumPy 数组作为输入的 Python 代码也可以接受memmap数组。然而，我们需要确保高效使用该数组。也就是说，数组绝不能一次性全部加载（否则会浪费系统内存，并失去该技术的优势）。

当你有一个包含已知数据类型和形状的原始二进制格式的大文件时，内存映射也很有用。在这种情况下，另一种解决方案是使用 NumPy 的fromfile()函数，并通过 Python 的原生open()函数创建文件句柄。使用f.seek()可以将光标定位到任何位置，并将指定数量的字节加载到 NumPy 数组中。

还有更多内容…

处理巨大的 NumPy 矩阵的另一种方法是通过 SciPy 的稀疏子包使用稀疏矩阵。当矩阵大部分包含零时，这种方法非常适用，这在偏微分方程的模拟、图算法或特定的机器学习应用中常常出现。将矩阵表示为稠密结构可能浪费内存，而稀疏矩阵提供了更高效的压缩表示。

在 SciPy 中使用稀疏矩阵并不简单，因为存在多种实现方法。每种实现最适合特定类型的应用。以下是一些参考资料：

关于稀疏矩阵的 SciPy 讲义，见scipy-lectures.github.io/advanced/scipy_sparse/index.html
稀疏矩阵的参考文档，请见docs.scipy.org/doc/scipy/reference/sparse.html
memmap 文档，请见docs.scipy.org/doc/numpy/reference/generated/numpy.memmap.html

另请参见

使用 HDF5 和 PyTables 操作大型数组食谱
使用 HDF5 和 PyTables 操作大型异构表格食谱

使用 HDF5 和 PyTables 操作大型数组

NumPy 数组可以通过 NumPy 内建函数如np.savetxt、np.save或np.savez持久化保存到磁盘，并使用类似的函数加载到内存中。当数组包含的点数少于几百万时，这些方法表现最好。对于更大的数组，这些方法面临两个主要问题：它们变得过于缓慢，并且需要将数组完全加载到内存中。包含数十亿个点的数组可能太大，无法适应系统内存，因此需要替代方法。

这些替代方法依赖于内存映射：数组存储在硬盘上，只有当 CPU 需要时，数组的部分数据才会被有选择性地加载到内存中。这种技术在节省内存方面非常高效，但由于硬盘访问会带来轻微的开销。缓存机制和优化可以缓解这个问题。

上一个食谱展示了一种使用 NumPy 的基本内存映射技术。在这个食谱中，我们将使用一个名为PyTables的包，它专门用于处理非常大的数据集。它通过一个广泛使用的开放文件格式规范层次数据格式（HDF5）实现了快速的内存映射技术。一个 HDF5 文件包含一个或多个数据集（数组或异构表格），这些数据集按类 Unix 的 POSIX 层次结构组织。数据集的任何部分都可以高效且轻松地访问，而无需不必要地浪费系统内存。

正如我们在下面的食谱中将看到的，PyTables 的一个替代方案是h5py。在某些情况下，它比 PyTables 更加轻量且更适应需求。

在这个食谱中，我们将学习如何使用 HDF5 和 PyTables 操作大数组。下一个食谱将讲解类似 pandas 的异构表格。

准备工作

本食谱和下一个食谱需要 PyTables 3.0+ 版本。使用 Anaconda，可以通过 conda install tables 安装 PyTables。你还可以在 pytables.github.io/usersguide/installation.html 找到二进制安装程序。Windows 用户可以在 www.lfd.uci.edu/~gohlke/pyt… 上找到安装程序。

提示

在 3.0 版本之前，PyTables 使用驼峰命名法约定来命名属性和方法。最新版本使用更标准的 Python 约定，使用下划线。所以，例如，tb.open_file 在 3.0 版本之前是 tb.openFile。

如何操作…

首先，我们需要导入 NumPy 和 PyTables（该包的名称是 tables）：
```
In [1]: import numpy as np
        import tables as tb
```

让我们创建一个新的空 HDF5 文件：

In [2]: f = tb.open_file('myfile.h5', 'w')

我们创建了一个名为 experiment1 的新顶级组：

In [3]: f.create_group('/', 'experiment1')
Out[3]: /experiment1 (Group) u''
  children := []

让我们还向这个组添加一些元数据：

In [4]: f.set_node_attr('/experiment1', 'date', '2014-09-01')

在这个组中，我们创建了一个 1000*1000 的数组，命名为 array1：

In [5]: x = np.random.rand(1000, 1000)
        f.create_array('/experiment1', 'array1', x)
Out[5]: /experiment1/array1 (Array(1000L, 1000L))

最后，我们需要关闭文件以提交磁盘上的更改：
```
In [6]: f.close()
```
现在，让我们打开这个文件。我们也可以在另一个 Python 会话中完成这项操作，因为数组已经保存在 HDF5 文件中。
```
In [7]: f = tb.open_file('myfile.h5', 'r')
```

我们可以通过提供组路径和属性名称来检索属性：

In [8]: f.get_node_attr('/experiment1', 'date')
Out[8]: '2014-09-01'

我们可以使用属性访问文件中的任何项目，通过用点替代路径中的斜杠，并以 root 开头（对应路径 /）。IPython 的自动补全功能在互动式探索文件时特别有用。
```
In [9]: y = f.root.experiment1.array1  
        # /experiment1/array1
        type(y)
Out[9]: tables.array.Array
```
这个数组可以像 NumPy 数组一样使用，但一个重要的区别是它存储在磁盘上，而不是系统内存中。在这个数组上执行计算时，会自动将请求的数组部分加载到内存中，因此只访问数组的视图会更高效。
```
In [10]: np.array_equal(x[0,:], y[0,:])
Out[10]: True
```

也可以通过绝对路径获取节点，这在路径只在运行时已知时非常有用：

In [11]: f.get_node('/experiment1/array1')
Out[11]: /experiment1/array1 (Array(1000, 1000))

本食谱完成了，接下来让我们做一些清理工作：

In [12]: f.close()
In [13]: import os
         os.remove('myfile.h5')

它是如何工作的…

在这个食谱中，我们将一个单独的数组存储到文件中，但当需要在一个文件中保存多个数组时，HDF5 特别有用。HDF5 通常用于大型项目，特别是当大数组必须在层次结构中组织时。例如，它在 NASA、NOAA 和许多其他科学机构中被广泛使用（请参见 www.hdfgroup.org/users.html）。研究人员可以在多个设备、多个试验和多个实验中存储记录的数据。

在 HDF5 中，数据是以树的形式组织的。树的节点可以是组（类似于文件系统中的文件夹）或数据集（类似于文件）。一个组可以包含子组和数据集，而数据集只包含数据。组和数据集都可以包含具有基本数据类型（整数或浮点数、字符串等）的属性（元数据）。HDF5 还支持内部和外部链接；某个路径可以引用同一文件中的另一个组或数据集，或者引用另一个文件中的组或数据集。如果需要为同一实验或项目使用不同的文件，这一特性可能会很有用。

能够在不将整个数组及文件加载到内存中的情况下访问单个数组的一个块是非常方便的。此外，加载后的数组可以使用标准的 NumPy 切片语法进行多态访问。接受 NumPy 数组作为参数的代码原则上也可以接受 PyTables 数组对象作为参数。

还有更多...

在本食谱中，我们创建了一个 PyTables Array 对象来存储我们的 NumPy 数组。其他类似的对象类型包括 CArrays，它们将大数组分块存储并支持无损压缩。此外，EArray 对象在最多一个维度上是可扩展的，这在创建文件中的数组时，如果数组的维度未知时，非常有用。一个常见的使用场景是在在线实验中记录数据。

另一个主要的 HDF5 对象类型是 Table，它在一个二维表中存储异构数据类型的表格数据。在 PyTables 中，Table 与 Array 的关系就像 pandas 的 DataFrame 与 NumPy 的 ndarray 之间的关系。我们将在下一个食谱中看到它们。

HDF5 文件的一个有趣特点是它们不依赖于 PyTables。由于 HDF5 是一个开放的格式规范，因此几乎所有语言（C、FORTRAN、MATLAB 等）都有相关的库，因此可以轻松地在这些语言中打开 HDF5 文件。

在 HDF5 中，数据集可以存储在连续的内存块中，或者分块存储。块（Chunks）是原子对象，HDF5/PyTables 只能读取和写入整个块。块在内部通过一种树形数据结构组织，称为B 树。当我们创建一个新的数组或表时，可以指定块形状。这是一个内部细节，但它在写入和读取数据集的部分时可能会极大影响性能。

最优的块形状取决于我们计划如何访问数据。小块的数量较多（由于管理大量块导致较大的开销）与大块较少（磁盘 I/O 效率低）之间存在权衡。一般来说，建议块的大小小于 1 MB。块缓存也是一个重要参数，可能会影响性能。

相关地，在我们创建 EArray 或 Table 对象时，应该指定一个可选参数，预期的行数，以便优化文件的内部结构。如果你打算在大型 HDF5 文件（超过 100 MB）上进行任何稍复杂的操作，PyTables 用户指南中的优化部分（参见以下参考链接）是必须阅读的。

最后，我们应该提到另一个名为 h5py 的 Python HDF5 库。这个轻量级库提供了一个简单的接口来操作 HDF5 文件，侧重于数组而非表格。它为从 NumPy 访问 HDF5 数组提供了非常自然的方式，如果你不需要 PyTables 中类似数据库的功能，它可能已经足够。有关 h5py 的更多信息，请参见 www.h5py.org。

下面是一些参考资料：

HDF5 分块，详见 www.hdfgroup.org/HDF5/doc/Ad…
PyTables 优化指南，详见 pytables.github.io/usersguide/optimization.html
PyTables 和 h5py 的区别，从 h5py 的角度来看，详见 github.com/h5py/h5py/wiki/FAQ#whats-the-difference-between-h5py-and-pytables
PyTables 和 h5py 的区别，从 PyTables 的角度来看，详见 www.pytables.org/moin/FAQ#Ho…

另请参见

处理巨大 NumPy 数组的内存映射 配方
使用 HDF5 和 PyTables 操作大型异构表格 配方
在第二章，交互计算的最佳实践 中的 进行可重复交互计算实验的十个提示 配方

使用 HDF5 和 PyTables 操作大型异构表格

PyTables 可以将同质数据块存储为类似 NumPy 数组的形式在 HDF5 文件中。它也可以存储异构表格，如我们将在这个配方中看到的那样。

准备工作

你需要 PyTables 来实现这个配方（有关安装说明，请参见前面的配方）。

如何实现...

让我们导入 NumPy 和 PyTables：

In [1]: import numpy as np
        import tables as tb

让我们创建一个新的 HDF5 文件：

In [2]: f = tb.open_file('myfile.h5', 'w')

我们将创建一个 HDF5 表格，包含两列：城市名称（最多 64 个字符的字符串）和其人口（32 位整数）。我们可以通过创建一个复杂的数据类型，并使用 NumPy 来指定这些列：
```
In [3]: dtype = np.dtype([('city','S64'),
                          ('population', 'i4')])
```

现在，我们在 /table1 中创建表格：

In [4]: table = f.create_table('/', 'table1', dtype)

让我们添加几行数据：

In [5]: table.append([('Brussels', 1138854),
                      ('London',   8308369),
                      ('Paris',    2243833)])

添加行数据后，我们需要刷新表格以提交更改到磁盘：
```
In [6]: table.flush()
```

有多种方式可以访问表格中的数据。最简单但效率不高的方式是将整个表格加载到内存中，这将返回一个 NumPy 数组：

In [7]: table[:]
Out[7]: array([('Brussels', 1138854),
               ('London', 8308369), 
               ('Paris', 2243833)], 
               dtype=[('city', 'S64'), 
                      ('population', '<i4')])

也可以加载特定的列（所有行）：

In [8]: table.col('city')
Out[8]: array(['Brussels', 'London', 'Paris'],
              dtype='|S64')

在处理大量行时，我们可以在表格中执行类似 SQL 的查询，加载所有满足特定条件的行：

In [9]: [row['city'] for row in 
         table.where('population>2e6')]
Out[9]: ['London', 'Paris']

最后，如果已知索引，我们可以访问特定的行：
```
In [10]: table[1]
Out[10]: ('London', 8308369)
```

工作原理...

表格可以像本教程中那样从头开始创建，也可以从现有的 NumPy 数组或 pandas DataFrame创建。在第一种情况下，可以使用 NumPy 数据类型（如这里所示）、字典或继承自IsDescription的类来描述列。在第二种情况下，表格描述将根据给定的数组或表格自动推断。

可以通过table.append()高效地在表格末尾添加行。要添加一行，首先使用row = table.row获取一个新行实例，像操作字典一样设置行的字段，然后调用row.append()将新行添加到表格的末尾。在一组写操作后调用flush()可以确保这些更改同步到磁盘。PyTables 使用复杂的缓存机制来确保在表格中读写数据时获得最大性能，因此新行不会立即写入磁盘。

PyTables 支持高效的 SQL-like 查询，称为内核查询。包含查询表达式的字符串会在所有行上编译并评估。较低效的方法是使用table.iterrows()遍历所有行，并对行的字段使用if语句。

还有更多内容...

这里有一些参考资料：

内核查询，详见pytables.github.io/usersguide/condition_syntax.html。
PyTables 和 HDF5 的替代方案可能来自 Blaze 项目，该项目在撰写时仍处于早期开发阶段。有关 Blaze 的更多信息，请参考blaze.pydata.org。

另见

使用 HDF5 和 PyTables 操作大数组的教程

第五章：高性能计算

在本章中，我们将涵盖以下主题：

使用 Numba 和即时编译加速纯 Python 代码
使用 Numexpr 加速数组计算
使用 ctypes 将 C 库包装为 Python
使用 Cython 加速 Python 代码
通过编写更少的 Python 代码和更多的 C 代码来优化 Cython 代码
使用 Cython 和 OpenMP 释放 GIL，利用多核处理器
使用 CUDA 为 NVIDIA 图形卡（GPU）编写大规模并行代码
使用 OpenCL 为异构平台编写大规模并行代码
使用 IPython 将 Python 代码分配到多个核心
在 IPython 中与异步并行任务进行交互
在 IPython 中使用 MPI 并行化代码
在 notebook 中尝试 Julia 语言

介绍

上一章介绍了代码优化的技术。有时，这些方法不足以满足需求，我们需要诉诸于更高级的高性能计算技术。

在本章中，我们将看到三种广泛的方法类别，它们并不互相排斥：

即时编译（JIT）Python 代码
从 Python 转向低级语言，如 C
使用并行计算将任务分配到多个计算单元

通过即时编译，Python 代码会动态编译成低级语言。编译发生在运行时，而不是执行之前。由于代码是编译的而非解释的，因此运行速度更快。JIT 编译是一种流行的技术，因为它能够同时实现快速和高级语言，而这两个特点在过去通常是互相排斥的。

JIT 编译技术已经在如 Numba、Numexpr、Blaze 等包中实现。在本章中，我们将介绍前两个包。Blaze 是一个有前景的项目，但在写本书时，它仍处于起步阶段。

我们还将介绍一种新的高级语言，Julia，它使用 JIT 编译实现高性能。这种语言可以在 IPython notebook 中有效使用，得益于 IJulia 包。

注意

PyPy (pypy.org)，Psyco 的继任者，是另一个相关的项目。这个 Python 的替代实现（参考实现是 CPython）集成了 JIT 编译器。因此，它通常比 CPython 快得多。然而，在本书写作时，PyPy 还不完全支持 NumPy。此外，PyPy 和 SciPy 往往形成各自独立的社区。

求助于如 C 这样的低级语言是另一种有趣的方法。流行的库包括ctypes、SWIG 或 Cython。使用 ctypes 需要编写 C 代码并能够访问 C 编译器，或者使用编译过的 C 库。相比之下，Cython 允许我们在 Python 的超集语言中编写代码，然后将其转译为 C，并带来不同的性能结果。不幸的是，编写高效的 Cython 代码并不总是容易的。在本章中，我们将介绍 ctypes 和 Cython，并展示如何在复杂示例中实现显著的加速。

最后，我们将介绍两类并行计算技术：使用 IPython 利用多个 CPU 核心和使用大规模并行架构，如图形处理单元（GPU）。

这里有一些参考资料：

Travis Oliphant 关于 PyPy 和 NumPy 的博客文章，内容可以在technicaldiscovery.blogspot.com/2011/10/thoughts-on-porting-numpy-to-pypy.html找到
与 C 语言接口的 Python 教程，详细内容请参考scipy-lectures.github.io/advanced/interfacing_with_c/interfacing_with_c.html

CPython 和并发编程

Python 有时因为对多核处理器的本地支持较差而受到批评。让我们来解释一下原因。

Python 语言的主流实现是用 C 编写的CPython。CPython 集成了一种机制，称为全局解释器锁（GIL）。正如在wiki.python.org/moin/GlobalInterpreterLock中所提到的：

GIL 通过防止多个本地线程同时执行 Python 字节码来促进内存管理。

换句话说，通过禁用一个 Python 进程中的并发线程，GIL 大大简化了内存管理系统。因此，内存管理在 CPython 中并不是线程安全的。

一个重要的影响是，使用 CPython 时，一个纯 Python 程序无法轻松地在多个核心上并行执行。这是一个重要问题，因为现代处理器的核心数量越来越多。

我们有什么可能的解决方案，以便利用多核处理器？

移除 CPython 中的 GIL。这个解决方案曾经被尝试过，但从未被纳入 CPython。这会在 CPython 的实现中带来太多复杂性，并且会降低单线程程序的性能。
使用多个进程而非多个线程。这是一个流行的解决方案；可以使用原生的multiprocessing模块，或者使用 IPython。在本章中，我们将介绍后者。
用 Cython 重写代码的特定部分，并用 C 语言变量替换所有 Python 变量。这样可以暂时移除 GIL，使得在循环中能够使用多核处理器。我们将在释放 GIL 以利用 Cython 和 OpenMP 实现多核处理器的配方中讨论这一解决方案。
用一种对多核处理器有更好支持的语言实现代码的特定部分，并从你的 Python 程序中调用它。
让你的代码使用可以从多核处理器中受益的 NumPy 函数，如numpy.dot()。NumPy 需要用 BLAS/LAPACK/ATLAS/MKL 进行编译。

关于 GIL 的必读参考资料可以在www.dabeaz.com/GIL/找到。

与编译器相关的安装说明

在本节中，我们将提供一些使用编译器与 Python 配合的说明。使用场景包括使用 ctypes、使用 Cython 以及构建 Python 的 C 扩展。

Linux

在 Linux 上，你可以安装GCC（GNU 编译器集合）。在 Ubuntu 或 Debian 上，你可以通过命令sudo apt-get install build-essential安装 GCC。

Mac OS X

在 Mac OS X 上，你可以安装 Apple XCode。从 XCode 4.3 开始，你必须通过 XCode 菜单中的Preferences | Downloads | Command Line Tools手动安装命令行工具。

Windows

在 Windows 上，使用编译器与 Python 配合使用是著名的麻烦。通常很难在线找到所有必要的说明。我们在这里详细说明这些说明（你也可以在本书的 GitHub 仓库中找到它们）：

说明根据你使用的是 32 位还是 64 位版本的 Python，以及你使用的是 Python 2.x 还是 Python 3.x 有所不同。要快速了解这些信息，你可以在 Python 终端中键入以下命令：

import sys
print(sys.version)
print(64 if sys.maxsize > 2**32 else 32)

Python 32 位

首先，你需要安装一个 C 编译器。对于 Python 32 位，你可以从www.mingw.org下载并安装 MinGW，它是 GCC 的开源发行版。
根据你所使用的distutils库的版本，你可能需要手动修复其源代码中的一个 bug。用文本编辑器打开C:\Python27\Lib\distutils\cygwinccompiler.py（或者根据你的具体配置路径类似的文件），并将所有-mno-cygwin的出现位置替换为空字符串。
打开或创建一个名为distutils.cfg的文本文件，路径为C:\Python27\Lib\distutils\，并添加以下行：
```
[build]
compiler = mingw32
```

Python 64 位

对于 Python 2.x，你需要 Visual Studio 2008 Express。对于 Python 3.x，你需要 Visual Studio 2010 Express。
你还需要微软 Windows SDK（根据你的 Python 版本选择 2008 或 2010 版本）：
- Python 2.x：可在www.microsoft.com/en-us/download/details.aspx?id=3138下载 Windows 7 和.NET Framework 3.5 的微软 Windows SDK。
- Python 3.x：适用于 Windows 7 和.NET Framework 4 的 Microsoft Windows SDK，下载地址：www.microsoft.com/en-us/download/details.aspx?id=8279
确保包含cl.exe的文件夹路径在系统的PATH环境变量中。该路径应类似于C:\Program Files (x86)\Microsoft Visual Studio 9.0\VC\bin\amd64（使用适用于 Windows 7 和.NET Framework 3.5 的 Microsoft Windows SDK 中的 Visual Studio 2008 C 编译器）。
每次你想要在 Python 中使用编译器时（例如，在输入ipython notebook之前），需要在 Windows 的命令行终端中执行几个命令：
```
call "C:\Program Files\Microsoft SDKs\Windows\v7.0\Bin\SetEnv.Cmd" /x64 /release
set DISTUTILS_USE_SDK=1
```

DLL 地狱

使用编译包时，特别是那些在 Chris Gohlke 的网页上获取的包（www.lfd.uci.edu/~gohlke/pythonlibs/），可能会遇到一些晦涩的 DLL 相关错误。要解决这些问题，你可以使用 Dependency Walker 打开这些无效的 DLL，Dependency Walker 的下载地址为www.dependencywalker.com。该程序可以告诉你缺少了哪个 DLL 文件，你可以在电脑上搜索它并将其位置添加到PATH环境变量中。

参考资料

以下是一些参考资料：

在 Windows 上安装 Cython，参考wiki.cython.org/InstallingOnWindows
在 Windows 64 位上使用 Cython，参考github.com/cython/cython/wiki/64BitCythonExtensionsOnWindows
为 Windows 构建 Python wheel，参考cowboyprogrammer.org/building-python-wheels-for-windows/

使用 Numba 和即时编译加速纯 Python 代码

Numba（numba.pydata.org）是由 Continuum Analytics（www.continuum.io）创建的一个包。截止本文撰写时，Numba 仍然是一个较新的、相对实验性的包，但其技术前景可期。Numba 会自动（即时）将纯 Python 代码翻译成优化后的机器码。实际上，这意味着我们可以在纯 Python 中编写一个非矢量化的函数，使用for循环，并通过使用一个装饰器自动将该函数矢量化。与纯 Python 代码相比，性能的提升可以达到几个数量级，甚至可能超越手动矢量化的 NumPy 代码。

在本节中，我们将展示如何加速生成曼德尔布罗特分形的纯 Python 代码。

准备工作

安装 Numba 最简单的方法是使用 Anaconda 发行版（也是 Continuum Analytics 维护的），然后在终端输入conda install numba。在 Windows 上，另一种选择是从 Chris Gohlke 的网页下载二进制安装程序，网址为www.lfd.uci.edu/~gohlke/pythonlibs/#numba。在这种情况下，有一些依赖项（Numpy-MKL、LLVMPy、llvmmath 和 Meta），它们都可以在同一页面上找到。

如何实现…

让我们导入 NumPy 并定义几个变量：

In [1]: import numpy as np
In [2]: size = 200
        iterations = 100

以下函数使用纯 Python 生成分形。它接受一个空数组 m 作为参数。

In [3]: def mandelbrot_python(m, size, iterations):
            for i in range(size):
                for j in range(size):
                    c = -2 + 3./size*j + 1j*(1.5-3./size*i)
                    z = 0
                    for n in range(iterations):
                        if np.abs(z) <= 10:
                            z = z*z + c
                            m[i, j] = n
                        else:
                            break

让我们运行这个模拟并显示分形图：

In [4]: m = np.zeros((size, size))
        mandelbrot_python(m, size, iterations)
In [5]: import matplotlib.pyplot as plt
        %matplotlib inline
        plt.imshow(np.log(m), cmap=plt.cm.hot)
        plt.xticks([]); plt.yticks([])

如何实现…

曼德尔布罗特分形

现在，我们评估这个函数所用的时间：

In [6]: %%timeit m = np.zeros((size, size))
        mandelbrot_python(m, size, iterations)
1 loops, best of 1: 6.18 s per loop

让我们尝试使用 Numba 加速这个函数。首先，我们导入这个包：
```
In [7]: import numba
        from numba import jit, complex128
```

接下来，我们在函数定义上方添加@jit装饰器。Numba 会尝试自动推断局部变量的类型，但我们也可以显式地指定类型：

In [8]: @jit(locals=dict(c=complex128, z=complex128))
        def mandelbrot_numba(m, size, iterations):
            for i in range(size):
                for j in range(size):
                    c = -2 + 3./size*j + 1j*(1.5-3./size*i)
                    z = 0
                    for n in range(iterations):
                        if np.abs(z) <= 10:
                            z = z*z + c
                            m[i, j] = n
                        else:
                            break

这个函数与纯 Python 版本的功能完全相同。它到底快了多少呢？
```
In [10]: %%timeit m = np.zeros((size, size))
         mandelbrot_numba(m, size, iterations)
1 loops, best of 10: 44.8 ms per loop
```
这里 Numba 版本的速度比纯 Python 版本快超过 100 倍！

它是如何工作的…

Python 字节码通常在运行时由 Python 解释器（例如 CPython）解释执行。相比之下，Numba 函数会在执行前直接解析并转换为机器代码，使用名为LLVM（低级虚拟机）的强大编译器架构。引用官方文档：

Numba 能识别 NumPy 数组作为已类型化的内存区域，因此可以加速使用 NumPy 数组的代码。其他类型不太明确的代码将被转换为 Python C-API 调用，从而有效地去除了“解释器”，但并未去除动态间接性。

Numba 并不是能够编译所有 Python 函数。对于局部变量的类型也有一些微妙的限制。Numba 会尝试自动推断函数变量的类型，但并不总是成功。在这种情况下，我们可以显式地指定类型。

Numba 通常在涉及 NumPy 数组的紧密循环（如本食谱中的示例）时能提供最显著的加速。

注意

Blaze 是 Continuum Analytics 的另一个项目，它是 NumPy 的下一代版本。它将提供比 NumPy 数组更具灵活性的数据结构，并且支持外存计算。与 Numba 一起，Blaze 将形成一个高效的类似编译器的基础设施，用于大数据算法和复杂的数值模拟。我们可以预期 Blaze 在未来会发挥重要作用，因为它应该将 Python 简洁易用的语法与原生代码的性能和并行处理技术（特别是多核处理器和图形处理单元）相结合。其他值得一提的相关项目，虽然稍微比 Blaze 和 Numba 旧，但仍有价值，包括Theano和Numexpr（我们将在下一个食谱中看到它们）。

还有更多内容…

让我们比较 Numba 与使用 NumPy 手动向量化代码的性能，后者是加速纯 Python 代码的标准方法，例如本教程中的代码。实际上，这意味着将 i 和 j 循环中的代码替换为数组计算。在这里，这相对容易，因为这些操作紧跟 单指令多数据 (SIMD) 的范式：

In [1]: import numpy as np
        import matplotlib.pyplot as plt
        %matplotlib inline
In [2]: def initialize(size):
            x, y = np.meshgrid(np.linspace(-2, 1, size),
                               np.linspace(-1.5, 1.5, size))
            c = x + 1j*y
            z = c.copy()
            m = np.zeros((size, size))
            return c, z, m
In [3]: size = 200
        iterations = 100

        def mandelbrot_numpy(c, z, m, iterations):
            for n in range(iterations):
                indices = np.abs(z) <= 10
                z[indices] = z[indices]**2 + c[indices]
                m[indices] = n
In [4]: %%timeit -n1 -r10 c, z, m = initialize(size)
        mandelbrot_numpy(c, z, m, iterations)
1 loops, best of 10: 245 ms per loop

在这里，Numba 超过了 NumPy。然而，我们不能仅凭这一个实验得出明确的结论。Numba 或 NumPy 哪个更快，取决于算法的具体实现、仿真参数、机器特性等因素。

这里有一些参考资料：

Numba 文档，请见 numba.pydata.org/doc.html
Numba 支持的类型，请见 numba.pydata.org/numba-doc/dev/types.html
Numba 示例请见 numba.pydata.org/numba-doc/dev/examples.html
Blaze 可用，请见 blaze.pydata.org
Theano 可用，请见 deeplearning.net/software/theano/

另见

使用 Numexpr 加速数组计算 的方案

使用 Numexpr 加速数组计算

Numexpr 是一个改进了 NumPy 弱点的包；复杂数组表达式的评估有时会很慢。原因是，在中间步骤中会创建多个临时数组，这在处理大数组时不是最优的。Numexpr 评估涉及数组的代数表达式，对其进行解析、编译，最后比 NumPy 更快地执行它们。

这个原理与 Numba 有些相似，因为普通的 Python 代码是通过 JIT 编译器动态编译的。然而，Numexpr 只处理代数数组表达式，而不是任意的 Python 代码。我们将在本教程中看到它是如何工作的。

准备工作

你可以在 github.com/pydata/numexpr 的文档中找到安装 Numexpr 的说明。

如何实现……

让我们导入 NumPy 和 Numexpr：

In [1]: import numpy as np
        import numexpr as ne

然后，我们生成三个大向量：

In [2]: x, y, z = np.random.rand(3, 1000000)

现在，我们评估 NumPy 计算涉及我们向量的复杂代数表达式所需的时间：
```
In [3]: %timeit x + (y**2 + (z*x + 1)*3)
10 loops, best of 3: 48.1 ms per loop
```

让我们用 Numexpr 执行相同的计算。我们需要将表达式作为字符串给出：

In [4]: %timeit ne.evaluate('x + (y**2 + (z*x + 1)*3)')
100 loops, best of 3: 11.5 ms per loop

Numexpr 可以使用多个核心。在这里，我们有 2 个物理核心和 4 个虚拟线程，支持英特尔超线程技术。我们可以使用 set_num_threads() 函数指定希望 Numexpr 使用的核心数：

In [5]: ne.ncores
Out[5]: 4
In [6]: for i in range(1, 5):
            ne.set_num_threads(i)
            %timeit ne.evaluate('x + (y**2 + (z*x + 1)*3)')
10 loops, best of 3: 19.4 ms per loop
10 loops, best of 3: 14 ms per loop
10 loops, best of 3: 12.8 ms per loop
10 loops, best of 3: 11.5 ms per loop

它是如何工作的……

Numexpr 会分析数组表达式，对其进行解析，并将其编译成低级语言。Numexpr 知道 CPU 向量化指令和 CPU 缓存特性。因此，Numexpr 可以动态优化向量化计算。

Numexpr、Numba 和 Blaze 之间存在一定的重叠。我们可以预期将来这些项目之间会有一些交叉影响。

另见

使用 Numba 和即时编译加速纯 Python 代码 的食谱

使用 ctypes 在 Python 中包装 C 库

在 Python 中包装 C 库使我们能够利用现有的 C 代码，或在像 C 这样快速的语言中实现代码的关键部分。

使用 Python 调用外部编译的库相对简单。第一种方式是通过 os.system 命令调用命令行可执行文件，但这种方法不能扩展到已编译的库（在 Windows 上，动态链接库，或称 DLL）。更强大的方法是使用一个名为 ctypes 的本地 Python 模块。该模块允许我们从 Python 调用在编译库中定义的函数（这些库是用 C 编写的）。ctypes 模块负责 C 和 Python 之间的数据类型转换。此外，numpy.ctypeslib 模块提供了在外部库使用数据缓冲区的地方使用 NumPy 数组的功能。

在这个示例中，我们将用 C 重写 Mandelbrot 分形的代码，将其编译成共享库，并从 Python 调用它。

准备工作

这段代码是为 Windows 编写的，可以通过一些小修改适配到其他系统。

需要一个 C 编译器。你将在本章的介绍部分找到所有与编译器相关的指令。特别是，要使 C 编译器在 Windows 上工作，你需要在启动 IPython 笔记本之前，在 Windows 终端执行一系列指令。你将在书籍的仓库中找到包含本章代码的文件夹中的批处理脚本，里面有适当的指令。

如何操作…

第一步是用 C 编写并编译 Mandelbrot 示例。第二步是使用 ctypes 从 Python 访问库。如果你只关心如何访问一个已编译的库，可以直接跳到第 3 步，假设 mandelbrot.dll 是一个已编译的库，定义了名为 mandelbrot() 的函数。

让我们用 C 编写 Mandelbrot 分形的代码：

In [1]: %%writefile mandelbrot.c

        // Needed when creating a DLL.
        #define EXPORT __declspec(dllexport)

        #include "stdio.h"
        #include "stdlib.h"

        // This function will be available in the DLL.
        EXPORT void __stdcall mandelbrot(int size,
                                         int iterations,
                                         int *col)
        {
            // Variable declarations.
            int i, j, n, index;
            double cx, cy;
            double z0, z1, z0_tmp, z0_2, z1_2;

            // Loop within the grid.
            for (i = 0; i < size; i++)
            {
                cy = -1.5 + (double)i / size * 3;
                for (j = 0; j < size; j++)
                {
                    // We initialize the loop of the 
                    // system.
                    cx = -2.0 + (double)j / size * 3;
                    index = i * size + j;
                    // Let's run the system.
                    z0 = 0.0;
                    z1 = 0.0;
                    for (n = 0; n < iterations; n++)
                    {
                        z0_2 = z0 * z0;
                        z1_2 = z1 * z1;
                        if (z0_2 + z1_2 <= 100)
                        {
                            // Update the system.
                            z0_tmp = z0_2 - z1_2 + cx;
                            z1 = 2 * z0 * z1 + cy;
                            z0 = z0_tmp;
                            col[index] = n;
                        }
                        else
                        {
                            break;
                        }
                    }
                }
            }
        }

现在，让我们使用 Microsoft Visual Studio 的 cl.exe 将这个 C 源文件构建成一个 DLL。/LD 选项指定创建 DLL：
```
In [2]: !cl /LD mandelbrot.c
/out:mandelbrot.dll
Creating library mandelbrot.lib and object mandelbrot.exp
```

让我们使用 ctypes 访问库：

In [3]: import ctypes
In [4]: # Load the DLL file in Python.
        lb = ctypes.CDLL('mandelbrot.dll')
        lib = ctypes.WinDLL(None, handle=lb._handle)
        # Access the mandelbrot function.
        mandelbrot = lib.mandelbrot

NumPy 和 ctypes 使我们能够包装 DLL 中定义的 C 函数：

In [5]: from numpy.ctypeslib import ndpointer
In [6]: # Define the types of the output and arguments.
        mandelbrot.restype = None
        mandelbrot.argtypes = [ctypes.c_int, ctypes.c_int,
                               ndpointer(ctypes.c_int)]

要使用这个函数，我们首先需要初始化一个空数组，并将其作为参数传递给 mandelbrot() 包装函数：

In [7]: import numpy as np
        # We initialize an empty array.
        size = 200
        iterations = 100
        col = np.empty((size, size), dtype=np.int32)
        # We execute the C function.
        mandelbrot(size, iterations, col)
In [8]: %timeit mandelbrot(size, iterations, col)
100 loops, best of 3: 12.5 ms per loop

我们在脚本的最后释放库：

In [9]: from ctypes.wintypes import HMODULE
        ctypes.windll.kernel32.FreeLibrary.argtypes = [HMODULE]
        ctypes.windll.kernel32.FreeLibrary(lb._handle)

如何运作…

在 C 代码中，__declspec(dllexport) 命令声明该函数在 DLL 中可见。__stdcall 关键字声明 Windows 上的标准调用约定。

mandelbrot() 函数接受以下参数：

col 缓冲区的大小（col 值是对应点位于原点周围圆盘内的最后一次迭代）
迭代次数的数量
指针指向整数缓冲区

mandelbrot()不返回任何值；相反，它更新了通过引用传递给函数的缓冲区（它是一个指针）。

为了在 Python 中包装这个函数，我们需要声明输入参数的类型。ctypes 模块为不同的数据类型定义了常量。此外，numpy.ctypeslib.ndpointer()函数允许我们在 C 函数中需要指针的地方使用 NumPy 数组。传递给ndpointer()的参数数据类型需要与传递给函数的 NumPy 数组的数据类型相对应。

一旦函数正确地被包装，它就可以像标准 Python 函数一样调用。这里，在调用mandelbrot()之后，最初为空的 NumPy 数组被填充上了曼德布罗特分形。

还有更多……

SciPy 包含一个名为weave的模块，提供类似的功能。我们可以在 Python 字符串中编写 C 代码，然后让 weave 在运行时使用 C 编译器编译并执行它。这个模块似乎维护得不好，且似乎与 Python 3 不兼容。Cython 或 ctypes 可能是更好的选择。

ctypes 的一个更新替代品是 cffi（cffi.readthedocs.org），它可能会更快且更方便使用。你还可以参考eli.thegreenplace.net/2013/03/09/python-ffi-with-ctypes-and-cffi/。

使用 Cython 加速 Python 代码

Cython既是一种语言（Python 的超集），也是一个 Python 库。使用 Cython，我们从一个常规的 Python 程序开始，并添加有关变量类型的注释。然后，Cython 将这段代码翻译为 C，并将结果编译为 Python 扩展模块。最后，我们可以在任何 Python 程序中使用这个编译后的模块。

虽然 Python 中的动态类型会带来性能开销，但 Cython 中的静态类型变量通常会导致代码执行更快。

性能提升在 CPU 密集型程序中最为显著，特别是在紧密的 Python 循环中。相比之下，I/O 密集型程序预计从 Cython 实现中受益不大。

在这个示例中，我们将看到如何使用 Cython 加速曼德布罗特代码示例。

准备工作

需要一个 C 编译器。你将在本章的介绍部分找到所有与编译器相关的说明。

你还需要从www.cython.org安装 Cython。在 Anaconda 中，你可以在终端输入conda install cython。

如何做……

我们假设size和iterations变量已如前面的示例中定义。

要在 IPython 笔记本中使用 Cython，我们首先需要导入 IPython 提供的cythonmagic扩展：
```
In [6]: %load_ext cythonmagic
```
作为第一次尝试，我们只需在mandelbrot()函数定义之前添加%%cython魔法命令即可。内部，该单元魔法会将单元编译为独立的 Cython 模块，因此所有必需的导入操作都需要在同一单元中完成。此单元无法访问交互式命名空间中定义的任何变量或函数：
```
In [6]: %%cython
        import numpy as np
        def mandelbrot_cython(m, size, iterations):
            # The exact same content as in 
            # mandelbrot_python (first recipe of 
            # this chapter).
```

这个版本有多快？

In [7]: %%timeit -n1 -r1 m = np.zeros((size, size),
                                      dtype=np.int32)
        mandelbrot_cython(m, size, iterations)
1 loops, best of 1: 5.7 s per loop

这里几乎没有加速效果。我们需要指定 Python 变量的类型。

让我们使用类型化内存视图为 NumPy 数组添加类型信息（我们将在*如何工作……*部分进行解释）。我们还采用了稍微不同的方式来测试粒子是否已经逃离了域（if测试）：

In [8]: %%cython
        import numpy as np
        def mandelbrot_cython(int[:,::1] m,
                              int size,
                              int iterations):
            cdef int i, j, n
            cdef complex z, c
            for i in range(size):
                for j in range(size):
                    c = -2 + 3./size*j + 1j*(1.5-3./size*i)
                    z = 0
                    for n in range(iterations):
                        if z.real**2 + z.imag**2 <= 100:
                            z = z*z + c
                            m[i, j] = n
                        else:
                            break

这个新版本有多快？
```
In [9]: %%timeit -n1 -r1 m = np.zeros((size, size),
                                      dtype=np.int32)
        mandelbrot_cython(m, size, iterations)
1 loops, best of 1: 230 ms per loop
```
我们所做的只是指定了局部变量和函数参数的类型，并且在计算z的绝对值时绕过了 NumPy 的np.abs()函数。这些变化帮助 Cython 从 Python 代码生成了更优化的 C 代码。

如何工作……

cdef关键字将变量声明为静态类型的 C 变量。C 变量能加速代码执行，因为它减少了 Python 动态类型带来的开销。函数参数也可以声明为静态类型的 C 变量。

通常情况下，在紧密循环中使用的变量应当使用cdef声明。为了确保代码得到良好的优化，我们可以使用注解。只需在%%cython魔法命令后添加-a标志，未优化的行将以黄色渐变显示（白色行表示较快，黄色行表示较慢）。这一点可以通过以下截图看到。颜色的变化取决于每行相对的 Python API 调用次数。

如何工作……

Cython 中的注解

有两种方法可以使用 Cython 将 NumPy 数组声明为 C 变量：使用数组缓冲区或使用类型化内存视图。在这个示例中，我们使用了类型化内存视图。在下一个示例中，我们将介绍数组缓冲区。

类型化内存视图允许使用类似 NumPy 的索引语法高效地访问数据缓冲区。例如，我们可以使用int[:,::1]来声明一个按 C 顺序存储的二维 NumPy 数组，数组元素类型为整数，::1表示该维度是连续布局。类型化内存视图可以像 NumPy 数组一样进行索引。

然而，内存视图并不实现像 NumPy 那样的逐元素操作。因此，内存视图在紧密的for循环中充当便捷的数据容器。对于逐元素的 NumPy 类似操作，应使用数组缓冲区。

通过将调用np.abs替换为更快的表达式，我们可以显著提高性能。原因是np.abs是一个 NumPy 函数，具有一定的调用开销。它是为处理较大数组而设计的，而不是标量值。这种开销在紧密的循环中会导致性能大幅下降。使用 Cython 注解可以帮助发现这一瓶颈。

还有更多内容……

使用 Cython 从 IPython 中进行操作非常方便，尤其是通过%%cython单元魔法。然而，有时我们需要使用 Cython 创建一个可重用的 C 扩展模块。实际上，IPython 的%%cython单元魔法在后台正是这么做的。

第一步是编写一个独立的 Cython 脚本，保存在.pyx文件中。它应当与%%cython单元魔法的完整内容完全一致。

第二步是创建一个setup.py文件，我们将用它来编译 Cython 模块。以下是一个基本的setup.py文件，假设有一个mandelbrot.pyx文件：

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

setup(
    cmdclass = {'build_ext': build_ext},
    ext_modules = [Extension("mandelbrot",
                             ["mandelbrot.pyx"])]
)

第三步是用 Python 执行这个设置脚本：

In [3]: !python setup.py build_ext --inplace
running build_ext
cythoning mandelbrot.pyx to mandelbrot.c
building 'mandelbrot' extension

在构建过程中创建了两个文件：C 源文件和已编译的 Python 扩展。文件扩展名在 Windows 上是.pyd（DLL 文件），在 UNIX 上是.so：
```
In [4]: !dir mandelbrot.*
mandelbrot.c
mandelbrot.pyd
mandelbrot.pyx
```
最后，我们可以像往常一样加载已编译的模块（使用from mandelbrot import mandelbrot）。

使用这种技术，Cython 代码也可以集成到 Python 包中。以下是一些参考资料：

分发 Cython 模块，详见docs.cython.org/src/userguide/source_files_and_compilation.html
使用 Cython 进行编译，详见docs.cython.org/src/reference/compilation.html

另见

通过编写更少的 Python 代码和更多的 C 代码来优化 Cython 代码的食谱
释放 GIL 以利用多核处理器，使用 Cython 和 OpenMP的食谱

通过编写更少的 Python 代码和更多的 C 代码来优化 Cython 代码

在这个食谱中，我们将考虑一个更复杂的 Cython 示例。从一个纯 Python 中的慢实现开始，我们将逐步使用不同的 Cython 特性来加速它。

我们将实现一个非常简单的光线追踪引擎。光线追踪通过模拟光传播的物理属性来渲染场景。这种渲染方法能够生成照片级真实感的场景，但计算量非常大。

在这里，我们将渲染一个单一的球体，带有漫反射和镜面反射光照。首先，我们将给出纯 Python 版本的示例代码。然后，我们将逐步使用 Cython 加速它。

注意

代码很长，包含许多函数。我们将首先给出纯 Python 版本的完整代码。然后，我们将描述用 Cython 加速代码所需的更改。所有脚本可以在本书网站上找到。

如何操作…

首先，让我们实现纯 Python 版本：

In [1]: import numpy as np
        import matplotlib.pyplot as plt
In [2]: %matplotlib inline
In [3]: w, h = 200, 200  # Size of the window in pixels.

我们为向量创建一个归一化函数：

def normalize(x):
    # This function normalizes a vector.
    x /= np.linalg.norm(x)
    return x

现在，我们创建一个计算光线与球体交点的函数：

def intersect_sphere(O, D, S, R):
    # Return the distance from O to the intersection
    # of the ray (O, D) and the sphere (S, R), or
    # +inf if there is no intersection.
    # O and S are 3D points, D (direction) is a
    # normalized vector, R is a scalar.
    a = np.dot(D, D)
    OS = O - S
    b = 2 * np.dot(D, OS)
    c = np.dot(OS, OS) - R*R
    disc = b*b - 4*a*c
    if disc > 0:
        distSqrt = np.sqrt(disc)
        q = (-b - distSqrt) / 2.0 if b < 0 \
            else (-b + distSqrt) / 2.0
        t0 = q / a
        t1 = c / q
        t0, t1 = min(t0, t1), max(t0, t1)
        if t1 >= 0:
            return t1 if t0 < 0 else t0
    return np.inf

以下函数进行光线追踪：

 def trace_ray(rayO, rayD):
    # Find first point of intersection with the scene.
    t = intersect_sphere(rayO, rayD, position, radius)
    # No intersection?
    if t == np.inf:
        return
    # Find the point of intersection on the object.
    M = rayO + rayD * t
    N = normalize(M - position)
    toL = normalize(L - M)
    toO = normalize(O - M)
    # Ambient color.
    col = ambient
    # Diffuse color.
    col += diffuse * max(np.dot(N, toL), 0) * color
    # Specular color.
    col += specular_c * color_light * \
        max(np.dot(N, normalize(toL + toO)), 0) \
           ** specular_k
    return col

最后，主循环在以下函数中实现：

def run():
    img = np.zeros((h, w, 3))
    # Loop through all pixels.
    for i, x in enumerate(np.linspace(-1.,1.,w)):
        for j, y in enumerate(np.linspace(-1.,1.,h)):
            # Position of the pixel.
            Q[0], Q[1] = x, y
            # Direction of the ray going through the
            # optical center.
            D = normalize(Q - O)
            depth = 0
            rayO, rayD = O, D
            # Launch the ray and get the
            # color of the pixel.
            col = trace_ray(rayO, rayD)
            if col is None:
                continue
            img[h - j - 1, i, :] = np.clip(col, 0, 1)
    return img

现在，我们初始化场景并定义一些参数：

In [4]: # Sphere properties.
        position = np.array([0., 0., 1.])
        radius = 1.
        color = np.array([0., 0., 1.])
        diffuse = 1.
        specular_c = 1.
        specular_k = 50

        # Light position and color.
        L = np.array([5., 5., -10.])
        color_light = np.ones(3)
        ambient = .05

        # Camera.
        O = np.array([0., 0., -1.])  # Position.
        Q = np.array([0., 0., 0.])  # Pointing to.

让我们渲染场景：

In [5]: img = run()
In [6]: plt.imshow(img)
        plt.xticks([]); plt.yticks([])

如何操作…

使用 Python 和 Cython 进行光线追踪。左侧：这个食谱的结果。右侧：扩展版本的结果。

这个实现有多慢？

In [7]: %timeit run()
1 loops, best of 1: 3.58 s per loop

如果我们只使用 %%cython 魔法并在单元格顶部添加适当的 import numpy as np 和 cimport numpy as np 命令，我们只能获得适度的改进，仅比原来快了十分之一秒。
我们可以通过提供变量类型的信息来做得更好。由于我们在 NumPy 数组上执行向量化计算，我们不能轻易使用内存视图。因此，我们将使用数组缓冲区。首先，在 Cython 模块（或 %%cython 单元格）的最开始，我们声明 NumPy 数据类型，如下所示：
```
import numpy as np
cimport numpy as np
DBL = np.double
ctypedef np.double_t DBL_C
```
然后，我们声明一个 NumPy 数组，使用 cdef np.ndarray[DBL_C, ndim=1]（在此示例中，是一个一维的双精度浮点数数组）。这里有一个难点，因为 NumPy 数组只能在函数内部声明，不能在顶层声明。因此，我们需要稍微调整代码的整体架构，将一些数组作为函数参数传递，而不是使用全局变量。然而，即使声明了所有变量的类型，我们几乎没有获得任何加速。
在当前的实现中，由于大量对小数组（三个元素）的 NumPy 函数调用，我们遭遇了性能下降。NumPy 设计上是处理大数组的，对于这么小的数组使用它并没有多大意义。

在这个特定的情况下，我们可以尝试通过使用 C 标准库重写一些函数来绕过 NumPy。我们使用 cdef 关键字声明一个 C 风格的函数。这些函数可以带来显著的性能加速。在这里，通过将 normalize() Python 函数替换为以下 C 函数，我们获得了 2-3 倍的加速：
```
from libc.math cimport sqrt
cdef normalize(np.ndarray[DBL_C, ndim=1] x):
    cdef double n
    n = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2])
    x[0] /= n
    x[1] /= n
    x[2] /= n
    return x
```
为了获得最显著的加速，我们需要完全绕过 NumPy。那么，究竟在什么地方使用了 NumPy 呢？
- 许多变量是 NumPy 数组（主要是一维向量，包含三个元素）。
- 按元素操作会产生隐式的 NumPy API 调用。
- 我们还使用了一些 NumPy 内置函数，比如 numpy.dot()。
为了在我们的示例中绕过 NumPy，我们需要根据具体需求重新实现这些功能。第一个选择是使用原生 Python 类型（例如元组）来表示向量，并编写 C 风格的函数来实现对元组的操作（假设它们总是有三个元素）。例如，两个元组相加可以实现如下：
```
cdef tuple add(tuple x, tuple y):
    return (x[0]+y[0], x[1]+y[1], x[2]+y[2])
```
我们获得了一个有趣的加速（相比纯 Python 快了 30 倍），但是通过使用纯 C 数据类型，我们可以做到更好。
我们将定义一个纯 C 结构体，而不是使用 Python 类型来表示我们的向量。换句话说，我们不仅绕过了 NumPy，还通过转向纯 C 代码绕过了 Python。为了在 Cython 中声明表示 3D 向量的 C 结构体，我们可以使用以下代码：
```
cdef struct Vec3:
    double x, y, z
```
要创建一个新的 Vec3 变量，我们可以使用以下函数：
```
cdef Vec3 vec3(double x, double y, double z):
    cdef Vec3 v
    v.x = x
    v.y = y
    v.z = z
    return v
```
作为一个例子，以下是用于加法操作两个 Vec3 对象的函数：
```
cdef Vec3 add(Vec3 u, Vec3 v):
    return vec3(u.x + v.x, u.y + v.y, u.z + v.z)
```
代码可以更新，以利用这些快速的 C 风格函数。最后，图像可以声明为 3D 内存视图。经过这些改动，Cython 实现的速度达到了约 12 毫秒，比纯 Python 版本快了 300 倍！

总结来说，我们通过将整个实现基本重写为 C 语言，并结合改进的 Python 语法，获得了非常有趣的加速效果。

如何运作…

简要说明一下光线追踪代码是如何工作的。我们将三维场景建模为包含平面和球体等物体（此处仅有一个球体）。此外，还有一个相机和一个表示渲染图像的平面：

如何运作...

光线追踪原理（“光线追踪图”由 Henrik 提供，来自 Wikimedia Commons）

代码中有一个主循环，遍历图像的所有像素。对于每个像素，我们从相机中心发射一条光线，穿过当前像素并与场景中的物体相交，计算该光线与物体的第一个交点。然后，我们根据物体材质的颜色、光源的位置、物体在交点处的法线等因素，计算像素的颜色。这里有几个受物理启发的光照方程，描述了颜色如何依赖于这些参数。在这里，我们使用Blinn-Phong 着色模型，包括环境光、漫反射光和镜面反射光组件：

如何运作...

Blinn-Phong 着色模型（“Phong 组件”，来自 Wikimedia Commons）

当然，完整的光线追踪引擎远比我们在这个示例中实现的复杂。我们还可以模拟其他光学和光照特性，如反射、折射、阴影、景深等。也可以将光线追踪算法实现到显卡上，进行实时的真实感渲染。这里有几个参考资料：

Blinn-Phong 着色模型的 Wikipedia 页面，详见 en.wikipedia.org/wiki/Blinn-Phong_shading_model
光线追踪的 Wikipedia 页面，详见 en.wikipedia.org/wiki/Ray_tracing_(graphics)

还有更多...

尽管功能强大，Cython 仍然需要对 Python、NumPy 和 C 有较好的理解。最显著的性能提升来自于将动态类型的 Python 变量转换为静态类型的 C 变量，尤其是在紧凑的循环中。

这里有几个参考资料：

可用的 Cython 扩展类型，详情见 docs.cython.org/src/userguide/extension_types.html
我们的光线追踪引擎的扩展版本，详见 gist.github.com/rossant/6046463

另见：

使用 Cython 加速 Python 代码 配方
使用 Cython 和 OpenMP 释放 GIL 以利用多核处理器 配方

使用 Cython 和 OpenMP 释放 GIL，以充分利用多核处理器

如我们在本章介绍中所看到的，CPython 的 GIL 阻止了纯 Python 代码利用多核处理器。使用 Cython，我们可以在代码的某个部分临时释放 GIL，从而启用多核计算。这是通过OpenMP完成的，OpenMP 是一个多处理 API，大多数 C 编译器都支持它。

在本食谱中，我们将看到如何在多个核心上并行化前一个食谱的代码。

准备工作

要在 Cython 中启用 OpenMP，您只需为编译器指定一些选项。除了一个好的 C 编译器之外，您无需在计算机上安装任何特殊的软件。有关更多详细信息，请参阅本章介绍中的说明。

在本食谱中，我们使用的是微软的 Visual C++编译器（适用于 Windows），但代码可以轻松适应其他系统。

如何实现…

我们的简单光线追踪引擎实现是极其并行的；有一个遍历所有像素的主循环，在这个循环中，完全相同的函数会被重复调用。循环迭代之间没有串扰。因此，理论上可以在并行中执行所有迭代。

在这里，我们将并行执行一个循环（遍历图像中的所有列），使用 OpenMP。

您可以在书籍网站上找到完整的代码。我们这里只展示最重要的步骤：

我们在%%cython魔法命令中添加以下选项：--compile-args=/openmp --link-args=/openmp。具体语法可能依赖于您的编译器和/或系统。例如，/openmp应该替换为 GCC 中的-fopenmp。
我们导入了prange()函数：
```
from cython.parallel import prange
```
我们在每个函数定义后添加nogil，以移除 GIL。在标记为nogil的函数中，不能使用任何 Python 变量或函数。例如：
```
cdef Vec3 add(Vec3 x, Vec3 y) nogil:
    return vec3(x.x + y.x, x.y + y.y, x.z + y.z)
```
要使用 OpenMP 在多个核心上并行运行循环，我们使用prange()：
```
with nogil:
    for i in prange(w):
        # ...
```
在使用任何并行计算功能（如prange()）之前，需要释放 GIL。
通过这些更改，我们在四核处理器上达到了 4 倍的加速。

原理…

GIL（全局解释器锁）已经在本章的介绍中描述过。nogil关键字告诉 Cython，某个特定的函数或代码段应该在没有 GIL 的情况下执行。当 GIL 被释放时，无法进行任何 Python API 调用，这意味着只能使用 C 变量和 C 函数（用cdef声明的函数）。

另见

用 Cython 加速 Python 代码 食谱
通过编写更少的 Python 代码和更多的 C 代码来优化 Cython 代码 食谱
通过 IPython 将 Python 代码分配到多个核心 食谱

为 NVIDIA 显卡（GPU）编写大规模并行代码，使用 CUDA

图形处理单元（GPU）是专门用于实时渲染的强大处理器。我们几乎在任何计算机、笔记本、视频游戏主机、平板或智能手机中都能找到 GPU。它们的并行架构包含数十到数千个核心。视频游戏产业在过去二十年里一直在推动更强大 GPU 的发展。

GPUs（图形处理单元）已被广泛应用于现代超级计算机中（例如，位于橡树岭国家实验室的 Cray Titan，约 20 佩塔 FLOPS，约 20,000 个 CPU，以及大量的 NVIDIA GPUs）。今天，一个高端 $1000 的 GPU 大致相当于 2000 年的一个价值$ 1 亿的超级计算机（约几个 teraFLOPS）。

注意

FLOPS 指的是每秒浮点运算次数（FLoating-point Operations Per Second）。一个 1 teraFLOPS 的 GPU 每秒可以执行高达一万亿次浮点运算。

自 2000 年代中期以来，GPU 已不再局限于图形处理。我们现在可以在 GPU 上实现科学算法。唯一的条件是算法必须遵循 SIMD（单指令多数据）范式，即一系列指令并行执行多个数据。这被称为 图形处理单元上的通用计算（GPGPU）。GPGPU 已被应用于多个领域：气象学、数据挖掘、计算机视觉、图像处理、金融、物理学、生物信息学等多个领域。为 GPU 编写代码可能具有挑战性，因为它需要理解硬件的内部架构。

CUDA 是由 NVIDIA 公司于 2007 年创建的专有 GPGPU 框架，是主要的 GPU 制造商之一。用 CUDA 编写的程序仅能在 NVIDIA 显卡上运行。还有另一个竞争性的 GPGPU 框架，称为 OpenCL，这是一个由其他主要公司支持的开放标准。OpenCL 程序可以在大多数厂商的 GPU 和 CPU 上运行（特别是 NVIDIA、AMD 和英特尔）。

在这个示例中，我们将展示一个非常基础的 GPGPU 示例。我们将使用 CUDA 实现 Mandelbrot 分形的非常简单的并行计算。在下一个示例中，我们将使用 OpenCL 实现完全相同的例子。

提示

对于一个新项目，你应该选择 OpenCL 还是 CUDA？答案主要取决于你用户群体的硬件配置。如果你需要在实验室中针对所有配备 NVIDIA 显卡的计算机实现尽可能高的性能，并且发布你的程序到全球并不是优先考虑的事情，你可以选择 CUDA。如果你计划将程序分发给使用不同平台的多人，你应该选择 OpenCL。在功能上，这两者大体上是相当的。

我们可以通过 PyCUDA 在 Python 中使用 CUDA，这是由 Andreas Klöckner 编写的一个 Python 包 (documen.tician.de/pycuda/)。

准备工作

安装和配置 PyCUDA 通常不简单。首先，您需要一块 NVIDIA GPU。然后，您需要安装 CUDA SDK。最后，您需要安装并配置 PyCUDA。请注意，PyCUDA 依赖于一些外部包，特别是 pytools。

在 Windows 上，您应该使用 Chris Gohlke 的包。确保您的 CUDA 版本与 PyCUDA 包中使用的版本匹配。如果遇到 DLL 相关的问题，使用 Dependency Walker 检查 PyCUDA 安装文件夹中的*.pyd文件（在 Anaconda 下，它应该位于C:\anaconda\lib\site-packages\pycuda）。如果您使用的是 Windows 64 位，确保C:\Windows\SysWOW64在您的系统 PATH 中。最后，确保您有与 Python 版本相对应的 Visual Studio 版本（请参阅本章开头关于 C 编译器的说明）。

您可以在以下链接找到更多信息：

可在developer.nvidia.com/cuda-downloads下载 CUDA SDK
PyCUDA 维基可以在wiki.tiker.net找到

如何做……

让我们导入 PyCUDA：

In [1]: import pycuda.driver as cuda
        import pycuda.autoinit
        from pycuda.compiler import SourceModule
        import numpy as np

我们初始化一个将包含分形的 NumPy 数组：

In [2]: size = 200
        iterations = 100
        col = np.empty((size, size), dtype=np.int32)

我们为这个数组分配 GPU 内存：

In [3]: col_gpu = cuda.mem_alloc(col.nbytes)

我们将 CUDA 内核写入字符串中。mandelbrot()函数的参数包括：

图像的大小
迭代次数
指针指向内存缓冲区

这个函数在每个像素上执行。它更新col缓冲区中的像素颜色：

In [4]: code = """
        __global__ void mandelbrot(int size,
                                   int iterations,
                                   int *col) {
            // Get the row and column of the thread.
            int i = blockIdx.y * blockDim.y + threadIdx.y;
            int j = blockIdx.x * blockDim.x + threadIdx.x;
            int index = i * size + j;

            // Declare and initialize the variables.
            double cx, cy;
            double z0, z1, z0_tmp, z0_2, z1_2;
            cx = -2.0 + (double)j / size * 3;
            cy = -1.5 + (double)i / size * 3;

            // Main loop.
            z0 = z1 = 0.0;
            for (int n = 0; n < iterations; n++)
            {
                z0_2 = z0 * z0;
                z1_2 = z1 * z1;
                if (z0_2 + z1_2 <= 100)
                {
                    // Need to update both z0 and z1,
                    // hence the need for z0_tmp.
                    z0_tmp = z0_2 - z1_2 + cx;
                    z1 = 2 * z0 * z1 + cy;
                    z0 = z0_tmp;
                    col[index] = n;
                }
                else break;
            }
        }
        """

现在，我们编译 CUDA 程序：

In [5]: prg = SourceModule(code)
        mandelbrot = prg.get_function("mandelbrot")

我们定义了块大小和网格大小，指定线程如何根据数据进行并行化：

In [6]: block_size = 10
        block = (block_size, block_size, 1)
        grid = (size // block_size, size // block_size, 1)

我们调用编译后的函数：

In [7]: mandelbrot(np.int32(size),
                   np.int32(iterations),
                   col_gpu,
                   block=block, grid=grid)

一旦函数执行完毕，我们将 CUDA 缓冲区的内容复制回 NumPy 数组col：
```
In [8]: cuda.memcpy_dtoh(col, col_gpu)
```
现在，col数组包含了曼德尔布罗特分形。我们发现这个 CUDA 程序在移动 GeForce GPU 上执行时间为 0.7 毫秒。

它是如何工作的……

GPU 编程是一个丰富且高度技术性的主题，涉及 GPU 的低级架构细节。当然，我们这里只用最简单的范式（“极其并行”的问题）浅尝辄止。我们将在后续部分提供更多参考。

一个 CUDA GPU 有多个多处理器，每个多处理器有多个流处理器（也称为CUDA 核心）。每个多处理器与其他多处理器并行执行。在一个多处理器中，流处理器在相同时间执行相同的指令，但操作不同的数据位（SIMD 范式）。

CUDA 编程模型的核心概念包括内核、线程、块和网格：

内核是用类似 C 的语言编写的程序，运行在 GPU 上。
线程代表在一个流处理器上执行的一个内核。
一个块包含多个在一个多处理器上执行的线程。
一个网格包含多个块。

每个块的线程数受多处理器大小的限制，并且取决于显卡型号（在最近的型号中为 1024）。然而，一个网格可以包含任意数量的块。

在一个块内，线程通常以 warp（通常是 32 个线程）的形式执行。当内核中的条件分支组织成 32 个线程一组时，性能会更好。

块内的线程可以通过 __syncthreads() 函数在同步屏障处进行同步。此功能允许一个块内的线程之间进行通信。然而，块是独立执行的，因此来自不同块的两个线程不能同步。

在一个块内，线程被组织成 1D、2D 或 3D 结构，网格中的块也是如此，如下图所示。这个结构非常方便，因为它与现实世界问题中常见的多维数据集相匹配。

它是如何工作的…

CUDA 编程模型（显示线程、块和网格 — 图像来自 NVIDIA 公司）

内核可以检索块内的线程索引（threadIdx）以及网格内的块索引（blockIdx），以确定它应该处理的数据部分。在这个配方中，分形的 2D 图像被划分为 10 x 10 的块，每个块包含 100 个像素，每个像素对应一个线程。内核 mandelbrot 计算单个像素的颜色。

GPU 上有多个内存层次，从块内少数线程共享的小而快的本地内存，到所有块共享的大而慢的全局内存。我们需要调整代码中的内存访问模式，以匹配硬件约束并实现更高的性能。特别地，当 warp 内的线程访问全局内存中的连续地址时，数据访问效率更高；硬件会将所有内存访问合并为对连续 DRAM（动态随机存取内存）位置的单次访问。

PyCUDA 让我们可以将数据从 NumPy 数组上传/下载到驻留在 GPU 上的缓冲区。这项操作通常非常耗费资源。复杂的现实世界问题常常涉及在 CPU 和 GPU 上都进行迭代步骤，因此二者之间的通信成为常见的性能瓶颈。当这些交换较少时，性能可以得到提升。

在 (Py)CUDA 的 C/Python 端有一些模板代码，涉及初始化 GPU、分配数据、将数据上传/下载到/从 GPU、编译内核、执行内核等等。你可以在 CUDA/PyCUDA 文档中找到所有的详细信息，但作为一种初步方法，你也可以直接复制并粘贴这个配方或任何教程中的代码。

还有更多内容…

这里有一些参考资料：

官方 CUDA 门户：developer.nvidia.com/category/zone/cuda-zone
CUDA 的教育与培训，详情请见developer.nvidia.com/cuda-education-training
关于 CUDA 的推荐书籍，详见developer.nvidia.com/suggested-reading
关于选择 CUDA 还是 OpenCL，详见wiki.tiker.net/CudaVsOpenCL
关于 CUDA 和 OpenCL 的博客文章，详见streamcomputing.eu/blog/2011-06-22/opencl-vs-cuda-misconceptions/

另见

为异构平台编写大规模并行代码与 OpenCL食谱

为异构平台编写大规模并行代码与 OpenCL

在前一个食谱中，我们介绍了 CUDA，这是由 NVIDIA 公司创建的专有GPGPU 框架。在这个食谱中，我们介绍了 OpenCL，这是一个由苹果公司在 2008 年启动的开放框架。现在，主流公司，包括 Intel、NVIDIA、AMD、Qualcomm、ARM 等，均已采用。它们都属于非盈利的技术联盟Khronos Group（该联盟还维护着 OpenGL 实时渲染规范）。用 OpenCL 编写的程序可以在 GPU 和 CPU 上运行（异构计算）。

提示

在概念、语法和特性上，CUDA 和 OpenCL 相对相似。由于 CUDA 的 API 与硬件的匹配度比 OpenCL 的通用 API 更高，CUDA 有时能带来稍微更高的性能。

我们可以通过PyOpenCL在 Python 中使用 OpenCL，这是 Andreas Klöckner 编写的一个 Python 包(documen.tician.de/pyopencl/)。

在这个食谱中，我们将使用 OpenCL 实现曼德布罗特分形。OpenCL 内核与前一个食谱中的 CUDA 内核非常相似。用于访问 OpenCL 的 Python API 与 PyCUDA 略有不同，但概念是相同的。

准备工作

安装 PyOpenCL 通常并不简单。第一步是安装适用于你硬件（CPU 和/或 GPU）的 OpenCL SDK。然后，你需要安装并配置 PyOpenCL。在 Windows 上，你应使用 Chris Gohlke 的包。上一食谱中的一些安装说明在这里同样适用。

这里有一些参考资料：

PyOpenCL Wiki 可在wiki.tiker.net访问
PyOpenCL 的文档可在documen.tician.de/pyopencl/查看

这里是各个 OpenCL SDK 的链接：

Intel 的 SDK 可在software.intel.com/en-us/vcsource/tools/opencl-sdk下载
AMD 的 SDK 可在developer.amd.com/tools-and-sdks/heterogeneous-computing/下载
NVIDIA 的 SDK 可在developer.nvidia.com/opencl下载

如何操作…

让我们导入 PyOpenCL：

In [1]: import pyopencl as cl
        import numpy as np

以下对象定义了一些与设备内存管理相关的标志：
```
In [2]: mf = cl.mem_flags
```

我们创建一个 OpenCL 上下文和一个命令队列：

In [3]: ctx = cl.create_some_context()
        queue = cl.CommandQueue(ctx)

现在，我们初始化将包含分形的 NumPy 数组：

In [4]: size = 200
        iterations = 100
        col = np.empty((size, size), dtype=np.int32)

我们为这个数组分配 GPU 内存：

In [5]: col_buf = cl.Buffer(ctx,
                            mf.WRITE_ONLY,
                            col.nbytes)

我们将 OpenCL 内核写入字符串中：

In [6]: code = """
        __kernel void mandelbrot(int size,
                                 int iterations,
                                 global int *col)
        {
            // Get the row and column index of the thread.
            int i = get_global_id(1);
            int j = get_global_id(0);
            int index = i * size + j;

            // Declare and initialize the variables.
            double cx, cy;
            double z0, z1, z0_tmp, z0_2, z1_2;
            cx = -2.0 + (double)j / size * 3;
            cy = -1.5 + (double)i / size * 3;

            // Main loop.
            z0 = z1 = 0.0;
            for (int n = 0; n < iterations; n++)
            {
                z0_2 = z0 * z0;
                z1_2 = z1 * z1;
                if (z0_2 + z1_2 <= 100)
                {
                    // Need to update both z0 and z1.
                    z0_tmp = z0_2 - z1_2 + cx;
                    z1 = 2 * z0 * z1 + cy;
                    z0 = z0_tmp;
                    col[index] = n;
                }
                else break;
            }
        }
        """

现在，我们编译 OpenCL 程序：

In [7]: prg = cl.Program(ctx, code).build()
Build on <pyopencl.Device 'Intel(R) Core(TM) i3-2365M CPU @ 1.40GHz' on 'Intel(R) OpenCL' at 0x765b188> succeeded.

我们调用已编译的函数，并将命令队列、网格大小和缓冲区作为参数传递：

In [8]: prg.mandelbrot(queue, col.shape, None, np.int32(size), np.int32(iterations), col_buf).wait()

一旦函数执行完成，我们将 OpenCL 缓冲区的内容复制回 NumPy 数组 col：
```
In [9]: cl.enqueue_copy(queue, col, col_buf)
```
最后，我们可以通过 imshow() 显示 NumPy 数组 col 来检查函数是否成功。我们还可以使用 %timeit 进行快速基准测试，结果显示该函数在 Intel i3 双核 CPU 上大约需要 3.7 毫秒完成。

它是如何工作的……

前面食谱中详细介绍的原理同样适用于此。CUDA 和 OpenCL 之间存在术语上的差异：

CUDA 线程相当于 OpenCL 的工作项。
CUDA 块相当于 OpenCL 的工作组。
CUDA 网格相当于 OpenCL 的NDRange。
CUDA 流处理器相当于 OpenCL 的计算单元。

在内核中，我们可以使用 get_local_id(dim)、get_group_id(dim) 和 get_global_id(dim) 获取工作项的索引。函数参数中的 global 限定符表示某个变量对应于全局内存中的对象。

OpenCL 上下文是工作项执行的环境。它包括带有内存和命令队列的设备。命令队列是主机应用程序用来将工作提交到设备的队列。

该程序在 CPU 或 GPU 上的表现相同，取决于安装的 OpenCL SDK 和可用的 OpenCL 上下文。如果存在多个上下文，PyOpenCL 可能会要求用户选择设备。上下文也可以通过编程方式指定（参见 documen.tician.de/pyopencl/runtime.html#pyopencl.Context）。在 CPU 上，代码通过多核和使用 SSE 或 AVX 等向量指令进行并行化和向量化。

还有更多内容……

OpenCL 是一个相对年轻的标准，但我们应该预期它在未来会变得越来越重要。它得到了 GPU 行业内最大公司的支持。它支持与 OpenGL 的互操作性，OpenGL 是实时硬件加速计算机图形的行业标准（由同样的 Khronos Group 维护）。它正在逐步支持移动平台（智能手机和平板电脑），并且在浏览器中也在逐步支持，使用WebCL（在撰写本文时仍为草案阶段）。

这里是一些 OpenCL 资源：

OpenCL 教程可用：opencl.codeplex.com
可用的课程列表：developer.amd.com/partners/university-programs/opencl-university-course-listings/
关于 OpenCL 的书籍，见 streamcomputing.eu/knowledge/for-developers/books/

另见

*为 NVIDIA 显卡（GPU）编写大规模并行代码（CUDA）*的配方

使用 IPython 将 Python 代码分发到多个核心。

尽管 CPython 的 GIL 存在，仍然可以通过使用多个进程而不是多个线程，在多核计算机上并行执行多个任务。Python 提供了一个本地的multiprocessing模块。IPython 提供了一个更简单的界面，带来了强大的并行计算功能，并且支持交互式环境。我们将在这里描述这个工具。

如何操作…

首先，我们在单独的进程中启动四个 IPython 引擎。我们基本上有两个选项来做到这一点：
- 在系统 shell 中执行ipcluster start -n 4。
- 通过点击 Clusters 标签并启动四个引擎，使用 IPython 笔记本主页提供的网页界面。
然后，我们创建一个客户端，作为 IPython 引擎的代理。该客户端会自动检测正在运行的引擎：
```
In [2]: from IPython.parallel import Client
        rc = Client()
```
让我们检查一下正在运行的引擎数量：
```
In [3]: rc.ids
Out[3]: [0, 1, 2, 3]
```

要在多个引擎上并行运行命令，我们可以使用%px行魔法命令或%%px单元格魔法命令：

In [4]: %%px
        import os
        print("Process {0:d}.".format(os.getpid()))
[stdout:0] Process 2988.
[stdout:1] Process 5192.
[stdout:2] Process 4484.
[stdout:3] Process 1360.

我们可以使用--targets或-t选项指定在哪些引擎上运行命令：

In [5]: %%px -t 1,2
        # The os module has already been imported in
        # the previous cell.
        print("Process {0:d}.".format(os.getpid()))
[stdout:1] Process 5192.
[stdout:2] Process 4484.

默认情况下，%px魔法命令在阻塞模式下执行命令；只有当所有引擎上的命令完成时，单元格才会返回。也可以使用--noblock或-a选项运行非阻塞命令。在这种情况下，单元格会立即返回，任务的状态和结果可以从 IPython 的交互式会话中异步轮询：
```
In [6]: %%px -a
        import time
        time.sleep(5)
Out[6]: <AsyncResult: execute>
```
之前的命令返回了一个ASyncResult实例，我们可以用它来轮询任务的状态：
```
In [7]: print(_.elapsed, _.ready())
(0.061, False)
```

%pxresult会阻塞，直到任务完成：

In [8]: %pxresult
In [9]: print(_.elapsed, _.ready())
(5.019, True)

IPython 提供了适用于常见用例的方便函数，例如并行的map函数：

In [10]: v = rc[:]
         res = v.map(lambda x: x*x, range(10))
In [11]: print(res.get())
[0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

它是如何工作的…

将代码分发到多个核心的步骤如下：

启动多个 IPython 引擎（通常每个核心对应一个进程）。
创建一个作为这些引擎代理的Client。
使用客户端在引擎上启动任务并获取结果。

引擎是执行代码的 Python 进程，它们运行在不同的计算单元上。它们与 IPython 内核非常相似。

有两种主要的接口用于访问引擎：

使用直接接口，我们可以直接并显式地通过它们的标识符访问引擎。
使用负载均衡接口，我们通过一个接口访问引擎，该接口会自动并动态地将任务分配到合适的引擎。

我们也可以为替代的并行方式创建自定义接口。

在这个例子中，我们使用了直接接口；通过在%px魔法命令中明确指定引擎的标识符，我们显式地访问了各个引擎。

正如我们在这个示例中所看到的，任务可以同步或异步启动。%px* 魔法命令在笔记本中尤其方便，因为它们让我们能够在多个引擎上无缝并行工作。

还有更多…

IPython 的并行计算功能为我们提供了一种简单的方式，可以在多个核心上并行启动独立的任务。一个更高级的用例是任务之间有依赖的情况。

依赖并行任务

依赖有两种类型：

功能依赖：它决定了给定任务是否可以在特定引擎上执行，依据引擎的操作系统、特定 Python 模块的存在与否或其他条件。IPython 提供了一个 @require 装饰器，用于那些需要特定 Python 模块才能在引擎上运行的函数。另一个装饰器是 @depend，它让我们可以定义任意的条件，这些条件在 Python 函数中实现，并返回 True 或 False。
图依赖：它决定了给定任务是否可以在给定时间、特定引擎上执行。我们可能要求某个任务仅在一个或多个其他任务完成后才能执行。此外，我们可以在任何单独的引擎中强制执行这个条件；某个引擎可能需要先执行一组特定任务，才会执行我们的任务。例如，下面是如何要求任务 B 和 C（它们的异步结果分别为 arB 和 arC）在任务 A 启动前完成：
```
with view.temp_flags(after=[arB, arC]):
    arA = view.apply_async(f)
```

IPython 提供了选项，指定任务运行所需的依赖是否全部或部分满足。此外，我们还可以指定是否应将依赖于成功和/或失败的任务视为满足条件。

当一个任务的依赖未满足时，调度器会将任务重新分配到一个引擎，再到另一个引擎，依此类推，直到找到合适的引擎。如果无法在任何引擎上满足依赖条件，任务将引发 ImpossibleDependency 错误。

在 IPython.parallel 中，将数据传递给依赖任务并不特别容易。一种可能性是在交互式会话中使用阻塞调用；等待任务完成，获取结果并将其发送回下一个任务。另一种可能性是通过文件系统在引擎之间共享数据，但这种方案在多个计算机上效果不好。一个替代方案可以参考：nbviewer.ipython.org/gist/minrk/11415238。

替代的并行计算解决方案

除了 IPython，还有许多其他的 Python 并行计算框架，包括 ParallelPython、joblib 等等。

还有一些第三方（通常是商业）服务提供基于 Python 的云服务，如 PythonAnywhere 或 Wakari。它们通常有两种使用方式：

将大量计算任务分配到多个节点并行执行：我们可以使用数百或数千台服务器进行并行计算，而无需担心整个基础设施的维护（由公司负责处理）。
在线托管 Python 应用程序，通常带有 Web 界面：例如，使用 Wakari，IPython 笔记本可以在云端运行。一个有趣的应用场景是教学，学生可以通过连接到互联网的 Web 浏览器即时使用 IPython，而无需在本地安装任何软件。

参考资料

以下是关于 IPython.parallel 的一些参考资料：

IPython.parallel 文档，见ipython.org/ipython-doc/dev/parallel/
IPython 开发者提供的 IPython 并行教程，见nbviewer.ipython.org/github/minrk/IPython-parallel-tutorial/blob/master/Index.ipynb
IPython.parallel 中的依赖关系解释，见ipython.org/ipython-doc/dev/parallel/parallel_task.html#dependencies
DAG 依赖关系的描述，见ipython.org/ipython-doc/dev/parallel/dag_dependencies.html
有关 IPython.parallel 的高级技巧示例，见github.com/ipython/ipython/tree/master/examples/Parallel%20Computing

以下是关于 Python 中替代并行计算解决方案的一些参考资料：

Parallel Python，见www.parallelpython.com
Joblib，见pythonhosted.org/joblib/parallel.html
可用的并行计算包列表，见wiki.python.org/moin/ParallelProcessing
Python Anywhere，见www.pythonanywhere.com
Wakari，见wakari.io
在 Wakari 上使用 IPCluster 的介绍，见continuum.io/blog/ipcluster-wakari-intro
使用 Wakari 进行教学的介绍，见continuum.io/blog/teaching-with-wakari

另见

在 IPython 中与异步并行任务交互的示例
在 IPython 中使用 MPI 并行化代码的示例

在 IPython 中与异步并行任务交互

在这个示例中，我们将展示如何与在 IPython 中并行运行的异步任务交互。

准备工作

你需要启动 IPython 引擎（参见前面的步骤）。最简单的选项是从笔记本仪表板的 Clusters 标签页启动它们。在这个步骤中，我们使用了四个引擎。

如何做…

让我们导入一些模块：

In [1]: import time
        import sys
        from IPython import parallel
        from IPython.display import clear_output, display
        from IPython.html import widgets

我们创建一个 Client：
```
In [2]: rc = parallel.Client()
```
现在，我们创建一个基于 IPython 引擎的负载均衡视图：
```
In [3]: view = rc.load_balanced_view()
```

我们为并行任务定义一个简单的函数：

In [4]: def f(x):
            import time
            time.sleep(.1)
            return x*x

我们将在 100 个整数上并行运行此函数：
```
In [5]: numbers = list(range(100))
```
我们使用 map_async() 在所有引擎上并行执行 f() 函数，作用对象是我们的 numbers 列表。该函数会立即返回一个 AsyncResult 对象，允许我们交互式地检索任务的相关信息：
```
In [6]: ar = view.map_async(f, numbers)
```

该对象有一个 metadata 属性：所有引擎的字典列表。我们可以获取提交和完成日期、状态、标准输出和错误以及其他信息：

In [7]: ar.metadata[0]
Out[7]: {
 'execute_result': None,
 'engine_id': None,
...
 'submitted': datetime.datetime(2014, 1, 1, 10, 30, 38, 9487),
 'follow': None}

遍历 AsyncResult 实例时工作正常；迭代实时进行，任务在完成时进展：

In [8]: for _ in ar:
            print(_, end=', ')
0, 1, 4,..., 9409, 9604, 9801,

现在，我们为异步任务创建一个简单的进度条。其思路是创建一个循环，每秒轮询任务状态。IntProgressWidget 小部件实时更新，显示任务的进度：

In [9]: def progress_bar(ar):
            # We create a progress bar.
            w = widgets.IntProgressWidget()
            # The maximum value is the number of tasks.
            w.max = len(ar.msg_ids)
            # We display the widget in the output area.
            display(w)
            # Repeat every second:
            while not ar.ready():
                # Update the widget's value with the
                # number of tasks that have finished
                # so far.
                w.value = ar.progress
                time.sleep(1)
            w.value = w.max
In [10]: ar = view.map_async(f, numbers)
In [11]: progress_bar(ar)

以下截图显示了进度条：

如何做…

最后，调试引擎上的并行任务非常容易。我们可以通过在 %%px 单元魔法中调用 %qtconsole 来启动远程内核上的 Qt 客户端：
```
In [12]: %%px -t 1
         %qtconsole
```
Qt 控制台允许我们检查远程命名空间进行调试或分析，如下图所示：

用于调试 IPython 引擎的 Qt 控制台

它是如何工作的…

AsyncResult 实例由异步并行函数返回。它们实现了几个有用的属性和方法，特别是：

elapsed: 自提交以来的经过时间
progress: 到目前为止已完成的任务数量
serial_time: 所有并行任务计算时间的总和
metadata: 包含任务更多信息的字典
ready(): 返回调用是否已完成
successful(): 返回调用是否没有抛出异常（如果任务尚未完成则抛出异常）
wait(): 阻塞直到任务完成（有一个可选的超时参数）
get(): 阻塞直到任务完成并返回结果（有一个可选的超时参数）

还有更多…

这里是一些参考资料：

AsyncResult 类的文档请参阅 ipython.org/ipython-doc/dev/parallel/asyncresult.html
任务接口的文档请参阅 ipython.org/ipython-doc/dev/parallel/parallel_task.html
实时打印引擎输出，示例代码可在github.com/ipython/ipython/blob/master/examples/Parallel%20Computing/iopubwatcher.py中找到

另见

使用 IPython 跨多个核心分发 Python 代码的配方
使用 MPI 在 IPython 中并行化代码的配方

使用 MPI 在 IPython 中并行化代码

消息传递接口（MPI）是一个用于并行系统的标准化通信协议。它在许多并行计算应用中被用来在节点之间交换数据。MPI 的入门门槛较高，但它非常高效且功能强大。

IPython 的并行计算系统从底层设计时就考虑到了与 MPI 的兼容。如果你是 MPI 的新手，从 IPython 开始使用它是个不错的选择。如果你是有经验的 MPI 用户，你会发现 IPython 与你的并行应用无缝集成。

在这个配方中，我们将看到如何通过一个非常简单的例子，使用 MPI 与 IPython 结合。

准备工作

要在 IPython 中使用 MPI，你需要：

一个标准的 MPI 实现，例如OpenMPI或MPICH
mpi4py 包可以在mpi4py.scipy.org找到

例如，下面是 Ubuntu 和 Anaconda 上安装 MPI 的命令：

conda install mpich2
conda install mpi4py

你还可以使用pip install mpi4py来安装 mpi4py。MPI 也可以在 Windows 上使用。可以参考Python Tools for Visual Studio的网站，网址为pytools.codeplex.com，该网站提供了相关的安装指导。

如何实现……

我们首先需要创建一个 MPI 配置文件，命令为：
```
In [1]: !ipython profile create --parallel --profile=mpi
```
然后，我们打开~/.ipython/profile_mpi/ipcluster_config.py文件并添加一行：c.IPClusterEngines.engine_launcher_class = 'MPI'。
一旦 MPI 配置文件创建并配置好，我们可以通过在终端中输入：ipcluster start -n 2 --engines MPI --profile=mpi来启动引擎。

现在，为了真正使用引擎，我们在笔记本中创建一个客户端：

In [2]: import numpy as np
        from IPython.parallel import Client
In [3]: c = Client(profile='mpi')

让我们在所有引擎上创建一个视图：
```
In [4]: view = c[:]
```
在这个例子中，我们通过两个核心并行计算 0 到 15 之间所有整数的总和。我们首先将包含 16 个值的数组分配给各个引擎（每个引擎获取一个子数组）：
```
In [5]: view.scatter('a', np.arange(16., dtype='float'))
Out[5]: <AsyncResult: scatter>
```

我们使用 MPI 的allreduce()函数并行计算总和。每个节点执行相同的计算并返回相同的结果：

In [6]: %%px
        from mpi4py import MPI
        import numpy as np
        print MPI.COMM_WORLD.allreduce(np.sum(a), op=MPI.SUM)
[stdout:0] 120.0
[stdout:1] 120.0

提示

如果你得到了不同的结果，意味着引擎实际上并没有使用 MPI 启动（请参见stackoverflow.com/a/20159018/1595060）。

它是如何工作的……

在这个例子中，每个节点：

接收整数的一个子集
计算这些整数的局部和
将这个局部和发送到所有其他引擎
接收其他引擎的局部和
计算这些局部和的总和

这是allreduce()在 MPI 中的工作原理；其原理是首先分发数据到各个引擎，然后通过一个全局操作符（此处为MPI.SUM）减少本地计算。

IPython 的直接接口也本地支持 scatter/gather 范式，而不需要使用 MPI。然而，这些操作只能在交互式会话中启动，而不能在引擎本身启动。

MPI 中还有许多其他的并行计算范式。你可以在这里找到更多信息：

Wes Kendall 的 MPI 教程，可以在mpitutorial.com找到
Blaise Barney（劳伦斯·利弗莫尔国家实验室）的 MPI 教程，可以在computing.llnl.gov/tutorials/mpi/找到

另见

使用 IPython 将 Python 代码分布到多个核心的食谱
在 IPython 中与异步并行任务交互的食谱

在笔记本中尝试 Julia 语言

Julia（julialang.org）是一种年轻的高层次动态语言，专为高性能数值计算而设计。它的第一个版本在 2012 年发布，在 MIT 进行了三年的开发。Julia 借鉴了 Python、R、MATLAB、Ruby、Lisp、C 等语言的思想。它的主要优势是将高层次动态语言的表现力和易用性与 C 语言（几乎）的速度相结合。这是通过基于 LLVM 的即时编译器（JIT）实现的，目标是 x86-64 架构的机器代码。

在这个食谱中，我们将使用IJulia包在 IPython 笔记本中尝试 Julia，IJulia 包可以在github.com/JuliaLang/IJulia.jl找到。我们还将展示如何从 Julia 使用 Python 包（如 NumPy 和 matplotlib）。具体来说，我们将计算并显示一个 Julia 集。

这个食谱灵感来自 David P. Sanders 在 SciPy 2014 会议上提供的 Julia 教程（nbviewer.ipython.org/github/dpsanders/scipy_2014_julia/tree/master/）。

准备工作

首先，你需要安装 Julia。你可以在 Julia 官网julialang.org/downloads/上找到适用于 Windows、Mac OS X 和 Linux 的安装包。在 Ubuntu 上，你可以在终端输入sudo apt-get install julia。对于 IJulia，你还需要一个 C++编译器。在 Ubuntu 上，你可以输入sudo apt-get install build-essential。

然后，使用julia命令打开一个 Julia 终端，并在 Julia 终端中输入Pkg.add("IJulia")来安装 IJulia。这个包还会在你的 IPython 安装中创建一个julia配置文件。

最后，要启动 Julia 笔记本，请在终端中运行ipython notebook --profile=julia。你将看到 IPython 笔记本的仪表板。唯一的区别是，在笔记本中使用的是 Julia 语言，而不是 Python。

这个食谱已在 Ubuntu 14.04 和 Julia 0.2.1 版本上进行测试。

如何实现…

我们无法避免传统的Hello World示例。println()函数用于显示一个字符串，并在末尾添加换行符：
```
In [1]: println("Hello world!")
Hello world!
```
我们创建了一个多态函数f，计算表达式z*z+c。我们将在数组上评估这个函数，因此我们使用带点（.）前缀的逐元素操作符：
```
In [2]: f(z, c) = z.*z .+ c
Out[2]: f (generic function with 1 method)
```
让我们在标量复数上评估f（虚数* i *是1im）。
```
In [3]: f(2.0 + 1.0im, 1.0)
Out[3]: 4.0 + 4.0im
```
现在，我们创建一个(2, 2)矩阵。组件之间用空格隔开，行之间用分号（;）隔开。这个Array的类型会根据其组件自动推断出来。Array类型是 Julia 中的内建数据类型，类似但不完全等同于 NumPy 的ndarray类型：
```
In [4]: z = [-1.0 - 1.0im  1.0 - 1.0im;
             -1.0 + 1.0im  1.0 + 1.0im]
Out[4]: 2x2 Array{Complex{Float64},2}:
 -1.0-1.0im  1.0-1.0im
 -1.0+1.0im  1.0+1.0im
```
我们可以使用方括号[]对数组进行索引。与 Python 的一个显著不同之处在于，索引是从 1 开始的，而不是 0。MATLAB 也有相同的惯例。此外，关键字end表示该维度中的最后一个元素：
```
In [5]: z[1,end]
Out[5]: 1.0 - 1.0im
```

我们可以在矩阵z和标量c（多态）上评估f：

In [6]: f(z, 0)
Out[6]: 2x2 Array{Complex{Float64},2}:
 0.0+2.0im  0.0-2.0im
 0.0-2.0im  0.0+2.0im

现在，我们创建一个函数julia，用于计算 Julia 集。可选的命名参数与位置参数通过分号（;）分隔。Julia 的流程控制语法接近 Python，但省略了冒号，缩进不重要，并且end关键字是必需的：

In [7]: function julia(z, c; maxiter=200)
            for n = 1:maxiter
                if abs2(z) > 4.0
                    return n-1
                end
                z = f(z, c)
            end
            return maxiter
        end
Out[7]: julia (generic function with 1 method)

我们可以在 Julia 中使用 Python 包。首先，我们需要通过 Julia 的内置包管理器（Pkg）安装PyCall包。安装完成后，我们可以通过using PyCall在交互式会话中使用它：
```
In [8]: Pkg.add("PyCall")
        using PyCall
```
我们可以使用@pyimport 宏（Julia 中的元编程特性）导入 Python 包。这个宏相当于 Python 中的import命令：
```
In [9]: @pyimport numpy as np
```

np命名空间现在在 Julia 的交互式会话中可用。NumPy 数组会自动转换为 Julia 的Array对象：

In [10]: z = np.linspace(-1., 1., 100)
Out[10]: 100-element Array{Float64,1}:
 -1.0     
 -0.979798
...
  0.979798
  1.0

我们可以使用列表推导式来在多个参数上评估函数julia：

In [11]: m = [julia(z[i], 0.5) for i=1:100]
Out[11]: 100-element Array{Int64,1}:
 2
...
 2

让我们尝试一下 Gadfly 绘图包。这个库提供了一个高级绘图接口，灵感来源于 Leland Wilkinson 博士的教材*《图形语法》*。在笔记本中，借助D3.js库，图表是交互式的：
```
In [12]: Pkg.add("Gadfly")
         using Gadfly 
In [13]: plot(x=1:100, y=m, Geom.point, Geom.line)
Out[13]: Plot(...)
```
这是一个截图：

在 IPython 笔记本中用 Julia 绘制的 Gadfly 图：

现在，我们通过使用两个嵌套的循环来计算 Julia 集。通常，与 Python 不同，使用for循环而不是向量化操作并不会显著影响性能。高性能的代码可以通过向量化操作或for循环来编写：

In [14]: @time m = [julia(complex(r, i), complex(-0.06, 0.67)) 
                    for i = 1:-.001:-1,
                        r = -1.5:.001:1.5];
elapsed time: 0.881234749 seconds (48040248 bytes allocated)

最后，我们使用PyPlot包来绘制 Julia 中的 matplotlib 图形：

In [15]: Pkg.add("PyPlot")
         using PyPlot
In [16]: imshow(m, cmap="RdGy", 
                extent=[-1.5, 1.5, -1, 1]);

如何实现…

它是如何工作的……

过去，编程语言通常是低级的，使用起来困难但运行快速（例如 C）；或者是高级的，使用方便但运行较慢（例如 Python）。在 Python 中，解决这个问题的方法包括 NumPy 和 Cython 等。

Julia 开发者选择创建一种新的高层次但快速的语言，将两者的优点结合在一起。这本质上通过使用 LLVM 实现的现代即时编译技术来实现。

Julia 动态解析代码并生成低级代码，采用 LLVM 中间表示。这种表示具有独立于语言的指令集，然后编译为机器码。使用显式循环编写的代码会直接编译为机器码。这也解释了为什么在 Julia 中通常不需要进行性能驱动的向量化代码。

还有更多……

Julia 的优势包括：

基于多重分派的强大灵活的动态类型系统，用于参数化多态性
支持元编程的功能
简单的接口，用于从 Julia 调用 C、FORTRAN 或 Python 代码
内建对精细粒度并行和分布式计算的支持
内建的多维数组数据类型和数值计算库
基于 Git 的内建包管理器
用于数据分析的外部包，如 DataFrames（相当于 pandas）和 Gadfly（统计绘图库）
在 IPython notebook 中的集成

Python 相对于 Julia 的优势是什么？截至本文撰写时，Julia 比 Python 和 SciPy 更年轻且不够成熟。因此，Julia 的包和文档比 Python 少。Julia 语言的语法仍在变化。此外，Python 在生产环境中的使用比 Julia 更为普遍。因此，将数值计算代码引入生产环境时，Python 代码会更容易。

话虽如此，Julia 生态系统和社区正在快速增长。我们可以合理预期，Julia 在未来将变得越来越受欢迎。而且，既然两种语言都可以在 IPython notebook 中使用，我们不必在 Python 和 Julia 之间选择。我们可以从 Julia 调用 Python 代码并使用 Python 模块，反之亦然。

我们在这篇食谱中仅仅触及了 Julia 语言的表面。我们未能详细覆盖的一些有趣话题包括 Julia 的类型系统、元编程功能、并行计算支持和包管理器等。

这里有一些参考资料：

Wikipedia 上的 Julia 语言条目可在 en.wikipedia.org/wiki/Julia_%28programming_language%29 查阅
Julia 官方文档可在 docs.julialang.org/en/latest/ 阅读
我们为何创建 Julia 博客文章，可在 julialang.org/blog/2012/02/why-we-created-julia/ 阅读
用于从 Julia 调用 Python 的 PyCall.jl，可在 github.com/stevengj/PyCall.jl 获取
用于在 Julia 中使用 matplotlib 的 PyPlot.jl 可在 github.com/stevengj/PyPlot.jl 获取
Gadfly.jl，一个用于绘图的 Julia 库，访问地址：dcjones.github.io/Gadfly.jl/
DataFrames.jl，Julia 中类似 pandas 的库，访问地址：juliastats.github.io/DataFrames.jl/
Julia Studio，适用于 Julia 的集成开发环境，访问地址：forio.com/labs/julia-studio/

IPython-交互式计算和可视化秘籍第二版-二-