IPython-笔记本精要-全-

IPython 笔记本精要（全）

原文：annas-archive.org/md5/ad987c392fe72c3d7ebd38c7557410b3

译者：飞龙

协议：CC BY-NC-SA 4.0

前言

在过去的 30 年里，计算世界经历了难以置信的革命。就在不久前，高性能计算需要昂贵的硬件；专有软件的费用可能高达数百甚至数千美元；需要掌握 FORTRAN、C 或 C++等计算机语言；并且需要熟悉专门的库。即使获得了适当的硬件和软件，仅仅是为高级科学计算和数据处理建立一个工作环境也是一个巨大的挑战。许多工程师和科学家不得不成为操作系统专家，仅仅为了维护他们日常计算工作所需的工具集。

科学家、工程师和程序员很快就开始解决这个问题。随着性能的提升，硬件成本逐渐降低，并且出现了强烈的推动，开发能够通过多个平台集成不同库的脚本语言。正是在这种环境下，Python 在 1980 年代末期由 Guido Van Rossum 领导下开始开发。从一开始，Python 就被设计成一个前沿的高级计算机语言，结构足够简单，即使是非专家的程序员也能快速学习其基础。

Python 在快速开发中的一个吸引人的特点是它的交互式 Shell，通过这个 Shell，程序员可以在将概念纳入脚本之前进行交互式的实验。然而，原始的 Python Shell 功能有限，更好的交互性是必要的。从 2001 年开始，Fernando Perez 开始开发 IPython，一个专为科学计算设计的改进版交互式 Python Shell。

此后，IPython 已经发展成为一个建立在 Python 之上的完整计算环境。其中最令人兴奋的发展之一就是 IPython 笔记本，它是一个基于网页的 Python 计算界面。在本书中，读者将通过简单的步骤全面了解笔记本的功能。在学习笔记本界面的过程中，读者还将学习多个工具的基本特性，如用于高效数组计算的 NumPy、用于专业级图形绘制的 matplotlib、用于数据处理和分析的 pandas 以及用于科学计算的 SciPy。通过引入与每个主题相关的应用实例，内容变得既有趣又生动。最后但同样重要的是，我们还介绍了使用基于 GPU 的并行计算的高级方法。

我们生活在激动人心的计算时代。廉价但强大的硬件与通过 IPython 笔记本轻松获得的先进库相结合，提供了前所未有的计算能力。我们期待我们的读者能像我们一样，充满动力去探索这个勇敢的新计算世界。

本书内容

第一章，IPython Notebook 概览，展示了如何通过安装 Anaconda 发行版或通过 Wakari 在线连接快速访问 IPython Notebook。你将获得一个介绍性示例，突显 Notebook 界面的一些激动人心的功能。

第二章，Notebook 界面，深入探讨了 Notebook，涵盖了导航、与操作系统交互、运行脚本、加载和保存数据等内容。最后但同样重要的是，我们讨论了 IPython 的丰富显示系统，它允许在 Notebook 中包含各种媒体。

第三章，使用 matplotlib 绘制图形，展示了如何使用 matplotlib 库创建演示质量的图表。阅读本章后，你将能够在 Notebook 中绘制二维和三维数据图，并制作动画。

第四章，使用 pandas 处理数据，展示了如何使用 pandas 库进行数据处理和分析。本章详细研究了该库提供的主要数据结构，并展示了如何访问、插入和修改数据。本章还介绍了数据分析和数据图形显示。

第五章，使用 SciPy、Numba 和 NumbaPro 进行高级计算，介绍了通过 SciPy 访问的高级计算工具和算法。本章还涉及使用 Numba 和 NumbaPro 库的加速技术，包括使用 GPU 进行并行化。

附录 A，IPython Notebook 参考卡，讨论了如何启动 Notebook，编辑模式和命令模式中的快捷键，如何导入模块，以及如何访问各种帮助选项。

附录 B，Python 简要回顾，为读者提供了 Python 语法和特性概览，涵盖了基本类型、表达式、变量和赋值、基本数据结构、函数、对象和方法。

附录 C，NumPy 数组，为我们介绍了 NumPy 数组，并展示了如何创建数组和访问数组成员，最后介绍了索引和切片。

本书所需的内容

为了运行本书中的示例，以下是必需的：

• 操作系统：

  • Windows 7 及以上，32 位或 64 位版本。

  • Mac OS X 10.5 及以上，64 位版本。

  • 基于 Linux 的操作系统，如 Ubuntu 桌面版 14.04 及以上版本，32 位或 64 位版本。

  注意

  请注意，如果可用，推荐使用 64 位版本。

• 软件：

  • Anaconda Python 发行版，版本 3.4 及以上（可从continuum.io/downloads下载）

本书的读者对象

本书面向需要快速了解 IPython 笔记本的科学计算、数据处理与分析、图形显示创建以及高效计算的开发者、工程师、科学家和学生。

假设读者对 Python 编程有一定了解，但 Python 语法的基本内容已在附录中涵盖，所有编程概念也在正文中进行了解释。

如果你在寻找一本节奏适中的 IPython 笔记本入门书籍，其中包含大量应用实例和代码示例，那么本书适合你。

约定

在本书中，你将会看到多种文本风格，用以区分不同类型的信息。以下是这些风格的示例，并附上其含义的解释。

文本中的代码词汇、数据库表名、文件夹名、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 账号都按如下方式显示：“运行 IPython 的最简单方法是，在终端窗口中输入ipython命令。”

代码块如下所示：

temp_coffee = 185.
temp_cream = 40.
vol_coffee = 8.
vol_cream = 1.
initial_temp_mix_at_shop = temp_mixture(temp_coffee, vol_coffee, temp_cream, vol_cream)
temperatures_mix_at_shop = cooling_law(initial_temp_mix_at_shop, times)
temperatures_mix_at_shop

当我们希望特别强调代码块的某一部分时，相关的行或项目会以粗体显示：

[default]
temp_coffee = 185.
temp_cream = 40.
vol_coffee = 8.
vol_cream = 1.
initial_temp_mix_at_shop = temp_mixture(temp_coffee, vol_coffee, temp_cream, vol_cream)
temperatures_mix_at_shop = cooling_law(initial_temp_mix_at_shop, times)
temperatures_mix_at_shop

任何命令行输入或输出都如下所示：

ipython notebook

新术语和重要词汇以粗体显示。你在屏幕上看到的文字，例如在菜单或对话框中的文字，都会像这样出现在文本中：“简单地点击新建笔记本按钮来创建一个新的笔记本。”

注意

警告或重要提示会以框的形式显示，如下所示。

提示

小贴士和技巧如下所示：

读者反馈

我们始终欢迎读者的反馈。请告诉我们你对本书的看法——喜欢的部分或可能不喜欢的部分。读者的反馈对我们非常重要，它帮助我们开发出真正能让你最大程度受益的书籍。

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

如果你在某个领域有专长，并且有兴趣写书或参与编写书籍，请查看我们关于www.packtpub.com/authors的作者指南。

客户支持

现在，你已经成为一本 Packt 书籍的骄傲拥有者，我们提供了许多帮助你充分利用购买内容的资源。

勘误

尽管我们已尽最大努力确保内容的准确性，但错误仍然会发生。如果您在我们的书籍中发现错误——可能是文本或代码中的错误——我们将不胜感激，若您能报告给我们。通过这样做，您可以帮助其他读者避免困扰，并帮助我们改进本书的后续版本。如果您发现任何勘误，请访问www.packtpub.com/submit-errata，选择您的书籍，点击勘误 提交 表格链接，填写您的勘误详情。一旦您的勘误经过验证，您的提交将被接受，勘误将上传至我们的网站，或加入到该书籍的现有勘误列表中，列在该书籍的勘误部分。任何现有的勘误可以通过访问www.packtpub.com/support并选择您的书名查看。

盗版

互联网上的版权盗版问题是各类媒体中一个持续存在的问题。在 Packt，我们非常重视版权和许可的保护。如果您在互联网上发现任何形式的非法复制我们作品的行为，请立即向我们提供该位置地址或网站名称，以便我们采取补救措施。

如您发现涉嫌盗版的资料，请通过<copyright@packtpub.com>联系我们，并附上盗版材料的链接。

我们感谢您在保护作者权益方面的帮助，同时也感谢您帮助我们能够为您带来有价值的内容。

问题

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

第一章：IPython Notebook 概览

本章简要介绍了IPython notebook，并突出其一些特殊功能，这些功能使它成为科学计算和数据导向计算的优秀工具。IPython notebook 使用标准的文本格式，便于共享结果。

在快速安装说明之后，你将学习如何启动 notebook，并能够立即使用 IPython 进行计算。这一简单的初步设置已经足够利用 notebook 的诸多功能，比如交互式生成高质量图形、进行高级技术计算，以及使用专业库处理数据。

本书详细解释了所有示例，并且在线也可以获取。我们不要求读者具备深入的 Python 知识，但对 Python 语法不熟悉的读者可以参考附录 B，Python 简介，以便了解或复习相关内容。

本章将涵盖以下主题：

• 开始使用 Anaconda 或 Wakari

• 创建 notebook，并了解基本的编辑和执行语句

• 一个应用实例，重点展示 notebook 功能

开始使用 Anaconda 或 Wakari

设置 IPython notebook 环境有多种方法。我们建议你使用 Anaconda，这是一个免费的分发版，专为大规模数据处理、预测分析和科学计算设计。或者，你也可以使用 Wakari，它是 Anaconda 的 Web 版安装。Wakari 提供多个服务级别，但基本版是免费的，适合实验和学习。

提示

我们建议你同时设置 Wakari 账户和本地 Anaconda 安装。Wakari 具有便捷的共享和发布功能，而本地安装不需要 Internet 连接，响应速度可能更快。因此，你可以充分利用两者的优点！

安装 Anaconda

要在你的计算机上安装 Anaconda，请执行以下步骤：

1. 从store.continuum.io/cshop/anaconda/下载适合你平台的 Anaconda。

2. 文件完全下载后，安装 Anaconda：

   • Windows 用户可以双击安装程序，并按照屏幕上的指示进行操作。

   • Mac 用户可以双击 .pkg 文件，并按照屏幕上的指示进行操作。

   • Linux 用户可以运行以下命令：

     bash <downloaded file>
     
     

注意

Anaconda 支持多个不同版本的 Python。本书假设你使用的是版本 2.7，这是 Anaconda 默认的标准版本。最新版本的 Python，即版本 3.0，差异较大，且刚刚开始获得广泛使用。许多 Python 包仍然仅在版本 2.7 中得到完全支持。

运行 notebook

现在您已经准备好运行笔记本了。首先，我们创建一个名为my_notebooks的目录来存放您的笔记本，并在该目录中打开一个终端窗口。不同操作系统的操作步骤不同。

微软 Windows 用户需要执行以下步骤：

1. 打开资源管理器。

2. 导航到存储您笔记本的位置。

3. 点击新建文件夹按钮。

4. 重命名文件夹为my_notebooks。

5. 右键点击my_notebooks文件夹，从右键菜单中选择在此处打开命令窗口。

Mac OS X 和其他类 Unix 系统的用户需要执行以下步骤：

1. 打开终端窗口。

2. 运行以下命令：

   mkdir my_notebooks
   cd my_notebooks
   
   

3. 然后，在终端窗口中执行以下命令：

   ipython notebook
   
   

一段时间后，您的浏览器将自动加载笔记本仪表板，如下方的截图所示。仪表板是一个迷你文件系统，您可以在其中管理笔记本。仪表板中列出的笔记本与您在启动笔记本服务器时所在目录中的文件完全对应。

提示

Internet Explorer 不完全支持 IPython 笔记本中的所有功能。建议使用 Chrome、Firefox、Safari、Opera 或其他符合标准的浏览器。如果您的默认浏览器是其中之一，您就可以开始使用了。或者，您也可以关闭 Internet Explorer 笔记本，打开一个兼容的浏览器，并在命令窗口中输入您启动 IPython 时给出的笔记本地址。对于您打开的第一个笔记本，这个地址类似于http://127.0.0.1:8888。

创建 Wakari 账户

要访问 Wakari，只需前往www.wakari.io并创建一个账户。登录后，您将自动进入 Wakari 中使用笔记本界面的介绍页面。此界面在下方的截图中展示：

如上所示的界面元素描述如下：

• 标记为1的部分显示了您的笔记本和文件的目录列表。在该区域的顶部，有一个工具栏，包含创建新文件和目录的按钮，以及下载和上传文件的选项。

• 标记为2的部分显示了欢迎使用 Wakari笔记本。这是包含有关 IPython 和 Wakari 信息的初始笔记本。笔记本界面在第二章，笔记本界面中有详细讨论。

• 标记为3的部分显示了 Wakari 工具栏。此处有新建笔记本按钮以及其他工具的下拉菜单。

注意

本书重点介绍通过笔记本界面使用 IPython。但值得一提的是，还有两种其他的运行 IPython 的方式。

运行 IPython 最简单的方法是在终端窗口中输入ipython命令。这将启动一个 IPython 的命令行会话。正如读者可能猜到的，这是 IPython 的原始界面。

另外，IPython 可以通过 ipython qtconsole 命令启动。这将启动一个附加到 QT 窗口的 IPython 会话。QT 是一种流行的跨平台窗口系统，并随 Anaconda 发行版一起打包。这些替代方法在某些不支持笔记本界面的系统中可能会有用。

创建你的第一个笔记本

我们准备好创建第一个笔记本了！只需点击 新建 笔记本 按钮，即可创建一个新的笔记本。

• 在本地笔记本安装中，新建笔记本按钮出现在仪表盘的左上角。

• 在 Wakari 中，新建笔记本按钮位于仪表盘的顶部，颜色独特。不要使用 添加文件按钮。

请注意，Wakari 仪表盘左侧包含了一个目录列表。你可以利用这个目录以任何方便的方式来组织你的笔记本。Wakari 实际上提供了一个完全可用的 Linux shell。

我们现在可以开始计算了。以下截图展示了笔记本界面：

默认情况下，新建的笔记本命名为 UntitledX，其中 X 是一个数字。要更改名称，只需点击当前标题并编辑弹出的对话框。

在笔记本的顶部，你将看到一个空白框，左侧有 In [ ]: 的文字。这个框被称为代码单元格，是输入 IPython shell 命令的地方。通常，我们在新笔记本中发出的第一个命令是 %pylab inline。请在代码单元格中输入这一行命令，然后按 Shift + Enter（这是执行命令最常用的方式。仅按 Enter 会在当前单元格中创建一行新内容）。执行后，这个命令将显示如下消息：

Populating the interactive namespace from numpy and matplotlib

这个命令使得几种计算工具可以轻松访问，并且是进行交互式计算时推荐的方式。inline 指令告诉 IPython 我们希望图形嵌入到笔记本中，而不是通过外部程序渲染。

以 % 和 %% 开头的命令被称为魔法命令，用于设置配置选项和特殊功能。%pylab 魔法命令会将大量名称导入到 IPython 命名空间中。这个命令通常被认为会引起命名空间污染。推荐的使用库的方法是使用以下命令：

import numpy as np

然后，例如，要访问 NumPy 包中的 arange() 函数，可以使用 np.arange()。这种方法的问题在于，使用常见的数学函数变得繁琐，例如 sin()、cos() 等等。这些函数必须分别输入为 np.sin()、np.cos() 等，导致笔记本的可读性大大降低。

在本书中，我们采用了以下中庸的惯例：在进行交互式计算时，我们将使用%pylab指令，以便更轻松地输入公式。然而，在使用其他库或编写脚本时，我们将采用推荐的最佳实践来导入库。

示例——咖啡冷却问题

假设你在咖啡店买了一杯咖啡。你应该在店里混合奶油，还是等到回到办公室后再混合？我们假设目标是让咖啡尽可能热。因此，主要问题是咖啡在你走回办公室时如何冷却。

混合奶油的两种策略之间的区别是：

• 如果你在店里倒入奶油，在你走回办公室之前，咖啡的温度会突然下降在咖啡开始冷却之前。

• 如果你在回到办公室后再倒奶油，温度骤降发生在步行冷却期之后。

我们需要一个冷却过程的模型。最简单的模型是牛顿冷却定律，它指出冷却速率与咖啡杯中咖啡温度和环境温度之间的温差成正比。这反映了一个直观的概念：例如，如果外部温度是 40°F，咖啡的冷却速度会比外部温度是 50°F 时更快。这个假设得出了一个著名的公式来描述温度变化的方式：

常数r是介于 0 和 1 之间的一个数字，表示咖啡杯与外部环境之间的热交换。这个常数取决于多个因素，可能在没有实验的情况下很难估算。我们在这个初步示例中只是随意选择了一个值。

我们将从设置变量来表示外部温度和冷却速率开始，并定义一个计算液体冷却过程中的温度变化的函数。然后，在一个代码单元中输入表示冷却定律的代码行。按Enter或点击Return来向单元格添加新行。

如前所述，我们首先定义表示外部温度和冷却速率的变量：

temp_out = 70.0
r = 0.9

在单元格中输入上述代码后，按Shift + Enter执行单元格。请注意，在单元格执行后，会创建一个新的单元格。

注释

请注意，尽管此时temp_out的值是一个整数，我们仍然将其输入为70.0。这在此案例中并非严格必要，但它被认为是良好的编程实践。某些代码可能会根据是操作整数变量还是浮动小数变量而表现不同。例如，在 Python 2.7 中，20/8的结果是2，这是 20 除以 8 的整数商。另一方面，20.0/8.0会计算出浮动小数值2.5。通过强制将变量temp_out设为浮动小数值，我们避免了这种意外的行为。

另一个原因是为了提高代码的清晰度和可读性。当阅读者看到70.0这个值时，会很容易理解变量temp_out表示一个实数。因此，可以清楚地知道例如70.8这样的值也会被接受作为外部温度。

接下来，我们定义表示冷却定律的函数：

def cooling_law(temp_start, walk_time):
 return temp_out + (temp_start - temp_out) * r ** walk_time 

注意

请注意缩进的方式，因为 Python 通过缩进来定义代码块。再次按Shift + Enter来执行单元格。

cooling_law()函数接受起始温度和行走时间作为输入，返回咖啡的最终温度。注意，我们只是在定义函数，因此没有输出。在我们的示例中，我们将始终为变量选择有意义的名称。为了保持一致性，我们使用 Google Python 编码风格中的约定，详情见google-styleguide.googlecode.com/svn/trunk/pyguide.html#Python_Language_Rules。

注意

请注意，Python 中的指数运算符是**，而不是其他数学软件中的^。如果在尝试计算幂时遇到以下错误，很可能是你打算使用**运算符：

TypeError: unsupported operand type(s) for ^: 'float' and 'float'

现在我们可以计算在任何起始温度和行走时间下冷却的效果。例如，要计算咖啡在 10 分钟后的温度，假设初始温度为 185°F，请在单元格中运行以下代码：

cooling_law(185.0, 10.0)

相应的输出是：

110.0980206115

如果我们想知道多个行走时间下的最终温度该怎么办？例如，假设我们想每 5 分钟计算一次直到 20 分钟的最终温度。这正是NumPy让事情变得简单的地方：

times = arange(0.,21.,5.)
temperatures = cooling_law(185., times)
temperatures

我们首先通过使用arange()函数将times定义为一个NumPy数组。这个函数有三个参数：范围的起始值、范围的结束值和增量。

注意

你可能会想，为什么范围的结束值是21而不是20。这是计算机科学中常见的约定，Python 也遵循这一约定。当指定一个范围时，右端点永远不属于该范围。因此，如果我们将20指定为右端点，则该范围只包含0、5、10和15这几个值。

定义了times数组后，我们可以简单地用times作为第二个参数调用cooling_law()函数。这将计算在给定时间点的温度。

注意

你可能注意到这里有些奇怪的地方。第一次调用cooling_law()函数时，第二个参数是一个浮点数。第二次调用时，它是一个NumPy数组。这是得益于 Python 的动态类型和多态性。NumPy重新定义了算术运算符，以智能的方式处理数组。因此，我们不需要为这种情况特别定义新的函数。

一旦我们得到了温度，就可以将它们展示在图表中。要展示图表，请在单元格中执行以下命令：

plot(times, temperatures, 'o')

上述命令行会生成以下图表：

plot()函数是matplotlib包的一部分，我们将在第三章中详细学习，使用 matplotlib 绘图。在这个例子中，plot()的前两个参数是NumPy数组，分别指定水平和垂直坐标轴的数据。第三个参数指定绘图符号为实心圆。

我们现在准备好解决最初的问题：我们是应该在咖啡店里混合奶油，还是等回到办公室再混合？当我们混合奶油时，温度会突然下降。混合液体的温度是两种液体温度的加权平均值，权重由体积决定。以下代码定义了一个函数，用于计算混合后的温度：

def temp_mixture(t1, v1, t2, v2):
 return (t1 * v1 + t2 * v2) / (v1 + v2)

函数中的参数是每种液体的温度和体积。使用这个函数，我们现在可以计算在咖啡店混合奶油时的温度变化：

temp_coffee = 185.
temp_cream = 40.
vol_coffee = 8.
vol_cream = 1.
initial_temp_mix_at_shop = temp_mixture(temp_coffee, vol_coffee, temp_cream, vol_cream)
temperatures_mix_at_shop = cooling_law(initial_temp_mix_at_shop, times)
temperatures_mix_at_shop

提示

请注意，我们在单元格的末尾重复了变量temperatures_mix_at_shop。这并不是一个拼写错误。默认情况下，IPython 笔记本假设单元格的输出是单元格中计算的最后一个表达式。通常的做法是在单元格末尾列出希望显示的变量。稍后我们将看到如何显示更精美、格式化的输出。

像往常一样，把所有命令输入一个代码单元格中，然后按Shift + Enter运行整个单元格。我们首先设置咖啡和奶油的初始温度和体积。然后，我们调用temp_mixture()函数来计算混合物的初始温度。最后，我们使用cooling_law()函数来计算不同步行时间下的温度，将结果存储在temperatures_mix_at_shop变量中。上述命令行会生成以下输出：

array([ 168.88888889,  128.3929    ,  104.48042352,   90.36034528,
 82.02258029])

记住，times数组指定了从0到20的时间，间隔为 5 分钟。因此，上述输出给出了这些时刻的温度，假设奶油是在咖啡店混合的。

要计算考虑到奶油是在回到办公室后混合时的温度，请在单元格中执行以下命令：

temperatures_unmixed_coffee = cooling_law(temp_coffee, times)
temperatures_mix_at_office = temp_mixture(temperatures_unmixed_coffee, vol_coffee, temp_cream, vol_cream)
temperatures_mix_at_office

我们再次使用cooling_law()函数，但将未混合奶油的初始咖啡温度temp_coffee作为第一个输入变量。我们将结果存储在temperatures_unmixed_coffee变量中。

为了计算走路后混合奶油的效果，我们调用了temp_mixture()函数。请注意，两次计算的主要区别在于cooling_law()和temp_mixture()函数调用的顺序。上述命令行会生成以下输出：

array([ 168.88888889,  127.02786667,  102.30935165,   87.71331573,
 79.09450247])

现在让我们绘制两个温度数组。请在一个单元中执行以下命令行：

plot(times, temperatures_mix_at_shop, 'o')
plot(times, temperatures_mix_at_office, 'D', color='r')

第一次调用 plot() 函数与之前相同。第二次调用类似，但我们希望绘制的符号是填充的菱形，用参数 'D' 来表示。color='r' 选项使标记为红色。这将生成以下图表：

请注意，默认情况下，在一个代码单元中创建的所有图表都会绘制在同一坐标轴上。总的来说，我们可以得出结论，对于本示例中使用的数据参数，不管步行时间如何，在咖啡店搅拌奶油总是更好。读者可以自由更改参数并观察不同情况下会发生什么。

科学绘图应当清晰地表示所展示的内容、绘制的变量，以及所使用的单位。这可以通过向图表添加注释来很好地处理。如下代码所示，matplotlib 中添加注释是非常简单的：

plot(times, temperatures_mix_at_shop, 'o')
plot(times, temperatures_mix_at_office, 'D', color='r')
title('Coffee temperatures for different walking times')
xlabel('Waking time (min)')
ylabel('Temperature (F)')
legend(['Mix at shop', 'Mix at office']) 

在再次绘制数组之后，我们调用适当的函数来添加标题（title()）、横轴标签（xlabel()）、纵轴标签（ylabel()）和图例（legend()）。所有这些函数的参数都是字符串或字符串列表，例如 legend() 函数的情况。以下图表是我们根据前面的命令行得到的输出：

我们在进行此分析时有些不满意的地方；我们的办公室据称与咖啡店有固定的距离。此情况中的主要因素是外部温度。在夏季和冬季是否应该使用不同的策略？为了探讨这个问题，我们首先定义一个函数，接受奶油温度和外部温度作为输入。该函数的返回值是我们回到办公室时的最终温度差。

该函数定义如下：

temp_coffee = 185.
vol_coffee = 8.
vol_cream = 1.
walk_time = 10.0
r = 0.9
def temperature_difference(temp_cream, temp_out):
 temp_start = temp_mixture(temp_coffee,vol_coffee,temp_cream,vol_cream)
 temp_mix_at_shop = temp_out + (temp_start - temp_out) * r ** walk_time
 temp_start = temp_coffee
 temp_unmixed =  temp_out + (temp_start - temp_out) * r ** walk_time
 temp_mix_at_office = temp_mixture(temp_unmixed, vol_coffee, temp_cream, vol_cream)
 return temp_mix_at_shop - temp_mix_at_office

在前面的函数中，我们首先设置了在分析中将被视为常数的变量值，即咖啡的温度、咖啡和奶油的体积、步行时间和降温速率。然后，我们使用之前讨论过的相同公式定义了 temperature_difference 函数。现在，我们可以使用这个函数计算一个矩阵，表示不同奶油温度和外部温度下的温差：

cream_temperatures = arange(40.,51.,1.)
outside_temperatures = arange(35.,60.,1.)
cream_values, outside_values = meshgrid(cream_temperatures, outside_temperatures)
temperature_differences = temperature_difference(cream_values, outside_values)

单元格中的前两行使用arange()函数定义数组，这些数组包含等间距的奶油温度和外部温度值。然后，我们调用便捷函数meshgrid()。该函数返回两个数组，便于计算三维图表的数据，并将它们存储在变量cream_values和outside_values中。最后，我们调用temperature_difference()函数，并将返回值存储在temperature_differences数组中。这将是一个二维数组（即矩阵）。

现在我们准备绘制温差的三维图：

from mpl_toolkits.mplot3d import Axes3D
fig = figure()
fig.set_size_inches(7,5)
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(cream_values, outside_values, temperature_differences, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0)
xlabel('Cream temperature')
ylabel('Outside temperature')
title('Temperature difference') 

在前面的代码段中，我们首先通过以下代码行导入Axes3D类：

from mpl_toolkits.mplot3d import Axes3D

这个类位于mpl_toolkits.mplot3d模块中，并不会自动加载。因此，必须显式导入它。

然后，我们创建一个fig对象，它是figure类的一个实例，设置它的大小，并生成一个ax对象，它是Axes3D类的一个实例。最后，我们调用ax.plot_surface()方法生成图表。最后三行命令设置坐标轴标签和标题。

注意

在这个解释中，我们使用了一些面向对象编程中常见的术语。Python 中的对象仅仅是可以以某种特定方式处理的数据结构。每个对象都是类的实例，该类定义了对象的数据。类还定义了方法，这些方法是专门用于操作属于该类的对象的函数。

前述命令行的输出将生成以下图表：

请注意在调用ax.plot_surface()时的cmap=cm.coolwarm参数。它将图表的颜色映射设置为cm.coolwarm。这个颜色映射便捷地使用了蓝色到红色的渐变来表示函数值。因此，负温差用蓝色显示，正温差用红色显示。请注意，似乎有一条直线定义了温差从负到正的过渡位置。实际上，这对应的是外部温度与奶油温度相等的值。结果表明，如果奶油温度低于外部温度，我们应该在咖啡店将奶油混入咖啡中；否则，奶油应当在办公室倒入咖啡中。

练习

以下是一些练习题，帮助你理解并应用本章所学的概念：

• 在我们的例子中，我们讨论了如何确定冷却速率r。请修改该示例，绘制不同r值下的温度变化图，保持其他变量不变。

• 查阅matplotlib文档，网址为matplotlib.org，以了解如何生成温差的等高线图。

  我们的分析忽略了一个事实，即在行走过程中，奶油的温度也会发生变化。请修改笔记本，以便考虑到这一因素。

总结

在本章中，我们通过 Anaconda 设置了一个 IPython 环境，通过 Wakari 在线访问了 IPython notebook，创建了一个笔记本，并学习了编辑和执行命令的基础知识，最后，我们通过一个广泛应用的示例展示了笔记本的基本功能。

在下一章中，我们将更深入地探讨笔记本界面提供的功能——包括笔记本导航和编辑功能、与操作系统的交互、加载和保存数据以及运行脚本。

第二章。笔记本界面

IPython 笔记本具有广泛的用户界面，适合创建丰富格式的文档。本章将全面探索笔记本的功能。我们还将考虑使用笔记本时的陷阱和最佳实践。

本章将涵盖以下主题：

• 笔记本的编辑与浏览，包括单元格类型；添加、删除和移动单元格；加载和保存笔记本；以及键盘快捷键

• IPython 魔法命令

• 与操作系统交互

• 运行脚本、加载数据和保存数据

• 使用 IPython 的丰富显示系统嵌入图片、视频和其他媒体

编辑和浏览笔记本

当我们打开一个笔记本（通过点击仪表板中的名称或创建一个新笔记本）时，浏览器窗口中将显示如下内容：

在上面的截图中，从上到下，我们看到以下组件：

• 标题栏（标记为1的区域），其中包含笔记本的名称（在上面的例子中，我们可以看到第二章）以及有关笔记本版本的信息

• 菜单栏（标记为2的区域）看起来像一个常规的应用程序菜单

• 工具栏（标记为3的区域）用于快速访问最常用的功能

• 在标记为4的区域，显示一个空的计算单元格

从 IPython 版本 2.0 开始，笔记本有两种操作模式：

• 编辑：在此模式下，单个单元格被聚焦，我们可以在该单元格中输入文本、执行代码并执行与该单元格相关的任务。通过点击单元格或按Enter键激活编辑模式。

• 命令：在此模式下，我们执行与整个笔记本结构相关的任务，如移动、复制、剪切和粘贴单元格。可以使用一系列键盘快捷键来提高这些操作的效率。通过点击笔记本中的任何地方，单元格之外，或按Esc键来激活命令模式。

当我们打开一个笔记本时，它处于命令模式。让我们进入新笔记本的编辑模式。为此，可以点击空白单元格或按Enter键。笔记本的外观会发生轻微变化，如下图所示：

请注意所选单元格周围的粗边框以及笔记本菜单右上角的小铅笔图标。这些表示笔记本处于编辑模式。

在接下来的子章节中，我们将详细探讨每种笔记本模式。

获取帮助和中断计算

笔记本是一个复杂的工具，集成了多种不同的技术。新用户（甚至是有经验的用户）不太可能记住所有命令和快捷键。笔记本中的帮助菜单提供了相关文档的链接，应该根据需要经常参考。

注意

新手可能想访问笔记本界面导览，可以在nbviewer.ipython.org/github/ipython/ipython/blob/2.x/examples/Notebook/User%20Interface.ipynb上找到，作为入门帮助。

同样，获取任何对象（包括函数和方法）的帮助也很容易。例如，要访问sum()函数的帮助，可以在单元格中运行以下代码：

sum?

在对象名称后附加??将提供更详细的信息。顺便提一下，只在单元格中运行?会显示有关 IPython 特性的相关信息。

另一个需要从一开始就了解的重要事项是如何中断计算。这可以通过内核菜单完成，在这里可以中断并重新启动运行笔记本代码的内核进程。也可以通过点击工具栏上的停止按钮来中断内核。

编辑模式

编辑模式用于在单元格中输入文本并执行代码。让我们在新创建的笔记本中输入一些代码。像往常一样，我们希望将NumPy和matplotlib导入到当前命名空间，因此我们在第一个单元格中输入以下魔法命令：

%pylab inline

按Shift + Enter 或点击工具栏上的播放按钮来执行代码。注意，无论选择哪种方式，都会在当前单元格下方添加一个新单元格。

为了有一些具体的内容可以操作，假设我们想计算一项投资所积累的利息。请在连续三个单元格中输入以下代码：

• 在单元格 1 中，输入以下命令行：

  def return_on_investment(principal, interest_rate, number_of_years):
   return principal * e ** (interest_rate * number_of_years)
  
  

• 在单元格 2 中，输入以下命令行：

  principal = 250
  interest_rate = .034
  tstart = 0.0
  tend = 5.0
  npoints = 6
  
  

• 在单元格 3 中，输入以下命令行：

  tvalues = linspace(tstart, tend, npoints)
  amount_values = return_on_investment(principal, interest_rate, tvalues)
  plot(tvalues, amount_values, 'o')
  title('Return on investment, years {} to {}'.format(tstart, tend))
  xlabel('Years')
  ylabel('Return')
  tstart += tend
  tend += tend
  
  

现在，执行以下步骤：

1. 通过按Shift + Enter 以常规方式运行单元格 1 和单元格 2。

2. 而是通过按Ctrl + Enter 来运行单元格 3。

请注意，在执行后，单元格 3 继续被选中。保持按住Ctrl + Enter，并确保选中单元格 3。每次都会更新图表，以显示不同 5 年期的投资回报。

代码的工作原理如下：

• 在单元格 1 中，我们定义了一个函数，用于计算给定本金、利率和年数的投资回报。

• 在单元格 2 中，我们设置了本金和利息的实际值，并初始化了变量，以定义我们希望进行计算的时间段。

• 单元格 3 计算了 5 年期的回报金额并绘制了结果。

• 然后，变量tstart和tend被更新。命令行如下：

  tstart += tend
  tend += tend
  
  

其效果是，下次更新单元格时，时间会推进到下一个 5 年期。因此，通过反复按Ctrl + Enter，我们可以快速查看投资在连续 5 年期中的增长情况。

还有第三种方法可以运行单元格中的命令。再次点击单元格 2，然后在 Windows 中按下Alt + Enter，或者在 Mac 上按Option + Enter。这将运行单元格 2，并在其下插入一个新单元格。暂时忽略这个新单元格。我们不需要这个单元格，下一节中我们将学习如何删除它。

所以，运行单元格内容有三种方法：

• 按下Shift + Enter或工具栏上的播放按钮。这将运行单元格并选择下一个单元格（如果在笔记本的末尾，则创建一个新单元格）。这是执行单元格的最常用方法。

• 按下Ctrl + Enter。这将运行单元格并保持选择当前单元格。我们想要反复执行同一个单元格时这非常有用。例如，如果我们想对现有代码进行修改。

• 按下Alt + Enter。这将运行单元格并立即在其下插入一个新单元格。

编辑模式的另一个有用特点是Tab 自动补全。选择一个空单元格并输入以下命令：

print am

然后，按下Tab键。将会出现一个建议补全列表。使用键盘的箭头键或鼠标，我们可以选择amount_values，然后按Enter接受补全。

IPython 的一个非常重要的特点是可以轻松访问帮助信息。点击一个空单元格并输入：

linspace

然后，按Shift + Tab。将会出现一个包含linspace函数信息的工具提示。点击工具提示窗口右上角的+符号，可以获得更多信息。点击^符号，信息会显示在笔记本底部的信息区域。

提示

Tab和Shift + Tab功能是笔记本中最有用的功能；一定要经常使用它们！

命令模式

在命令模式中可用的快捷键远远多于编辑模式中的快捷键。幸运的是，并不需要一次性记住所有快捷键，因为命令模式中的大部分操作在菜单中也有提供。在本节中，我们只会描述一些命令模式的常见功能。以下表格列出了编辑单元格的一些有用快捷键；其他快捷键将稍后介绍：

快捷键	动作
Enter	激活编辑模式
Esc	激活命令模式
H	显示键盘快捷键列表
S 或 Ctrl + S	保存笔记本
A	在当前单元格上方插入一个单元格
B	在当前单元格下方插入一个单元格
D（按两次）	删除单元格
Z	撤销上一次删除
Ctrl + K	将单元格向上移动
Ctrl + J	将单元格向下移动
X	剪切单元格内容
C	复制单元格的内容
V	将当前单元格下方的内容粘贴到此单元格
Shift + V	将当前单元格上方的内容粘贴到此单元格

注意

使用笔记本时，最常见（也是最让人沮丧）的错误之一是输入错误的模式。记得使用Esc切换到命令模式，使用Enter切换到编辑模式。同时，记得点击单元格会自动进入编辑模式，所以需要按Esc回到命令模式。

继续尝试示例笔记本中的一些编辑快捷键。以下是一个你可以尝试的示例：

1. 按Esc进入命令模式。

2. 使用箭头键移动到我们在上一节中在单元格 2 和单元格 3 之间创建的空单元。

3. 按D两次，这会删除该单元。要恢复单元，按Z。

注意

注意，一些快捷键与其他软件中的常规快捷键不同。例如，剪切、复制和粘贴单元格的快捷键并不以Ctrl键为前缀。

单元类型

到目前为止，我们只使用笔记本单元输入代码。然而，我们也可以使用单元格输入解释性文本并为笔记本提供结构。笔记本使用Markdown语言来轻松地在单元格中插入丰富文本。Markdown 是 John Gruber 为 HTML 的纯文本编辑而创建的。请参阅项目页面daringfireball.net/projects/markdown/basics，了解语法的基本知识。

让我们在笔记本中看看它是如何工作的。如果你在上一节中创建了其他单元来尝试快捷键，现在删除它们，以便笔记本中只剩下%pylab inline单元和三个进行兴趣计算的单元。

点击%pylab inline单元并在其下方插入一个新单元。你可以使用菜单，或者进入命令模式（使用Esc键），然后使用快捷键 B。

我们现在想将新的单元格类型转换为 Markdown。有三种方法可以做到这一点。首先点击选中单元格，然后执行以下步骤之一：

• 点击笔记本菜单项单元格，选择单元格类型，然后点击Markdown，如下图所示。

• 从笔记本工具栏的下拉框中选择Markdown。

• 按Esc进入命令模式，然后按M。

注意，一旦单元转换为 Markdown，它会自动进入编辑模式。现在，在新的 Markdown 单元中输入以下内容（小心在指示的地方留出额外的空行）：

We want to know how an investment grows with a fixed interest.

The *compound interest* formula states that:
$$R = Pe^{rt}$$
where:

- $P$ is the principal (initial investment).
- $r$ is the annual interest rate, as a decimal.
- $t$ is the time in years.
- $e$ is the base of natural logarithms.
- $R$ is the total return after $t$ years (including principal)

For details, see the [corresponding Wikipedia entry](http://en.wikipedia.org/wiki/Compound_interest).

We start by defining a Python function that implements this formula.

输入文本后，按Shift + Enter执行该单元。笔记本将通过 Markdown 解释器而不是 IPython 解释器来运行该单元，单元会通过 HTML 呈现，产生以下截图中显示的输出：

在这个示例中，我们使用了以下 Markdown 功能：

• 文本通常输入，并且通过在文本中留出额外的空行来指示新段落。

• 斜体通过将文本包裹在星号之间表示，如 *复利*。

• 用双美元符号（$$）包围的公式，如 $$R = Pe^{rt}$$，会居中显示在页面上。

• 无序列表通过以破折号（-）开头的行表示。重要的是，在列表前后留出空行。

• 单个美元符号（$）会导致公式内联排版。

• 超链接的格式如下：对应的维基百科条目。

在 Markdown 单元格中，数学公式可以使用LaTeX输入，LaTeX 是一种广泛用于技术排版的语言，超出了本书的范围。幸运的是，我们不需要使用 LaTeX 的完整排版功能，而只需使用其公式编辑功能。可以在en.wikibooks.org/wiki/LaTeX/Mathematics找到 LaTeX 的简明入门教程。学习一些 LaTeX 是非常有用的，因为它也被用于其他 Python 库。例如，matplotlib 允许在图表标题和坐标轴标签中使用 LaTeX。在笔记本中，LaTeX 由 MathJax 渲染，MathJax 是一个由 Davide Cervone 用 JavaScript 实现的 LaTeX 解释器。详细信息请访问 www.mathjax.org/。

注意

要编辑已经显示的 Markdown 单元格的内容，只需双击该单元格。编辑完成后，使用 Shift + Enter 运行该单元格以重新渲染它。

为了给笔记本添加结构，我们可以添加不同大小的标题。让我们为我们的笔记本添加一个全局标题：

1. 在笔记本的最顶部添加一个新单元格，并将其类型更改为 标题 1。回想一下，有三种方法可以做到这一点：

   • 通过导航到 单元格 | 单元格类型

   • 使用工具栏上的 单元格类型 下拉菜单

   • 在命令模式下使用快捷键 1

2. 为笔记本输入一个标题，并使用 Shift + Enter 运行该单元格。

该笔记本支持六种标题大小，从标题 1（最大）到标题 6（最小）。

注意

Markdown 语言还允许插入标题，使用井号（#）符号。尽管这种方式可以节省输入，但我们建议使用标题 1 到标题 6 单元格。将标题放在单独的单元格中可以保持笔记本的结构，尤其是在保存时。这个结构会被 nbconvert 工具使用。

以下表格总结了我们目前考虑过的单元格类型：

单元格类型	命令模式快捷键	用途
代码	Y	这允许你编辑并写入新代码到 IPython 解释器。默认语言是 Python。
Markdown	M	这允许你编写解释性文本。
标题 1 到 标题 6	键位 1 到 6	这允许你对文档进行结构化
Raw NBConvert	R	当笔记本转换为其他格式时，该单元格的内容保持不变

IPython 魔法命令

魔法是 IPython 解释器的特殊指令，用于执行特定的操作。魔法有两种类型：

• 面向行的：这种类型的魔法以一个百分号（%）开头。

• 面向单元格的：这种类型的魔法以两个百分号（%%）开头。

我们已经熟悉其中一个魔法命令，即 %pylab inline。这个特定的魔法执行了以下两件事：它导入了 NumPy 和 matplotlib，并将笔记本设置为内联图表。要查看其他选项，将单元格更改为 %pylab。

运行这个单元格，然后再运行生成图表的单元格。IPython 不再直接在内联中绘制图形，而是会打开一个新窗口，显示如下截图中的图表：

这个窗口是交互式的，你可以调整图表的大小、移动它，并从这里保存到文件。

另一个有用的魔法是 %timeit，它记录运行一行 Python 代码所花费的时间。在一个新的单元格中运行以下代码：

%timeit return_on_investment(principal, interest_rate, tvalues)

输出将是如下所示：

100000 loops, best of 3: 3.73 µs per loop

为了获得更精确的估算，命令会运行 10,000 次，并计算平均运行时间。这会进行三次，并报告最好的结果。

%timeit 魔法在编辑模式下也可以使用。为了演示这一点，运行以下命令：

principal = 250
interest_rates = [0.0001 * i for i in range(100000)]
tfinal = 10

在下一个单元格中，运行以下命令：

%%timeit
returns = []
for r in interest_rates:
 returns.append(return_on_investment(principal, r, tfinal))

上述代码计算了 100,000 个不同利率值的返回值列表，但仅使用了普通的 Python 代码。该代码的运行时间将在以下输出中显示：

10 loops, best of 3: 31.6 ms per loop

现在，让我们使用 NumPy 数组重写相同的计算。运行以下命令：

principal = 250
interest_rates = arange(0, 10, 0.0001)
tfinal = 10

在下一个单元格中，运行以下命令：

%%timeit
returns = return_on_investment(principal, interest_rates, tfinal)

现在，运行时将在以下输出中显示：

100 loops, best of 3: 5.53 ms per loop

比较这两个输出后，我们可以看到使用 NumPy 获得的速度提升。

要列出所有可用的魔法，运行以下命令：

%lsmagic

一旦你拥有了所有魔法的列表，你可以通过运行一个带有魔法名称后加问号的单元格来查询特定的魔法：%pylab?。

这将把有关 %pylab 魔法的信息显示在浏览器底部的一个单独窗口中。

另一个有趣的功能是能够运行用其他语言编写的代码。为了说明这些可能性，我们将展示如何使用 Cython 加速代码，因为 Cython 可以将 Python 代码编译成 C 代码。让我们编写一个函数，计算由正弦曲线所围成的区域的近似值。下面是如何在纯 Python 中定义这个函数：

import math
def sin_area(a, b, nintervals):
 dx = (b-a)/nintervals
 sleft = 0.0
 sright = 0.0
 for i in range(nintervals):
 sleft += math.sin(a + i * dx)
 sright += math.sin(a + (i + 1) * dx)
 return dx * (sright + sleft) / 2 

我们将通过取左端点法和右端点法的平均值来近似计算面积（这等同于梯形法则）。代码确实效率低下且不符合 Python 风格。特别要注意的是，我们使用的是 Python 库版本的 sin() 函数，而不是 NumPy 的实现。在这个例子中，NumPy 的实现实际上会导致更慢的代码，因为在列表和数组之间的转换频繁发生。

要运行一个简单的测试，在单元格中执行以下命令：

sin_area(0, pi, 10000)

运行前面的单元格后，我们得到了以下输出：

1.9999999835506606

输出是合理的，因为面积的实际值是 2。现在，我们来使用以下命令计算执行时间：

%timeit sin_area(0, pi, 10000)

我们将得到以下输出：

100 loops, best of 3: 3.7 ms per loop

现在我们来用 Cython 实现相同的功能。由于 Cython 魔法是通过扩展模块实现的，我们需要先加载这个模块。我们将使用以下命令加载扩展模块：

%load_ext cythonmagic

现在，我们将定义 Cython 函数。我们不会详细讨论语法，但请注意，它与 Python 非常相似（在这个例子中，主要的区别是我们必须声明变量以指定它们的 C 类型）：

%%cython
cimport cython
from libc.math cimport sin

@cython.cdivision(True)
def sin_area_cython(a, b, nintervals):
 cdef double dx, sleft, sright
 cdef int i
 dx = (b-a)/nintervals
 sleft = 0.0
 sright = 0.0
 for i in range(nintervals):
 sleft += sin(a + i * dx)
 sright += sin(a + (i + 1) * dx)
 return dx * (sright + sleft) / 2

使用以下命令测试前面的函数：

sin_area_cython(0, pi, 10000)

运行前面的函数后，我们得到的输出与之前相同：

1.9999999835506608

现在，我们来使用以下命令计算该函数的执行时间：

%timeit sin_area_cython(0, pi, 10000)

运行时显示在以下输出中：

1000 loops, best of 3: 1.12 ms per loop

我们看到 Cython 代码的执行时间大约是 Python 代码总时间的 30%。需要强调的是，这不是加速这段代码的推荐方法。一种更简单的解决方案是使用 NumPy 来矢量化计算：

def sin_area_numpy(a, b, nintervals):
 dx = (b - a) / nintervals
 xvalues = arange(a, b, dx)
 sleft = sum(sin(xvalues))
 sright = sum(sin(xvalues + dx))
 return dx * (sleft + sright) / 2

运行前面的代码后，以下输出显示了时间：

1000 loops, best of 3: 248 µs per loop

这里有一个经验教训：当我们试图加速代码时，首先要尝试的是总是使用 NumPy 数组编写代码，利用矢量化函数的优势。如果需要进一步加速，我们可以使用像 Numba 和 NumbaPro（本书后面会讨论）这样的专用库来加速代码。事实上，这些库提供了一种比直接使用 Cython 更简单的将代码编译为 C 的方法。

与操作系统交互

任何具有一定复杂性的代码，在必须打开和关闭文件、运行脚本或访问在线数据时，都会与计算机的操作系统交互。普通的 Python 已经有很多工具可以访问这些功能，而 IPython 和笔记本则提供了更多的功能和便利性。

保存笔记本

笔记本会定期自动保存。默认的间隔时间是 2 分钟，但可以在配置文件中修改或使用 %autosave 魔法进行更改。例如，要将自动保存间隔改为 5 分钟，可以运行以下命令：

%autosave 300

请注意，时间单位为秒。要禁用自动保存功能，运行以下命令：

%autosave 0

我们还可以使用File菜单或点击工具栏上的磁盘图标来保存笔记本。这将创建一个检查点。检查点存储在一个隐藏文件夹中，可以通过File菜单恢复。请注意，只有最新的检查点会被提供。

笔记本以.ipynb扩展名保存为纯文本文件，使用 JSON 格式。JSON 是一种广泛用于 Web 应用程序信息交换的格式，相关文档可以参考www.json.org/。这使得与他人交换笔记本变得非常简单：只需将.ipynb文件给他们，他们就可以将其复制到适当的工作目录中。下次在该目录中打开笔记本服务器时，新笔记本将可用（或者可以从仪表板刷新目录列表）。此外，由于 JSON 是纯文本格式，它也可以轻松进行版本管理。

将笔记本转换为其他格式

笔记本可以使用nbconvert工具转换为其他格式。该工具是一个命令行工具。因此，要使用它，请在包含笔记本文件的目录中打开终端窗口。

提示

Windows 用户可以按Shift键并右键单击包含笔记本文件的目录名称，然后选择在此处打开命令窗口。

打开一个 shell 窗口并输入以下命令：

ipython nbconvert "Chapter 2.ipynb"

当然，你必须将Chapter 2.ipynb替换为包含你的笔记本的文件名。文件名必须用引号括起来。

上述命令将把笔记本转换为静态 HTML 页面，可以直接放置在 Web 服务器上。

注意

另一种在 Web 上发布笔记本的方法是使用网站nbviewer.ipython.org/。

还可以使用以下命令，通过nbconvert创建幻灯片放映：

ipython nbconvert "Chapter 2.ipynb" --to slides

然而，为了获得良好的展示效果，我们必须首先在笔记本中指定其结构。为此，请进入笔记本并在Cell工具栏下拉列表中选择幻灯片放映。然后，确定每个单元格是否应该是幻灯片、子幻灯片或片段。

要查看幻灯片放映，需要将reveal.js文件安装到与包含展示文档的网页相同的目录中。你可以从github.com/hakimel/reveal.js下载该文件。如有必要，将包含所有文件的目录重命名为reveal.js。然后，我们就可以打开包含展示内容的 HTML 文件。

还可以将笔记本转换为 LaTeX 和 PDF 格式。但是，这需要安装一些 Anaconda 中没有包含的包。

运行 shell 命令

我们可以通过在单元格前加上感叹号（!）来运行任何 shell 命令。例如，要在 Windows 中获取目录列表，可以在单元格中运行以下命令：

!dir

在 Linux 或 OS X 中的等效命令如下：

!ls

你可以在单元格中输入任何复杂的命令行。例如，以下命令将编译每个 C 语言学生都要尝试的经典"Hello, world!"程序：

!cc hello.c –o hello

当然，除非你安装了 C 编译器cc并且有适当代码的hello.c文件，否则在你的电脑上运行时不会正常工作。

代替直接使用 Shell 命令，许多相同的功能可以通过魔法命令实现。例如，获取目录列表（在任何操作系统中）可以通过运行以下命令：

%ls

以下表格展示了一些常用的魔法命令，用于与系统交互：

Magic	用途
%cd	更改目录
%pwd	打印当前目录
%ls	列出当前目录的内容
%mkdir	创建一个新目录
%rmdir	删除目录
%echo	打印字符串
%alias	创建别名

%echo魔法命令常用于打印环境变量的值。例如，要打印 Windows 中PATH环境变量的内容，可以运行以下命令：

%echo %PATH%

在 Linux 或 OS X 中，使用以下命令：

%echo $PATH

%alias魔法命令可以为常用命令创建别名。例如，要定义一个显示 Windows 中PATH值的宏，可以执行以下命令：

%alias show_path echo %PATH%

在 Linux 或 OS X 中，使用以下命令：

%alias show_path echo $PATH

在定义了前面的命令后，我们可以运行以下命令来显示路径：

show_path

为了让输入命令更加便捷，名为automagic的功能允许直接输入面向行的魔法命令，而无需使用%符号（如前面所示的命令）。例如，要创建一个目录，我们可以直接输入以下命令：

mkdir my-directory

如果我们改变主意，可以使用以下命令删除该目录：

rmdir my-directory

automagic 功能由%automagic魔法控制。例如，使用以下命令关闭 automagic：

%automagic off

运行脚本、加载数据和保存数据

在处理一些复杂项目时，通常需要运行其他人编写的脚本。加载数据和保存结果也是必不可少的。在本节中，我们将描述 IPython 为这些任务提供的功能。

运行 Python 脚本

以下 Python 脚本生成一个洛伦兹方程解的图形，这是混沌理论中的一个著名例子。如果你正在输入代码，请不要在笔记本的单元格中输入，而是使用文本编辑器，将文件保存为lorenz.py，并放在包含笔记本文件的同一目录中。代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D

def make_lorenz(sigma, r, b):
    def func(statevec, t):
        x, y, z = statevec
        return [ sigma * (y - x),
                 r * x - y - x * z,
                 x * y - b * z ]
    return func

lorenz_eq = make_lorenz(10., 28., 8./3.)

tmax = 50
tdelta = 0.005
tvalues = np.arange(0, tmax, tdelta) 
ic = np.array([0.0, 1.0, 0.0])
sol = odeint(lorenz_eq, ic, tvalues)

x, y, z = np.array(zip(*sol))

fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z, lw=1, color='red')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')
plt.show()

现在，转到笔记本并使用以下命令运行一个单元格：

%run lorenz.py

这将运行脚本并生成解的图形，如下图所示：

%run魔法命令在笔记本的命名空间中执行脚本，因此脚本中定义的所有全局变量、函数和类都会在当前笔记本中可用。

也可以使用 %load 魔法来实现相同的目的：

%load lorenz.py

不同之处在于，%load 不会立即运行脚本，而是将其代码放入一个单元格中。然后可以从插入的单元格中运行该代码。%load 魔法的一个稍显恼人的行为是，即使之前已经使用 %load 插入了一个单元格，它仍会插入一个包含脚本代码的新单元格。笔记本无法知道用户是否希望覆盖现有单元格中的代码，因此这是最安全的行为。然而，不需要的代码必须手动删除。

%load 魔法还允许通过提供 URL 作为输入直接从网上加载代码：

%load http://matplotlib.org/mpl_examples/pylab_examples/boxplot_demo2.py

这将从 matplotlib 网站加载一个箱线图示例的代码到一个单元格中。要显示图像，必须在该单元格中运行脚本。

运行其他语言脚本

我们还可以直接在笔记本中运行用其他语言编写的脚本。以下表格列出了一些支持的语言：

单元魔法	语言
%%HTML 或 %%html	HTML
%%SVG 或 %%svg	缩放矢量图形语言（SVGL）
%%bash	Bash 脚本语言，适用于类 Unix 系统，如 Ubuntu 和 Mac OS X
%%cmd	MS Windows 命令行语言
%%javascript	JavaScript
%%latex	LaTeX，面向科学的文档准备语言
%%perl	PERL 脚本语言
%%powershell	MS Windows PowerShell 语言
%%python2 或 %%python3	运行一个与笔记本运行的 Python 版本不同的脚本
%%ruby	Ruby 脚本语言

现在，让我们看一些这些语言的脚本示例。在一个单元格中运行以下代码：

%%SVG
<svg width="400" height="300">
 <circle  cx="200" cy="150" r="100"
 style="fill:Wheat; stroke:SteelBlue; stroke-width:5;"/>
 <line x1="10" y1="10" x2="250" y2="85"
 style="stroke:SlateBlue; stroke-width:4"/>
 <polyline points="20,30 50,70 100,25 200,120"
 style="stroke:orange; stroke-width:3; 
 fill:olive; opacity:0.65;"/>
 <rect x="30" y="150" width="120" height="75"
 style="stroke:Navy; stroke-width:4; fill:LightSkyBlue;"/>
 <ellipse cx="310" cy="220" rx="55" ry="75"
 style="stroke:DarkSlateBlue; stroke-width:4;
 fill:DarkOrange; fill-opacity:0.45;"/>
 <polygon points="50,50 15,100 75,200 45,100"
 style="stroke:DarkTurquoise; stroke-width:5; fill:Beige;"/>
</svg>

这将显示一个由基本形状组成的图形，使用 SVG 描述。SVG 是 HTML 标准，因此此代码将在支持该标准的现代浏览器中运行。

为了演示 JavaScript 的使用，让我们首先在一个计算单元中定义一个可以轻松访问的 HTML 元素：

%%html
<h1 id="hellodisplay">Hello, world!</h1>

运行此单元格。大小为 h1 的 "Hello, world!" 消息将被显示。然后在另一个单元格中输入以下命令：

%%javascript
element = document.getElementById("hellodisplay")
element.style.color = 'blue'

当第二个单元格运行时，"Hello, world!" 消息的文本颜色会从黑色变为蓝色。

该笔记本实际上可以运行系统中已安装的任何脚本语言。这是通过使用 %%script 单元魔法实现的。举个例子，让我们运行一些 Julia 脚本语言的代码。Julia 是一种新的技术计算语言，可以从 julialang.org/ 下载。以下示例假设已安装 Julia，并且可以通过 julia 命令访问（这需要语言解释器的可执行文件在操作系统的路径中）。在一个单元格中输入以下代码并运行：

%%script julia
function factorial(n::Int)
 fact = 1
 for k=1:n
 fact *= k
 end
 fact
end

println(factorial(10))

上面的代码定义了一个函数（用julia编写），该函数计算一个整数的阶乘，并打印出 10 的阶乘。以下是输出结果：

factorial (generic function with 1 method)
3628800

第一行是来自julia解释器的消息，第二行是 10 的阶乘。

数据的加载和保存

数据的加载或保存方式取决于数据的性质以及使用数据的应用程序所期望的格式。由于无法涵盖所有数据结构和应用程序的组合，本节将仅介绍使用NumPy加载和保存数据的最基本方法。Python 中推荐的加载和保存结构化数据的方式是使用针对每种特定数据类型优化的专业库。例如，在处理表格数据时，我们可以使用pandas，正如在第四章中所描述的，使用 pandas 处理数据。

可以通过调用NumPy的save()函数来保存单个数组。以下是一个示例：

A = rand(5, 10)
print A
save('random_array.npy', A)

这段代码生成一个包含五行十列的随机值数组，打印出来，并将其保存到名为random_array.npy的文件中。.npy格式是专门用于NumPy数组的。现在，让我们使用以下命令删除包含该数组的变量：

del A
A

运行包含上述命令的单元格会产生一个错误，因为我们请求显示变量A时，它已经被删除。要恢复该数组，请在一个单元格中运行以下命令：

A = load('random_array.npy')
A

也可以将多个数组保存到一个压缩文件中，下面是一个示例：

xvalues = arange(0.0, 10.0, 0.5)
xsquares = xvalues ** 2
print xvalues
print xsquares
savez('values_and_squares.npz', values=xvalues, squares=xsquares)

注意如何使用关键字参数来指定保存到磁盘的数组的名称。这些数组现在已经以.npz格式保存到文件中。可以使用load()函数从磁盘恢复数据，该函数能够读取NumPy使用的两种格式的文件：

my_data = load('values_and_squares.npz')

如果传递给load()的文件是.npz类型，那么返回值是一个NpzFile类型的对象。这个对象不会立即读取数据。数据读取会延迟，直到需要使用数据时为止。要查看文件中存储了哪些数组，请在一个单元格中执行以下命令：

my_data.files

在我们的例子中，上述命令会产生以下输出：

['squares', 'values']

要将数组分配给变量，请使用 Python 字典访问符号，方法如下：

xvalues = my_data['values']
xsquares = my_data['squares']
plot(xvalues, xsquares)

上述代码生成了半个抛物线的图像：

丰富的显示系统

在一个令人兴奋的发展中，IPython 的最新版本包括了在笔记本中直接显示图像、视频、声音和其他媒体的功能。支持显示系统的类位于IPython.display模块中。在本节中，我们将讨论一些支持的格式。

图片和 YouTube 视频

图像可以从本地文件系统或网络加载。例如，要显示 character.png 文件中的图像，可以在单元格中运行以下命令：

from IPython.display import Image
Image('character.png')

也可以将图像存储在变量中，以便稍后显示：

img = Image('character.png')

character.png 文件可以从本书的网页下载。

要显示图像，我们可以使用 img 或 display(img)。下面的图像会被显示：

要从网络加载图像，只需将其 URL 作为参数提供：

Image('http://www.imagesource.com/Doc/IS0/Media/TR5/7/7/f/4/IS09A9H4K.jpg')

前面的命令行加载了一张花朵的图像：

默认情况下，图像被嵌入在笔记本中，以便可以离线查看。要插入指向图像的链接，请使用以下命令：

Image('http://www.imagesource.com/Doc/IS0/Media/TR5/7/7/f/4/IS09A9H4K.jpg', embed=False)

图像将像前面的示例那样显示，但这次只会在笔记本中插入指向图像的链接。这会使笔记本文件的大小变小，但有一个警告！如果图像在网上发生变化，变化将会反映到笔记本中。

也可以很容易地直接从 YouTube 嵌入视频。以下代码展示了曼德尔布罗集合的美丽动画：

from IPython.display import YouTubeVideo
YouTubeVideo('G_GBwuYuOOs')

HTML

为了完成这一部分，我们展示一个使用 IPython 显示 HTML 功能的扩展示例。这个例子的目标是构建并显示一个数学曲线的 HTML 表格。我们从生成图形并将其保存到磁盘开始：

%matplotlib
xvalues = linspace(-pi,pi,200)
fcts = [('sin', sin), ('cos', cos), ('exp', exp)]
for fctname, fct in fcts:
 yvalues = fct(xvalues)
 fig=figure()
 ax = fig.add_subplot(1,1,1)
 ax.plot(xvalues, yvalues, color='red')
 ax.set_xlabel('$x$')
 strname = '$\\%s(x)$' % fctname
 ax.set_ylabel(strname)
 fig.savefig(fctname + '.png')

该单元格从没有参数的 %matplotlib 魔法命令开始，因为我们不希望将图形嵌入在线（它们仍然会在外部查看器中生成）。接着，我们定义了包含我们要绘制的曲线的 fcts 列表。每条曲线由一个包含两个元素的元组指定：一个是函数的名称，另一个是表示该函数的实际对象。然后，图表将在循环中生成。这里的 Python 代码比我们迄今为止看到的要复杂一些，matplotlib 中使用的函数将在下一章中详细解释。现在，只需注意单元格结尾的调用：

fig.savefig(fctname + '.png')

前面的命令将图形文件保存到磁盘，使用 .png 格式。

接下来，我们生成 HTML 来创建表格，并将其存储在 html_string 变量中，如下所示：

html_string = '<table style="padding:20px">\n'
for fctname, fct in fcts:
 strname = strname = '$\\%s(x)$' % fctname
 filename = fctname + '.png'
 html_string += '<tr>\n'
 html_string += '<td style="width:80px;">%s</td>\n' % strname
 html_string += '<td style="width:500px;">'
 html_string += '<img src="img/%s">' % filename
 html_string += '</td>\n'
 html_string += '</tr>\n'
 html_string += '</table>\n'

HTML 代码是逐步生成的。我们从在单元格的第一行添加 <table> 元素开始。然后，在循环中，每次迭代生成一行表格。为了使代码更易读，我们在每一行代码中只添加一个 HTML 元素。

然后，我们可以打印我们生成的 HTML 以检查它是否正确：

print html_string

前面的命令产生了以下输出：

<table style="padding:20px">
<tr>
<td style="width:80px;">$\sin(x)$</td>
<td style="width:500px;"><img src="img/sin.png"></td>
</tr>
<tr>
<td style="width:80px;">$\cos(x)$</td>
<td style="width:500px;"><img src="img/cos.png"></td>
</tr>
<tr>
<td style="width:80px;">$\exp(x)$</td>
<td style="width:500px;"><img src="img/exp.png"></td>
</tr>
</table>

这似乎是正确的，因此我们准备渲染 HTML：

from IPython.display import HTML
HTML(html_string)

如果一切正确，以下曲线表格将会显示：

在示例结束时，不要忘记运行以下命令以恢复内联图形：

%matplotlib inline

IPython HTML 显示对象功能非常强大，如前面的示例所示。HTML5 丰富媒体，如声音和视频，也可以嵌入其中；不过，当前对所有格式的支持因浏览器而异。

总结

在这一章中，我们全面介绍了 IPython Notebook 界面。我们涵盖了日常使用的功能，如导航、魔法命令、与操作系统的交互、运行脚本以及加载和保存数据。最后，我们讨论了如何在笔记本中显示丰富格式的数据。

在下一章，你将学习如何使用matplotlib库生成高质量的科学图表和数据展示，重点介绍交互式图表。

第三章：matplotlib 图形

本章介绍了 matplotlib，这是一个用于生成出版质量图形的 IPython 库。在本章中，将讨论以下主题：

• 使用 plot 函数绘制二维图，并设置线宽、颜色和样式

• 绘图配置和注释

• 三维图形

• 动画

作为一个 IPython 库，matplotlib 由一系列类构成，且可以按照常规面向对象的方式进行编码。然而，matplotlib 也支持交互式模式。在这种模式下，图形是一步步构建的，每次添加和配置一个组件。我们强调第二种方法，因为它旨在快速生成图形。当面向对象的方式能提供更好的效果时，我们会在必要时进行解释。

注意

在此上下文中，交互式 一词的含义与今天的理解略有不同。由 matplotlib 绘制的图形并不具备交互式特性，即用户不能在图形渲染到笔记本后进行操作。相反，这一术语源于 matplotlib 最初主要用于命令行模式时，每一行新代码都会修改已有的图形。有趣的是，最初启发 matplotlib 的软件仍然使用基于命令行的界面。

绘图函数

plot() 函数是 matplotlib 库的主力功能。在这一部分，我们将探讨该函数包含的线条绘制和格式化功能。

为了让事情更具具体性，我们来看看逻辑增长的公式，如下所示：

该模型常用于表示经历初期指数增长阶段后，最终受某种因素限制的增长。例如，有限资源环境中的人口增长、新产品和/或技术创新的市场，最初吸引一个小而迅速增长的市场，但最终会达到饱和点。

理解数学模型的一个常见策略是研究当定义模型的参数发生变化时，模型如何变化。假设我们想看看当参数 b 变化时，曲线的形状会发生什么变化。

为了更高效地实现我们的目标，我们将使用函数工厂。这样，我们可以快速创建具有任意 r、a、b 和 c 值的逻辑模型。在一个单元格中运行以下代码：

def make_logistic(r, a, b, c):
 def f_logistic(t):
 return a / (b + c * exp(-r * t))
 return f_logistic

函数工厂模式利用了函数在 Python 中是一等对象的事实。这意味着函数可以像普通对象一样被处理：它们可以赋值给变量、存储在列表或字典中，并且可以作为参数或/和返回值传递给其他函数。

在我们的例子中，我们定义了make_logistic()函数，其输出本身是一个 Python 函数。注意，f_logistic()函数是在make_logistic()的函数体内定义的，然后在最后一行返回。

现在我们使用函数工厂创建三个表示逻辑曲线的函数，如下所示：

r = 0.15
a = 20.0
c = 15.0
b1, b2, b3 = 2.0, 3.0, 4.0
logistic1 = make_logistic(r, a, b1, c)
logistic2 = make_logistic(r, a, b2, c)
logistic3 = make_logistic(r, a, b3, c)

在上述代码中，我们首先固定了r、a和c的值，并为不同的b值定义了三条逻辑曲线。需要注意的重要点是，logistic1、logistic2和logistic3是函数。因此，例如，我们可以使用logistic1(2.5)来计算时间为 2.5 时的第一个逻辑曲线值。

现在我们可以使用以下代码绘制函数图像：

tmax = 40
tvalues = linspace(0, tmax, 300)
plot(tvalues, logistic1(tvalues)) 
plot(tvalues, logistic2(tvalues)) 
plot(tvalues, logistic3(tvalues))

上述代码的第一行将最大时间值tmax设置为40。接着，我们通过赋值定义了希望评估函数的时间集合，如下所示：

tvalues = linspace(0, tmax, 300)

linspace()函数非常方便用来生成绘图的点。上述代码创建了一个包含 300 个均匀分布点的数组，区间从0到tmax。注意，与range()和arange()等其他函数不同，区间的右端点默认是包括的。（要排除右端点，请使用endpoint=False选项。）

在定义了时间值的数组后，调用plot()函数来绘制曲线。在最基本的形式中，它绘制一条默认颜色和线型的单条曲线。在此用法中，两个参数是两个数组。第一个数组提供被绘制点的横坐标，第二个数组提供纵坐标。一个典型的例子是以下函数调用：

plot(x,y)

变量x和y必须是NumPy数组（或任何可以转换为数组的 Python 可迭代值），并且必须具有相同的维度。绘制的点具有以下坐标：

x[0], y[0]
x[1], y[1]
x[2], y[2]
…

上述命令将产生以下图形，显示三条逻辑曲线：

你可能已经注意到，在显示图形之前，会有一行看起来像这样的文本输出：

[<matplotlib.lines.Line2D at 0x7b57c50>]

这是最后一次调用plot()函数的返回值，它是一个包含Line2D类型对象的列表（或者只有一个元素）。阻止输出显示的一种方式是，在单元格的最后一行输入None。另外，我们还可以将单元格中最后一次调用的返回值赋值给一个虚拟变量：

_dummy_ = plot(tvalues, logistic3(tvalues))

plot()函数支持在同一个函数调用中绘制多条曲线。我们需要更改下面代码单元格中的内容，并再次运行它：

tmax = 40
tvalues = linspace(0, tmax, 300)
plot(tvalues, logistic1(tvalues), 
 tvalues, logistic2(tvalues), 
 tvalues, logistic3(tvalues))

这种形式可以节省一些输入，但在定制线条选项时，灵活性稍差。请注意，现在产生的文本输出是一个包含三个元素的列表：

[<matplotlib.lines.Line2D at 0x9bb6cc0>,
 <matplotlib.lines.Line2D at 0x9bb6ef0>,
 <matplotlib.lines.Line2D at 0x9bb9518>]

这种输出在某些情况下可能有用。目前，我们将继续为每条曲线使用一次plot()调用，因为它生成的代码更清晰且更灵活。

现在，让我们更改图中的线条选项并设置图形的边界。将单元格内容更改为如下所示：

plot(tvalues, logistic1(tvalues), 
 linewidth=1.5, color='DarkGreen', linestyle='-') 
plot(tvalues, logistic2(tvalues), 
 linewidth=2.0, color='#8B0000', linestyle=':') 
plot(tvalues, logistic3(tvalues), 
 linewidth=3.5, color=(0.0, 0.0, 0.5), linestyle='--')
axis([0, tmax, 0, 11.])
None

运行上述命令行将生成以下图形：

在前面的代码中设置的选项如下：

• 第一条曲线的线宽为1.5，颜色为 HTML 颜色DarkGreen，并且使用实线样式。

• 第二条曲线的线宽为2.0，颜色由十六进制字符串#8B0000指定，并采用点线样式。

• 第三条曲线的线宽为3.0，颜色由 RGB 组件(0.0, 0.0, 0.5)指定，且采用虚线样式。

注意，有不同的方式来指定固定颜色：HTML 颜色名称、十六进制字符串或浮点值元组。在最后一种情况下，元组中的条目分别代表红色、绿色和蓝色的强度，并且必须是介于0.0和1.0之间的浮点值。完整的 HTML 颜色名称列表可以在www.w3schools.com/html/html_colornames.asp找到。

线条样式由符号字符串指定。允许的值如下表所示：

符号字符串	线条样式
'-'	实线（默认）
'--'	虚线
':'	点线
'-.'	虚线点线
'None'、'' 或 ''	不显示

在调用plot()之后，我们使用以下函数来设置图形的边界：

axis([0, tmax, 0, 11.])

axis()的参数是一个包含四个元素的列表，按照顺序指定水平坐标的最大值和最小值，以及垂直坐标的最大值和最小值。

看起来在绘制图形之后再设置变量的边界似乎不直观。在交互模式下，matplotlib 会记住正在构建的图形的状态，并且在每次发出命令后，图形对象会在后台更新。图形仅在单元格中的所有计算完成后才会渲染，以便所有先前指定的选项生效。请注意，开始一个新单元格会清除所有图形数据。这种交互行为是matplotlib.pyplot模块的一部分，而pylab会导入该模块。

除了绘制连接数据点的线条外，还可以在指定的点上绘制标记。更改以下代码片段中所示的绘图命令，然后重新运行该单元格：

plot(tvalues, logistic1(tvalues), 
 linewidth=1.5, color='DarkGreen', linestyle='-',
 marker='o', markevery=50, markerfacecolor='GreenYellow',
 markersize=10.0) 
plot(tvalues, logistic2(tvalues), 
 linewidth=2.0, color='#8B0000', linestyle=':',
 marker='s', markevery=50, markerfacecolor='Salmon',
 markersize=10.0) 
plot(tvalues, logistic3(tvalues), 
 linewidth=2.0, color=(0.0, 0.0, 0.5), linestyle='--',
 marker = '*', markevery=50, markerfacecolor='SkyBlue',
 markersize=12.0)
axis([0, tmax, 0, 11.])
None

现在，图形将如下面的图所示：

与之前的代码唯一的区别是，现在我们添加了绘制标记的选项。以下是我们使用的选项：

• marker选项指定了标记的形状。形状由符号字符串表示。在之前的示例中，我们使用'o'表示圆形标记，'s'表示方形标记，'*'表示星形标记。可用标记的完整列表可以在matplotlib.org/api/markers_api.html#module-matplotlib.markers查看。

• markevery选项指定了数据点之间标记的间隔。在我们的示例中，我们每 50 个数据点放置一个标记。

• markercolor选项指定了标记的颜色。

• markersize选项指定了标记的大小，大小以像素为单位。

matplotlib 中有大量其他选项可以应用于线条。完整的列表可以在matplotlib.org/api/artist_api.html#module-matplotlib.lines查看。

添加标题、标签和图例

下一步是添加标题和坐标轴标签。在None行之前，向创建图表的代码单元格中添加以下三行代码：

title('Logistic growth: a={:5.2f}, c={:5.2f}, r={:5.2f}'.format(a, c, r))
xlabel('$t$')
ylabel('$N(t)=a/(b+ce^{-rt})$')

在第一行，我们调用title()函数来设置图表的标题。参数可以是任何 Python 字符串。在我们的示例中，我们使用了格式化字符串：

title('Logistic growth: a={:5.2f}, b={:5.2f}, r={:5.2f}'.format(a, c, r)) 

我们使用字符串类的format()方法。格式放置在大括号之间，例如{:5.2f}，它指定一个浮点格式，具有五个空间和两位精度。每个格式说明符依次与方法的一个数据参数关联。字符串格式化的部分细节可以在附录 B 中查看，完整的文档可以在docs.python.org/2/library/string.html查看。

坐标轴标签通过以下调用设置：

xlabel('$t$')
ylabel('$N(t)=a/(b+ce^{-rt})$')

与title()函数一样，xlabel()和ylabel()函数接受任何 Python 字符串。注意，在'$t$'和'$N(t)=a/(b+ce^{-rt}$'字符串中，我们使用 LaTeX 格式化数学公式。公式中的美元符号$...$表示 LaTeX 格式。

添加标题和标签后，我们的图表如下所示：

接下来，我们需要一种方法来识别图中的每条曲线。一种方法是使用legend，如下面所示：

legend(['b={:5.2f}'.format(b1),
 'b={:5.2f}'.format(b2),
 'b={:5.2f}'.format(b3)])

legend()函数接受一个字符串列表。每个字符串与曲线相关联，顺序是根据它们添加到图表的顺序。注意，我们再次使用了格式化字符串。

不幸的是，之前的代码并没有产生很好的效果。图例默认放置在图表的右上角，这在此情况下会遮挡图表的一部分。我们可以通过在legend函数中使用loc选项来轻松修复，如下所示：

legend(['b={:5.2f}'.format(b1),
 'b={:5.2f}'.format(b2),
 'b={:5.2f}'.format(b3)], loc='upper left')

运行此代码后，我们得到了最终版本的逻辑增长图，如下所示：

图例位置可以是以下字符串中的任何一个：'best'、'upper right'、'upper left'、'lower left'、'lower right'、'right'、'center left'、'center right'、'lower center'、'upper center' 和 'center'。还可以通过bbox_to_anchor选项精确指定图例的位置。为了查看其工作原理，请按如下方式修改图例的代码：

legend(['b={:5.2f}'.format(b1),
 'b={:5.2f}'.format(b2),
 'b={:5.2f}'.format(b3)],  bbox_to_anchor=(0.9,0.35))

请注意，bbox_to_anchor选项默认使用的坐标系统与我们为绘图指定的坐标系统不同。在前面的例子中，框的x和y坐标分别被解释为整个图形宽度和高度的比例。为了精确地放置图例框，需要做一些反复试验。请注意，图例框可以放置在绘图区域之外。例如，可以尝试坐标(1.32,1.02)。

legend()函数非常灵活，具有许多其他选项，相关文档可以参考matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.legend。

文本和注释

在这一小节中，我们将展示如何在 matplotlib 中向图形添加注释。我们将构建一个图形，演示曲线的切线在最高点和最低点必须是水平的这一事实。我们首先定义与曲线相关的函数，以及我们希望绘制该曲线的一组值，具体代码如下：

f = lambda x:  (x**3 - 6*x**2 + 9*x + 3) / (1 + 0.25*x**2)
xvalues = linspace(0, 5, 200)

前面代码的第一行使用 lambda 表达式定义了f()函数。我们在这里采用这种方法，因为函数的公式是简单的一行表达式。lambda 表达式的一般形式如下：

lambda <arguments> : <return expression>

这个表达式本身创建了一个匿名函数，可以在任何需要函数对象的地方使用。请注意，返回值必须是单一的表达式，不能包含任何语句。

这个函数的公式看起来可能有些不寻常，但它是通过试错法和一些微积分选择的，以便在区间0到5内生成一个漂亮的图形。xvalues数组被定义为包含该区间内 200 个均匀分布的点。

让我们先创建一个初始的曲线图，如下所示的代码：

plot(xvalues, f(xvalues), lw=2, color='FireBrick')
axis([0, 5, -1, 8])
grid()
xlabel('$x$')
ylabel('$f(x)$')
title('Extreme values of a function')
None # Prevent text output

这一段代码中的大部分内容在前一节已经解释过。唯一新增的部分是我们使用grid()函数来绘制网格。没有参数时，网格与图中的刻度线重合。与 matplotlib 中的其他内容一样，网格是高度可定制的。请参阅文档matplotlib.org/1.3.1/api/pyplot_api.html#matplotlib.pyplot.grid。

当执行前面的代码时，会生成如下图形：

注意，曲线有一个最高点（最大值）和一个最低点（最小值）。这些点合称为函数的极值（在所显示的区间内，实际上当x变大时，函数会无限增大）。我们希望在图表上标出这些点并添加注释。首先，我们将相关点存储如下：

x_min = 3.213
f_min = f(x_min)
x_max = 0.698
f_max = f(x_max)
p_min = array([x_min, f_min])
p_max = array([x_max, f_max])
print p_min
print p_max

变量 x_min 和 f_min 定义为图表中最低点的（大约）坐标。类似地，x_max 和 f_max 代表最高点。无需关心这些点是如何找到的。对于绘图而言，即使是通过试错法得到的粗略近似值也足够了。在第五章中，我们将看到如何通过 SciPy 精确地解决这类问题。现在，在绘制图表的单元格中，紧接着 title() 命令，添加以下代码，如下所示：

arrow_props = dict(facecolor='DimGray', width=3, shrink=0.05, 
 headwidth=7)
delta = array([0.1, 0.1])
offset = array([1.0, .85])
annotate('Maximum', xy=p_max+delta, xytext=p_max+offset,
 arrowprops=arrow_props, verticalalignment='bottom',
 horizontalalignment='left', fontsize=13)
annotate('Minimum', xy=p_min-delta, xytext=p_min-offset,
 arrowprops=arrow_props, verticalalignment='top',
 horizontalalignment='right', fontsize=13)

运行单元格以生成下图所示的图表：

在代码中，首先赋值变量 arrow_props、delta 和 offset，它们将用于设置 annotate() 函数的参数。annotate() 函数将文本注释添加到图表上，箭头可选地指示被注释的点。该函数的第一个参数是注释文本。接下来的两个参数给出了箭头和文本的位置：

• xy：这是被注释的点，并将对应于箭头的尖端。我们希望这个点是最大值/最小值点 p_min 和 p_max，但我们会加上/减去 delta 向量，以便箭头的尖端略微偏离实际点。

• xytext：这是文本的位置和箭头的基点。我们通过 offset 向量指定该点相对于 p_min 和 p_max 的偏移量。

annotate() 的所有其他参数都是格式化选项：

• arrowprops：这是一个包含箭头属性的 Python 字典。我们预定义了字典 arrow_props 并在此使用它。箭头在 matplotlib 中可以非常复杂，具体细节请参考文档。

• verticalalignment 和 horizontalalignment：这些指定了箭头与文本的对齐方式。

• fontsize：这表示文本的大小。文本也可以高度自定义，读者可参考文档获取详细信息。

annotate() 函数有大量的选项；欲了解所有可用的详细信息，用户应参考文档 matplotlib.org/1.3.1/api/pyplot_api.html#matplotlib.pyplot.annotate 获取完整资料。

现在，我们想要通过添加一个解释性文本框来对图形展示的内容进行注释。在调用annotate()之后，向单元格中添加以下代码：

bbox_props = dict(boxstyle='round', lw=2, fc='Beige')
text(2, 6, 'Maximum and minimum points\nhave horizontal tangents', 
 bbox=bbox_props, fontsize=12, verticalalignment='top')

text()函数用于将文本放置在图形的任意位置。前两个参数是文本框的位置，第三个参数是一个包含要显示文本的字符串。注意使用'\n'表示换行。其他参数是配置选项。bbox参数是一个字典，用于设置文本框的选项。如果省略，文本将不显示任何边框。在示例代码中，文本框是一个具有圆角的矩形，边框宽度为 2 像素，面色为米色。

作为最后的细节，让我们在极端点处添加切线。添加以下代码：

plot([x_min-0.75, x_min+0.75], [f_min, f_min],
 color='RoyalBlue', lw=3)
plot([x_max-0.75, x_max+0.75], [f_max, f_max],
 color='RoyalBlue', lw=3)

由于切线是直线段，我们只需要给出端点的坐标。之所以将切线的代码放在单元格顶部，是因为这样可以确保切线先被绘制，从而使得函数图形绘制在切线之上。以下是最终结果：

到目前为止，我们看到的例子仅仅触及了 matplotlib 的冰山一角。读者应查阅 matplotlib 文档以获得更多示例。

三维图形

在本节中，我们介绍了展示三维图形的方法，即空间中数学对象的图形。示例包括不局限于平面上的曲面和线条。

matplotlib 对三维图形提供了出色的支持。在本节中，我们将展示一个表面图和相应的等高线图的示例。三维库中可用的图形类型包括框架图、线图、散点图、三角形网格表面图、多边形图等。以下链接将帮助你了解这里没有介绍的图形类型：matplotlib.org/1.3.1/mpl_toolkits/mplot3d/tutorial.html#mplot3d-tutorial

在开始之前，我们需要使用以下命令行导入所需的三维库对象：

from mpl_toolkits.mplot3d import axes3d

现在，让我们通过在单元格中运行以下代码来绘制表面图：

def dist(x, y):
 return sqrt(x**2 + y**2) 
def fsurface(x, y):
 d = sqrt(x**2 + y**2)
 c = 5.0
 r = 7.5
 return c - (d**4 - r * d**2)

xybound = 2.5
fig = figure(figsize=(8,8))
ax = subplot(1, 1, 1, projection='3d')
X = linspace(-xybound, xybound, 25)
Y = linspace(-xybound, xybound, 25)
X, Y = meshgrid(X, Y)
Z = fsurface(X,Y)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
 antialiased=True, linewidth=0.2)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$f(x,y)$')
None #Prevent text output

我们从指定fsurface()函数开始，该函数定义了表面。函数定义的细节并不重要；我们只需要注意这是一个旋转曲面，中心有一个低点，四周被一个高峰包围。然后，我们开始通过以下几行代码来设置图形：

fig = figure(figsize=(8,8))
ax = subplot(1, 1, 1, projection='3d')

这一次，我们特别构造Figure对象，因为我们想要明确指定其大小。这里的大小定义为一个8 x 8英寸的正方形，但实际图形的大小取决于显示的分辨率和浏览器的放大倍数。然后我们创建一个子图并将其投影类型设置为'3d'。subplot()函数将在本节中详细介绍。

接下来，我们将定义计算函数的点网格：

xybound = 2.5
x = linspace(-xybound, xybound, 25)
y = linspace(-xybound, xybound, 25)
X, Y = meshgrid(x, y)

这里最重要的一点是使用meshgrid()函数，它是NumPy包的一部分。此函数接受两个一维数组，带有x和y值，并计算定义平面上相应点的网格的两个矩阵。要理解其工作原理，请运行以下代码：

xx = [1,2,3]
yy = [4,5,6]
XX, YY = meshgrid(xx, yy)
print XX
print YY

生成的两个矩阵，XX和YY，如下所示：

• XX矩阵：

  [[1 2 3]
   [1 2 3]
   [1 2 3]]
  

• YY矩阵：

  [[4 4 4]
   [5 5 5]
   [6 6 6]]
  

请注意，如果我们取XX的元素和YY中对应的条目，我们得到一组点(1,4)、(1,5)、(1,6)、(2,4)...、(3,5)、(3,6)，这些点在平面上是规则间隔的网格。

现在我们可以调用计算表面并绘制它的函数：

Z = fsurface(X,Y)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
 antialiased=True, linewidth=0.2)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$f(x,y)$')

第一行计算Z数组，其中包含表面的z坐标。这个调用在背景中大量使用了NumPy称为广播的特性。这是一组规则，告诉我们NumPy如何处理具有不同大小的数组的操作。更多信息请参见docs.scipy.org/doc/numpy/user/basics.broadcasting.html。

下一步是调用plot_surface()方法，它实际上进行绘图。前三个参数定义了正在绘制的数据，即数组X、Y和Z。cstride和rstride选项可用于跳过数据数组中的点。将这些值设置为大于 1 的值以跳过网格中的点，以防数据集过大。

我们正在使用一个由cmap=cm.coolwarm选项指定的色彩映射特性。色彩映射功能告诉 matplotlib 如何为绘图中的每个高度分配颜色。提供了大量内置的色彩映射。要查看完整列表，请在单元格中运行以下代码：

for key, value in cm.__dict__.items():
 if isinstance(value, matplotlib.colors.Colormap):
 print key

请注意，默认情况下，三维表面图不是抗锯齿的，因此我们在代码中设置了antialiased=True选项以生成更好的图像。

现在让我们在图中添加等高线图。我们希望三维表面图和等高线图并排显示。为了实现这一点，请将单元格中的代码修改如下：

fig = figure(figsize(20,8))
ax1 = subplot(1, 2, 1, projection='3d')
X = linspace(-xybound, xybound, 100)
Y = linspace(-xybound, xybound, 100)
X, Y = np.meshgrid(X, Y)
Z = fsurface(X,Y)
ax1.plot_surface(X, Y, Z, rstride=5, cstride=5, cmap=cm.coolwarm,
 antialiased=True, linewidth=0.2)
ax1.set_xlabel(r'$x$')
ax1.set_ylabel(r'$y$')
ax1.set_zlabel(r'$f(x,y)$')
ax1.set_title('A surface plot', fontsize=18)
ax2 = subplot(1, 2, 2)
ax2.set_aspect('equal')
levels = arange(5, 20, 2.5)
cs = ax2.contour(X, Y, Z,
 levels,
 cmap=cm.Reds,
 linewidths=1.5)
cs.clabel(levels[1::2], fontsize=12)
ax2.set_title('Contour Plot', fontsize=18)

运行代码的结果如下图所示：

首先，我们集中讨论contours()方法。第一个参数levels指定了绘制等高线的值（高度）。这个参数可以省略，matplotlib 会尝试选择合适的高度。其他参数是控制如何显示等高线的选项。我们在这个示例中指定了色图和线宽。

clabel()方法为等高线添加标签。第一个参数levels[1::2]指定每隔一个等高线就添加标签。

注意用于将两个坐标轴放置在同一图形中的代码。坐标轴通过以下命令行定义：

ax1 = subplot(1, 2, 1, projection='3d')
ax2 = subplot(1, 2, 2)

subplot()函数的一般形式如下：

subplot(nrows, ncols, axis_position, **kwargs)

这指定了一个具有nrows行和ncols列的数组中的Axes对象。坐标轴的位置是一个从 1 到nrows*ncols的整数。下图演示了在 3 x 2 数组情况下坐标轴的编号方式：

上图是通过以下命令行生成的：

fig = figure(figsize=(5,6))
nrows = 3
ncols = 2
for i in range(nrows*ncols):
 ax = subplot(nrows, ncols, i+1, axisbg='Bisque')
 axis([0,10,0,5])
 text(1, 2.2, 'subplot({},{},{})'.format(nrows, ncols, i+1),
 fontsize=14, color='Brown')

在按照常规方式定义图形大小后，我们设置了我们希望生成的坐标轴网格的行数和列数。然后，在循环中创建并配置每个Axes对象。仔细查看如何确定坐标轴的位置。还请注意，我们展示了如何设置坐标轴的背景颜色。

动画

我们将在本章结束时给出一个更复杂的示例，展示 matplotlib 的强大功能。我们将创建一个强迫摆锤的动画，这是一个著名且广泛研究的动态系统示例，表现出确定性混沌。

由于本节涉及更复杂的代码，我们将避免使用pylab，并采用通常推荐的模块导入方式。这使得代码在我们需要时更容易导出为脚本。我们还给出了 matplotlib 的一些面向对象功能的示例。

动画摆锤（或任何物理过程）的过程实际上非常简单：我们在有限的时间点计算摆锤的位置，并快速展示相应的图像。因此，代码自然会分解为以下三部分：

• 一个显示摆锤在任意位置的函数

• 设置在任意时刻计算摆锤位置的代码

• 实际计算摆锤位置并显示相应图像的代码

我们首先导入设置动画所需的所有模块和函数：

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.patches as patches
import matplotlib.lines as lines
from scipy.integrate import ode
from IPython.display import display, clear_output
import time

在以下代码中，我们将定义一个绘制简单摆锤草图的函数：

def draw_pendulum(ax, theta, length=5, radius=1):
 v = length * np.array([np.sin(theta), -np.cos(theta)])
 ax.axhline(0.0, color='Brown', lw=5, zorder=0)
 rod = lines.Line2D([0.0, v[0]], [0.0, v[1]], 
 lw=4, color='DarkGray', zorder=1)
 bob = patches.Circle(v, radius, 
 fc='DodgerBlue', ec='DarkSlateGray', 
 lw=1.5, zorder=3)
 peg = patches.Circle([0.0, 0.0], 0.3, 
 fc='black', zorder=2)
 ax.add_patch(bob)
 ax.add_patch(peg)
 ax.add_line(rod)
 return ax

该函数的第一个参数是一个Axes对象。其他参数如下：

• 摆锤与垂直面之间的角度theta

• 杆的length

• 摆锤的radius

上述量在下图中进行了指示：

然后，我们定义一个NumPy向量v，它表示摆锤相对于原点的位置。以下语句定义了要绘制的对象：

• ax.axhline()：此函数在图中绘制一条水平线

• rod：这是一个lines.Line2D对象（顺便提一下，这是用于绘制大多数 matplotlib 图表的对象）

• bob和peg：这些是patches.Circle类型的对象；matplotlib 的 patches 代表任何可以放置在图中的对象

以下代码行可用于测试绘图代码：

fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
ax.set_xlim(-10,10)
ax.set_ylim(-20,0.5)
draw_pendulum(ax, np.pi / 10, length=15, radius=0.5)
ax.set_xticks([])
ax.set_yticks([])

运行上一单元格中的代码将生成以下图像：

以下注释说明了前面示例中代码的工作原理：

1. 前两行定义了保存Figure和Axes对象的变量fig和ax。在 matplotlib 中，Figure对象是一个容器，包含所有其他绘图对象。每个Figure可以包含多个Axes，每个Axes包含单独的图表。注意使用figsize=(5,5)参数来设置图表的大小。

2. 接下来，我们设置坐标轴的限制。ax对象的set_aspect()方法用于将纵横比设置为相同。如果没有这个设置，圆形会被绘制成椭圆。然后，set_xlim()和set_ylim()方法分别指定坐标轴的范围。

3. 然后，我们调用draw_pendulum()函数，它完成所有的绘制工作。

4. 最后，我们使用ax.set_xticks([])和set_yticks([])从坐标轴中去除刻度标记。

下一步是设置所需内容，以找到摆锤的轨迹。摆锤的动力学由一组微分方程给出，定义在以下代码行中：

def pendulum_eq(t, x, c=1, f=lambda t:0.0):
 theta, omega = x
 return np.array([omega,
 -np.sin(theta) - c * omega + f(t)], 
 dtype=np.float64)
solver = ode(pendulum_eq)
solver.set_integrator('lsoda')
c = 0.3
f = lambda t: 2.0 * np.sin(3 * t)
solver.set_f_params(c, f)

这段代码首先定义了pendulum_eq()函数，该函数规定了摆锤的微分方程。方程的推导超出了本书的范围。单元格中的其余代码配置了一个ode()类型的对象，这是scipy.integrate模块的一部分。我们在这里不讨论细节，但该模块在第五章中有介绍，使用 SciPy、Numba 和 NumbaPro 的高级计算。

现在，我们准备通过在单元格中执行以下代码来运行动画：

tmax = 20.0
dt = 0.2

fig = plt.figure(1,figsize=(5,5))
ax = plt.subplot(1,1,1)
ax.set_aspect('equal')
ax.set_xlim(-10,10)
ax.set_ylim(-20,0.5)
xtext = -9
ytext = -2
dytext = -1.0

ic = np.array([0.0, 0.3])
solver.set_initial_value(ic, 0.0)
while True:
 clear_output(wait=True)
 time.sleep(1./20)
 t = solver.t
 if t > tmax:
 break
 theta, omega = solver.integrate(t + dt)
 if not solver.successful():
 print 'Solver interrupted'
 break
 ax.clear()
 ax.set_xticks([])
 ax.set_yticks([])
 ax.text(xtext, ytext, r'$t={:5.2f}$'.format(t))
 ax.text(xtext, ytext + dytext, 
 r'$\theta={:5.2f}$'.format(theta))
 ax.text(xtext, ytext + 2 * dytext, 
 r'$\dot{{\theta}}={:5.2f}$'.format(omega))
 draw_pendulum(ax, theta=theta, length=15, radius=0.5)
 display(fig)
 plt.close()

这可能是到目前为止书中最复杂的代码段，但大部分内容已经涉及。变量tmax和dt分别保存动画的最大时间和时间增量。然后，我们设置绘图的Figure和Axes对象。

接下来是一个while循环，在这个循环中，动画实际上被执行。以下是循环的一般框架：

while True:
 clear_output(wait=True)
 time.sleep(1./20)
 t = solver.t
 if t > tmax:
 break

 ... Code to determine the position of the pendulum...

 ... Code to draw pendulum ...

 display(fig)
 plt.close()

注意

我们不会详细讨论用来解微分方程的代码，因为它将在第五章中详细介绍，使用 SciPy、Numba 和 NumbaPro 进行高级计算。

循环条件是True，因此这可能是一个无限循环。然而，在循环内，我们检查当前时间是否大于动画的最大时间，如果是这种情况，则退出循环：

 t = solver.t
 if t > tmax:
 break

在循环中，我们做的第一件事是调用clear_output()函数。正如函数名所示，它会移除当前单元格的输出，这是在笔记本中进行简单动画的核心。wait=True参数告诉该函数在下一个图像完全绘制完毕后再清除输出，从而防止闪烁。

time.sleep(1./20)参数将计算暂停一段短时间，以防止动画运行过快。然后，计算并绘制摆的当前位置。接着，调用display(fig)显示图形。这里需要这样做，因为与静态图形不同，我们不希望图形仅在单元格的最后显示出来。

最后的细节是在循环结束时调用plt.close()。这可以防止在退出循环时绘制摆的图像多绘制一次。将此调用放在循环内还帮助避免闪烁现象。

鼓励读者调整动画的参数，特别是时间间隔dt、最大时间tmax和time.sleep()参数。需要一些反复试验才能获得令人满意的动画效果。

总结

在本章中，我们学习了如何使用 matplotlib 制作演示质量的图形。我们介绍了二维图形，如何设置图形选项、注释以及配置图形。你还学习了如何添加标签、标题和图例。我们还学习了如何绘制三维曲面图以及如何创建简单的动画。

在下一章中，我们将探索如何使用 pandas 库处理笔记本中的数据。

第四章：使用 pandas 处理数据

在本章中，我们将介绍 pandas，一个功能强大且多用途的 Python 库，提供了数据处理和分析的工具。我们将详细介绍 pandas 中用于存储数据的两个主要结构：Series 和 DataFrame 对象。你将学习如何创建这些结构以及如何访问和插入数据。我们还将介绍一个重要的主题：切片，即如何使用 pandas 提供的不同索引方法访问数据的部分内容。接下来，我们将讨论 pandas 提供的计算和图形工具，并通过展示如何使用实际数据集来结束本章。

pandas 是一个广泛的数据处理包，涵盖的内容超出了本书的范围。我们将仅介绍一些最有用的数据结构和功能。特别地，我们不会覆盖 Panel 数据结构和多重索引。然而，我们将为希望通过查阅官方文档扩展知识的读者打下坚实的基础。在本章中，我们假设进行了以下导入：

%pylab inline
from pandas import Series, DataFrame
import pandas as pd

Series 类

Series 对象表示一个一维的、带索引的数据系列。它可以看作是一个字典，主要的不同在于：Series 类中的索引是有序的。以下示例构造了一个 Series 对象并显示它：

grades1 = Series([76, 82, 78, 100],
 index = ['Alex', 'Robert', 'Minnie', 'Alice'],
 name = 'Assignment 1', dtype=float64)
grades1

这会产生以下输出：

Alex       76
Robert     82
Minnie     78
Alice     100
Name: Assignment 1, dtype: float64

请注意构造函数调用的格式：

Series(<data>, index=<indexes>, name=<name>, dtype=<type>)

data 和 indexes 通常是列表或 NumPy 数组，但也可以是任何 Python 可迭代对象。列表必须具有相同的长度。name 变量是一个描述序列中数据的字符串。type 变量是一个 NumPy 数据类型。indexes 和 name 变量是可选的（如果省略 indexes，则默认为从 0 开始的整数）。数据类型也是可选的，在这种情况下，它将从数据中推断出来。

Series 对象支持标准字典接口。举个例子，在一个单元格中运行以下代码：

print grades1['Minnie']
grades1['Minnie'] = 80
grades1['Theo'] = 92
grades1

前面的命令行的输出如下：

78.0
Alex       76
Robert     82
Minnie     80
Alice     100
Theo       92
Name: Assignment 1, dtype: float64

这是另一个有趣的示例：

for student in grades1.keys():
 print '{} got {} points in {}'.format(student, grades1[student], grades1.name)

前面的命令行产生以下输出：

Alex got 76.0 points in Assignment 1
Robert got 82.0 points in Assignment 1
Minnie got 80.0 points in Assignment 1
Alice got 100.0 points in Assignment 1
Theo got 92.0 points in Assignment 1

请注意，输出的顺序与元素在序列中插入的顺序完全相同。与标准的 Python 字典不同，Series 对象会跟踪元素的顺序。事实上，元素可以通过整数索引进行访问，如以下示例所示：

grades1[2]

前面的命令返回以下输出：

80.0

实际上，Python 的所有列表访问接口都被支持。例如，我们可以使用切片，它返回 Series 对象：

grades1[1:-1]

前面的命令给出以下输出：

Robert     82
Minnie     80
Alice     100
Name: Assignment 1, dtype: float64

索引功能更加灵活；以下示例展示了这一点：

grades1[['Theo', 'Alice']]

前面的命令返回以下输出：

Theo      92
Alice    100
dtype: float64

也可以通过以下命令将新数据附加到序列中：

grades1a = grades1.append(Series([79, 81], index=['Theo', 'Joe']))
grades1a

前面命令的输出如下：

Alex       76
Robert     82
Minnie     80
Alice     100
Theo       92
Kate       69
Molly      74
Theo       79
Joe        81
dtype: float64

请注意，序列现在包含与键Theo相关的两个条目。这是有意义的，因为在现实数据中，可能会有多个数据值与相同的索引相关联。在我们的例子中，一个学生可能提交了多版本的作业。当我们尝试访问这些数据时会发生什么呢？pandas 方便地返回一个Series对象，因此没有丢失任何数据：

grades1a['Theo']

前面命令的输出如下：

Theo    92
Theo    79
dtype: float64

注意

请注意，append()方法并不会将值附加到现有的Series对象中。相反，它会创建一个新的对象，该对象由原始Series对象和附加的元素组成。这种行为与向 Python 列表中附加元素的行为不同。Series类的许多方法表现出的行为与它们相应的列表方法不同。需要进行一些实验（或阅读文档）才能理解 pandas 所采用的约定。

我们来定义一个新的序列，使用以下命令行：

grades2 = Series([87, 76, 76, 94, 88],
 index = ['Alex', 'Lucy', 'Robert', 'Minnie', 'Alice'],
 name='Assignment 2',
 dtype=float64)
grades2

前面的命令行给出以下输出：

Alex      87
Lucy      76
Robert    76
Minnie    94
Alice     88
Name: Assignment 2, dtype: float64

如果我们想要计算每个学生在两个作业中的平均分，可以使用以下命令：

average = 0.5 * (grades1 + grades2)
average

运行前面的代码后，我们得到以下输出：

Alex      81.5
Alice     94.0
Lucy       NaN
Minnie    87.0
Robert    79.0
Theo       NaN
dtype: float64

值NaN代表非数值，它是一个特殊的浮点值，用于表示无效操作的结果，例如零除以零。在 pandas 中，它用于表示缺失的数据值。我们可以使用isnull()方法定位Series中的缺失值。例如，在单元格中运行以下代码：

averages.isnull()

运行前面的命令行会产生以下输出：

Alex      False
Alice     False
Lucy       True
Minnie    False
Robert    False
Theo       True
dtype: bool

如果我们决定可以安全地从序列中删除缺失数据，可以使用dropna()方法：

average.dropna()

前面的命令行产生以下输出：

Alex      81.5
Alice     94.0
Minnie    87.0
Robert    79.0
dtype: float64

请注意，这是另一个原始序列未被修改的情况。

Series类为其实例提供了一系列有用的方法。例如，我们可以对值和索引进行排序。要就地排序值，我们使用sort()方法：

grades1.sort()
grades1

这会生成以下输出：

Alex       76
Minnie     80
Robert     82
Theo       92
Alice     100
Name: Assignment 1, dtype: float64

要排序序列的索引，使用sort_index()方法。例如，考虑以下命令：

grades1.sort_index()

这会产生以下输出：

Alex       76
Minnie     80
Robert     82
Theo       92
Alice     100
Name: Assignment 1, dtype: float64

注意

请注意，这次排序不是就地进行的，返回了一个新的序列对象。

对于接下来的示例，我们将使用来自作者所在地附近气象站的 6 月每日最高气温数据。以下命令行生成了 6 月 6 日至 6 月 15 日的温度序列：

temps = Series([71,76,69,67,74,80,82,70,66,80],
 index=range(6,16), 
 name='Temperatures', dtype=float64)
temps

前面的命令产生以下输出：

6     71
7     76
8     69
9     67
10    74
11    80
12    82
13    70
14    66
15    80
Name: Temperatures, dtype: float64

让我们首先使用以下命令计算温度的均值和标准差：

print temps.mean(), temps.std()

前面计算的结果如下：

73.5 5.77831194112

如果我们想快速了解系列数据的概况，可以使用describe()方法：

temps.describe()

上述命令产生以下输出：

count    10.000000
mean     73.500000
std       5.778312
min      66.000000
25%      69.250000
50%      72.500000
75%      79.000000
max      82.000000
Name: Temperatures, dtype: float64

请注意，返回的信息是一个Series对象，因此可以存储以备后续计算需要。

要绘制系列的图表，我们使用plot()方法。如果我们只需要快速的图形数据概览，可以运行以下命令：

temps.plot()

然而，使用pandas也可以生成格式化良好的生产级数据图，因为matplotlib的所有功能都被pandas所支持。以下代码展示了如何使用在第三章，使用 matplotlib 绘图中讨论的某些图形格式化选项：

temps.plot(style='-s', lw=2, color='green')
axis((6,15,65, 85))
xlabel('Day')
ylabel('Temperature')
title('Maximum daily temperatures in June')
None # prevent text output

上述命令行生成如下图表：

假设我们要找出温度超过 75 度的日期。可以通过以下表达式实现：

temps[temps > 75]

上述命令返回以下系列：

7     76
11    80
12    82
15    80
Name: Temperatures, dtype: float64

Series类提供了许多有用的方法。记住，为了查看所有可用的方法，我们可以使用 IPython 的代码补全功能。开始输入temps.，你将看到可用的方法列表。

然后按下Tab键，一个列出所有可用方法的窗口将弹出。你可以探索其中的可用项。

DataFrame 类

DataFrame类用于表示二维数据。为了说明它的使用，下面我们创建一个包含学生数据的DataFrame类：

grades = DataFrame(
 [['Alice',  80., 92., 84,],
 ['Bob',    78., NaN, 86,],
 ['Samaly', 75., 78., 88.]],
 index = [17005, 17035, 17028],
 columns = ['Name', 'Test 1', 'Test 2', 'Final']
 )

这段代码展示了构建DataFrame类的最直接方式之一。在前面的例子中，数据可以指定为任何二维的 Python 数据结构，比如列表的列表（如示例所示）或NumPy数组。index选项设置行名，这里是表示学生 ID 的整数。同样，columns选项设置列名。index和column参数都可以作为任何一维的 Python 结构给出，例如列表、NumPy数组或Series对象。

要显示DataFrame类的输出，在单元格中运行以下语句：

grades

上述命令显示一个格式化良好的表格，如下所示：

DataFrame类具有极为灵活的初始化接口。我们建议读者运行以下命令来了解更多信息：

DataFrame?

这将显示有关构造选项的信息。我们的目标不是覆盖所有可能性，而是给出其灵活性的概念。请在单元格中运行以下代码：

idx = pd.Index(["First row", "Second row"])
col1 = Series([1, 2], index=idx)
col2 = Series([3, 4], index=idx)
data = {"Column 1":col1, "Column2":col2}
df = DataFrame(data)
df

上述代码生成如下表格：

这个例子展示了一种有用的方式来理解DataFrame对象：它由一个Series对象的字典组成，Index对象标记表格的行。字典中的每个元素对应表格中的一列。请记住，这仅仅是用来概念化DataFrame对象的一种方式，并不是描述其内部存储。

回到我们的学生数据示例。我们来添加一个列，表示每个学生的总分，计算方法是所有成绩的平均值，期末成绩的权重是两倍。可以使用以下代码进行计算：

grades.loc[:,'Score'] = 0.25 * (grades['Test 1'] + grades['Test 2'] + 2 * grades['Final']) 
grades

前述命令行的输出如下：

在前面的命令行中，我们使用了以下推荐的DataFrame类元素访问方法：

• .loc：该方法基于标签，即元素位置被解释为表格中的标签（行或列）。在前面的示例中使用了该方法。

• .iloc：该方法基于整数索引。参数必须是整数，并作为表格行和列的零基索引进行解释。例如，grades.iloc[0,1]表示前面示例中第 0 行和第 1 列的数据，即 Alice 在测试 1 中的成绩。

• .ix：这种索引方法支持混合整数和标签索引。例如，grades.ix[17035, 4]和grades.ix[17035, 'Score']都指 Bob 在该课程中的成绩。请注意，pandas 足够智能，能够识别行标签是整数，因此索引17035指的是标签，而不是表格中的位置。实际上，尝试访问grades.ix[1, 4]元素会引发错误，因为没有标签为 1 的行。

要使用这些方法，DataFrame对象中必须已经存在相应的条目（或条目）。因此，这些方法不能用于插入或附加新数据。

请注意，Bob 在第二次测试中没有成绩，这通过NaN条目表示（他可能在测试当天生病了）。当他进行补考时，可以通过以下方式更新他的成绩：

grades.loc[17035,'Test 2'] = 98
grades

在输出中，你会注意到 Bob 的期末成绩没有自动更新。这并不令人惊讶，因为DataFrame对象并不是作为电子表格程序设计的。要执行更新，你必须显式地再次执行计算成绩的单元格。这样做后，表格将如下所示：

提示

也可以使用常规索引来访问DataFrame条目，但这种做法不被推荐。例如，要引用 Samaly 的期末成绩，我们可以使用链式****引用，即使用grades['Test 2'][17028]。 （注意索引的顺序！）我们将避免这种用法。

老师有点失望，因为没有学生获得 A（分数超过 90）。于是，学生们被布置了一个额外的学分作业。为了在Final列旁边添加一个新成绩组成的列，我们可以运行以下命令：

grades.insert(4, 'Extra credit', [2., 6., 10.])
grades

显然，我们也可以插入行。为了添加一个新学生，我们可以使用以下命令：

grades.loc[17011,:] = ['George', 92, 88, 91, 9, NaN]
grades

当然，必须根据以下方式更新成绩，以考虑额外学分：

grades.loc[:,'Score'] = 0.25 * (grades['Test 1'] + grades['Test 2'] + 2 * grades['Final']) + grades['Extra credit'] 
grades

现在，假设我们想找到所有在 Test 1 中获得 A 并且分数低于 78 的学生。我们可以通过使用布尔表达式作为索引来实现，如下代码所示：

grades[(grades['Score'] >= 90) & (grades['Test 1'] < 78)]

从上面的示例中应注意两点重要内容：

• 我们需要使用&操作符，而不是and操作符。

• 由于&操作符的优先级较高，括号是必须的。

这将返回一个子表格，包含满足布尔表达式条件的行。

假设我们想要获得分数至少为 80 但不到 90 的学生的名字和分数（这些学生可能代表“B”级学生）。以下命令将对我们有所帮助：

grades[(80 <= grades['Score']) & grades['Score'] < 90].loc[:,['Name', 'Score']]]

这段代码的作用如下：

• 表达式grades[(80 <= grades['Score']) & grades['Score'] < 90]创建一个DataFrame类，包含所有学生的数据，这些学生的分数至少为 80，但小于 90。

• 然后，.loc[:,'Name', 'Score']获取这个DataFrame类的一个切片，其中包含所有行及标签为Name和Score的列。

关于 pandas 数据结构的一个重要点是，每当引用数据时，返回的对象可能是原始数据的副本或视图。让我们创建一个带有伪随机数据的DataFrame类来查看一些示例。为了让事情更有趣，每列将包含具有给定均值和标准差的正态数据。代码如下：

means = [0, 0, 1, 1, -1, -1, -2, -2]
sdevs = [1, 2, 1, 2,  1,  2,  1,  2]
random_data = {}
nrows = 30
for mean, sdev in zip(means, sdevs):
 label = 'Mean={}, sd={}'.format(mean, sdev)
 random_data[label] = normal(mean, sdev, nrows)
row_labels = ['Row {}'.format(i) for i in range(nrows)]
dframe = DataFrame (random_data, index=row_labels)

上述命令行创建了我们需要的数据，用于示例。请执行以下步骤：

1. 定义 Python 列表means和sdevs，分别包含分布的均值和标准差值。

2. 然后，创建一个名为random_data的字典，字典的键为字符串类型，且对应将要创建的DataFrame类的列标签。

3. 字典中的每个条目都对应一个大小为nrows的列表，包含由NumPy的normal()函数生成的数据。

4. 创建一个名为row_labels的列表，包含DataFrame类的行标签。

5. 使用这两组数据，即random_data字典和row_labels列表，在DataFrame构造函数中。

上述代码将生成一个 30 行 8 列的表格。你可以通过在单元格中单独评估dframe来查看这个表格。请注意，尽管这个表格的大小适中，但 IPython 笔记本在显示时做得非常好。

现在让我们选择DataFrame类的一个切片。为了演示的目的，我们将使用混合索引的.ix方法：

dframe_slice = dframe.ix['Row 3':'Row 11', 5:]
dframe_slice

请注意，范围是如何指定的：

• 表达式'Row 3':'Row 11'表示一个由标签指定的范围。请注意，与 Python 中的常见假设相反，这个范围包括最后一个元素（在本例中是Row 11）。

• 表达式5:（数字 5 后跟一个冒号）表示一个数值范围，从第五列到表格的末尾。

现在，在一个单元格中运行以下命令行：

dframe_slice.loc['Row 3','Mean=1, sd=2'] = normal(1, 2)
print dframe_slice.loc['Row 3','Mean=1, sd=2']
print dframe.loc['Row 3','Mean=1, sd=2']

第一行对数据表中的一个单元格进行了重采样，接下来的两行打印了结果。请注意，打印的值是相同的！这表明没有发生复制，变量dframe_slice引用了已经存在于由dframe变量引用的DataFrame类中的相同对象（内存区域）。（这类似于 C 语言中的指针，多个指针可以引用同一块内存。实际上，这是 Python 中变量的标准行为：默认没有复制。）

如果我们真的想要一个副本怎么办？所有 pandas 对象都有一个copy()方法，所以我们可以使用以下代码：

dframe_slice_copy = dframe.ix['Row 3':'Row 11', 5:].copy()
dframe_slice_copy

上面的命令行将产生与前一个示例相同的输出。然而，注意如果我们修改dframe_slice_copy会发生什么：

dframe_slice_copy.loc['Row 3','Mean=1, sd=2'] = normal(1, 2)
print dframe_slice_copy.loc['Row 3','Mean=1, sd=2']
print dframe.loc['Row 3','Mean=1, sd=2']

现在打印的值不同，确认只有副本被修改。

注意

在某些情况下，了解数据是在切片操作中被复制还是仅仅被引用是很重要的。特别是在处理更复杂的数据结构时，要格外小心。本书无法全面覆盖这一话题。然而，使用.loc、.iloc和.ix，如前面的例子所示，足以避免问题。关于链式索引可能引发错误的示例，请参考pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy以获取更多信息。

如果你遇到与SettingWithCopy相关的警告，检查一下是否尝试使用链式索引修改DataFrame对象的条目，比如在dframe_object['a_column']['a_row']中。将对象访问改为使用.loc，例如，将消除警告。

为了结束这一节，我们来看一些更多的切片DataFrame的示例。在以下所有示例中，都没有发生复制；只是创建了数据的新引用。

• 使用列表作为索引进行切片可以通过以下命令行执行：

  dframe.ix[['Row 12', 'Row 3', 'Row 24'], [3, 7]]
  
  

• 使用切片重新排序列可以通过以下命令行执行：

  dframe.iloc[:,[-1::-1]]
  
  

  上面的示例颠倒了列的顺序。若要进行任意的重新排序，可以使用包含列位置排列的列表：

  dframe.iloc[:,[2, 7, 0, 1, 3, 4, 6, 5]]
  
  

  请注意，dframe对象中的列没有实际重新排序，因为数据没有被复制。

• 使用布尔运算进行切片可以通过以下命令行执行：

  dframe.loc[dframe.loc[:,'Mean=1, sd=1']>0, 'Mean=1, sd=1']
  
  

  上述命令行选择标记为Mean=1, sd=1（且为正数）的列元素，并返回一个Series对象（因为数据是一维的）。如果你对这种方式理解有困难，可以单独在一个单元格中运行以下命令行：

  dframe.loc[:,'Mean=1, sd=1']>0
  
  

  这个语句将返回一个包含布尔值的Series对象。之前的命令行选择了dframe中对应于Series对象中为True的位置的行。

• 切片通常会返回一个形状与原始对象不同的对象。在需要保持形状不变的情况下，可以使用where()方法，如下所示。

  dframe.where(dframe>0)
  
  

  上述命令行返回一个DataFrame类，其条目对应于原始dframe对象中非负值的位置，并包含缺失值（NaN）。

  我们还可以使用以下命令行指定一个值，用来替换不满足给定条件的值：

  dframe.where(dframe>0, other=0)
  
  

  这条命令行将把对应于非负值的条目替换为 0。

计算和图形工具

pandas 对象具有丰富的内建计算工具。为了说明这些功能，我们将使用前一节中定义的dframe对象中存储的随机数据。如果你丢弃了该对象，这里是如何重新构造它的：

means = [0, 0, 1, 1, -1, -1, -2, -2]
sdevs = [1, 2, 1, 2,  1,  2,  1,  2]
random_data = {}
nrows = 30
for mean, sdev in zip(means, sdevs):
 label = 'Mean={}, sd={}'.format(mean, sdev)
 random_data[label] = normal(mean, sdev, nrows)
row_labels = ['Row {}'.format(i) for i in range(nrows)]
dframe = DataFrame (random_data, index=row_labels)

让我们探索一些内建计算工具的功能。

• 要获取对象的可用方法列表，可以在单元格中输入以下命令：

  dframe.
  
  

• 然后，按下Tab键。自动完成弹出窗口允许我们通过双击选择一个方法。例如，双击mean。单元格文本会变为以下内容：

  dframe.mean
  
  

• 现在，在上述命令行后加上问号并运行该单元格：

  dframe.mean?
  
  

  这将显示关于mean方法的信息（该方法显而易见地计算数据的平均值）。

使用 Tab 自动完成和 IPython 的帮助功能是学习 pandas 功能的一个极好的方法。我建议你总是通过这种方式显示文档，至少在第一次使用方法时要这样做。了解 pandas 提供的功能可以节省大量时间。

现在，让我们继续探索更多的功能：

• 假设我们想计算随机数据的列均值。这可以通过评估以下命令来实现：

  dframe.mean()
  
  

• 标准差值可以通过以下命令计算：

  dframe.std()
  
  

请注意，所有紧接着的命令行结果都作为Series对象返回，这是 pandas 为一维数据使用的默认对象类型。特别是，列标签成为对象的索引。假设我们想创建一个包含均值和标准差的DataFrame对象，且有两行。pandas 通过内建转换和构造器使得这一任务变得非常容易。

mean_series = dframe.mean()
std_series = dframe.std()
mean_std = DataFrame([dict(mean_series), 
 dict(std_series)], 
 index=['mean', 'std'])
mean_std

在这段代码中，我们首先计算均值和标准差，并将它们赋值给变量以便清晰表达。然后，我们调用接受 Python 字典列表的DataFrame构造函数。这是因为 pandas 允许将Series对象便捷地转换为字典：dict(mean_series)返回mean_series的字典表示，使用Series对象的索引作为字典的键。

假设我们想要标准化所有列的数据，使它们都具有相同的均值 100 和标准差 20。这可以通过以下命令行实现：

dframe_stnd = 100 + 20 * (dframe - mean_std.iloc[0,:]) / mean_std.iloc[1,:] 
dframe_stnd

上述命令行只是实现了标准化的定义：我们从数据中减去均值，除以标准差，按期望的偏差值进行缩放，然后加上期望的均值。为了检查我们是否得到预期结果，请运行以下命令行：

print dframe_stnd.mean()
print dframe_stnd.std()

为了说明这些可能性，让我们进行一个双侧检验，假设每列的均值为 0。我们首先计算列的Z 分数。每列的 Z 分数就是该列均值与模型均值（此处为 0）之间的偏差，并按标准差正确缩放：

zscores = mean_std.iloc[0,:] / (mean_std.iloc[1,:] / sqrt(len(dframe)))
zscores 

缩放因子sqrt(len(dframe))是数据点数量的平方根，数据点数量由表格中的行数给出。最后一步是计算每列的p 值。p 值只是衡量数据偏离均值超过相应 Z 分数的概率，这个概率基于假定的分布。这些值是从正态分布中获得的（严格来说，我们应该使用t 分布，因为我们使用的是样本标准差，但在这个示例中，这并不会产生实质性差异，因为数据是正态生成的，而且样本量足够大）。以下命令行使用 SciPy 中的正态分布对象norm来计算 p 值，并转换为百分比：

from scipy.stats import norm
pvalues = 2 * norm.cdf(-abs(zscores)) * 100
pvalues_series = Series(pvalues, index = zscores.index)
pvalues_series

计算 p 值的代码行如下：

pvalues = 2 * norm.cdf(-abs(zscores)) * 100

我们使用cdf()方法，它计算来自norm对象的正态分布的累积分布函数。然后我们将其乘以2，因为这是一个双侧检验，再乘以100以获得百分比。

下一行将 p 值转换为Series对象。这不是必需的，但能使结果更容易可视化。

以下是获得的结果：

Mean=-1, sd=1    1.374183e-02
Mean=-1, sd=2    1.541008e-01
Mean=-2, sd=1    2.812333e-26
Mean=-2, sd=2    1.323917e-04
Mean=0, sd=1     2.840077e+01
Mean=0, sd=2     6.402502e+01
Mean=1, sd=1     2.182986e-06
Mean=1, sd=2     5.678316e-01
dtype: float64

注意

请注意，在前面的示例中，您将得到不同的数字，因为数据是随机生成的。

由于数据生成方式的原因，结果正如我们所预期：除了均值为0的列之外，所有列的 p 值都非常小。

现在，让我们探索 pandas 提供的一些图形功能。pandas 绘图是使用 matplotlib 生成的，因此基本接口已经在第三章，使用 matplotlib 进行图形绘制中讨论过。在接下来的示例中，我们假设正在使用魔法。请在单元格中运行以下命令：

%pylab inline

pandas 大多数的绘图功能是作为 Series 或 DataFrame 对象的方法实现的。

让我们在表格中定义以下数据，以包括更多的数据点：

means = [0, 0, 1, 1, -1, -1, -2, -2]
sdevs = [1, 2, 1, 2,  1,  2,  1,  2]
random_data = {}
nrows = 300
for mean, sdev in zip(means, sdevs):
 label = 'Mean={}, sd={}'.format(mean, sdev)
 random_data[label] = normal(mean, sdev, nrows)
row_labels = ['Row {}'.format(i) for i in range(nrows)]
dframe = DataFrame (random_data, index=row_labels)

要显示数据的直方图网格，我们可以使用以下命令：

dframe.hist(color='DarkCyan')
subplots_adjust(left=0.5, right=2, top=2.5, bottom=1.0)

我们使用 hist() 方法生成直方图，并使用 color 选项，该选项传递给实际绘图的 matplotlib 函数调用。代码的第二行通过在图表中添加空隙，避免了坐标轴标签重叠。你可能会发现有些直方图看起来不正常。要修复它们的外观，可以调整 hist() 方法的 bins 和 range 选项，如下所示：

dframe.loc[:,'Mean=0, sd=2'].hist(bins=40, range=(-10,10), color='LightYellow')
title('Normal variates, mean 0, standard deviation 2')

这将绘制一个列中数据的直方图，均值为 0，标准差为 2，在从 -10 到 10 的范围内有 40 个区间。换句话说，每个区间的宽度是 0.5。请注意，绘图可能不包括从 -10 到 10 的整个范围，因为 pandas 会将绘图限制在实际包含数据的范围内。

例如，让我们根据几何布朗运动（GBM）生成数据，这是一种在数学金融中用于表示股票价格演变的模型。（详情见 en.wikipedia.org/wiki/Geometric_Brownian_motion。）这个模型通过两个参数定义，表示股票的百分比漂移和百分比波动率。我们从在模型中定义这两个值开始，并定义股票的初始值：

mu = 0.15
sigma = 0.33
S0 = 150

模拟应该从时间 0.0 运行到最大时间 20.0，并且我们希望生成 200 个数据点。以下命令定义了这些参数：

nsteps = 200
tmax = 20.
dt = tmax/nsteps
times = arange(0, tmax, dt)

股票模型自然会通过时间序列（Series 对象）来表示。然而，为了简化模拟，我们将使用 DataFrame 对象，并逐列构建模拟。我们从一个非常简单的表格开始，该表格仅包含整数索引和模拟时间：

gbm_data = DataFrame(times, columns=['t'], index=range(nsteps))

要查看表格的前几行，我们可以使用以下命令：

gbm_data.loc[:5,:]

在每一列添加后，你可能想要运行这个命令，以便更好地了解模拟的进展。

GBM 模型的基础是（毫不奇怪）一种叫做布朗运动（BM）的随机过程。这个过程包含两个部分。一个确定性成分，叫做漂移，其计算方法如下：

gbm_data['drift'] = (mu - sigma**2/2) * gbm_data.loc[:,'t']

下一个组件添加了随机性。它通过增量来定义，这些增量是正态分布的，均值为零，标准差由时间间隔乘以百分比波动性给出：

gbm_data['dW'] = normal(0.0, sigma * dt, nsteps)

然后，BM 组件被定义为增量的累计和，如下命令行所示：

gbm_data['W'] = gbm_data.loc[:,'dW'].cumsum()
gbm_data.ix[0, 'W'] = 0.0

在前面的命令行中，我们添加了第二行，因为我们希望从0开始，这不是cumsum()方法采用的默认设置。

我们现在准备计算股票模拟。计算方法是取漂移成分，加上 BM 成分，然后对结果取指数，最后将其与股票的初始值相乘。所有这些操作可以通过以下命令完成：

gbm_data['S'] = S0 * exp(gbm_data.loc[:,'drift'] + gbm_data.loc[:,'W'])

我们现在准备使用以下命令行绘制模拟结果：

gbm_data.plot(x='t', y='S', lw=2, color='green',
 title='Geometric Brownian Motion')

前面的命令行生成了以下图表。显然，由于随机性，您得到的图表会有所不同。

使用现实数据集的示例

在这一部分，我们将使用一个中等规模的现实数据集。我们将使用世界发展指标数据集，该数据集由世界银行免费提供。这个数据集规模适中，不算太大或复杂，适合用来实验。

在任何实际应用中，我们都需要从某个来源读取数据，将其重新格式化以满足我们的需求，然后将重新格式化的数据保存回某个存储系统。pandas 提供了多种格式的数据检索和存储功能：

• 逗号分隔 值（CSV）在文本文件中

• Excel

• JSON

• SQL

• HTML

• Stata

• 文本格式的剪贴板数据

• Python 序列化数据

pandas 支持的格式列表随着每次库的更新而不断增长。请参考pandas.pydata.org/pandas-docs/stable/io.html获取当前的列表。

由于本书的范围限制，无法覆盖 pandas 支持的所有格式。我们将把示例限制为 CSV 文件，这是一种简单的文本格式，广泛使用。大多数软件包和数据源都可以将数据格式化为 CSV 文件。

有趣的是，CSV 并不是一个正式描述的存储格式。pandas 在提供足够的选项来读取大多数文件方面做得很好。然而，数据的格式可能会因数据源而异。幸运的是，由于 CSV 文件只是文本文件，我们可以在电子表格程序或文本编辑器中打开文件，查看其结构。

本节的数据集可以从data.worldbank.org/data-catalog/world-development-indicators下载，也可以在本书网站上获取。如果选择从原始网站下载数据，请确保选择 CSV 文件格式。该文件为压缩的 ZIP 格式，大小约为 40 MB。解压缩后，我们将获得以下文件。

• WDI_Country.csv

• WDI_CS_Notes.csv

• WDI_Data.csv

• WDI_Description.csv

• WDI_Footnotes.csv

• WDI_Series.csv

• WDI_ST_Notes.csv

如同任何现实数据集一样，总是有很多与数据相关的附加信息。这些信息称为元数据，用于提供关于数据集的信息，包括行和/或列使用的标签、数据收集的细节以及数据含义的解释。元数据包含在归档中的各个文件中。建议读者使用电子表格软件（或文本编辑器）打开不同的文件，以便了解可用的信息类型。对我们来说，最重要的元数据文件是WDI_Series.csv，它包含有关数据标签含义的信息，涉及数据中包含的多个时间序列。

实际的数据存储在WDI_Data.csv文件中。由于该文件包含一些元数据，我们可以仅使用该文件完成所有工作。

确保WDI_Data.csv文件与包含 IPython 笔记本文件的目录在同一目录中，并在一个单元格中运行以下命令：

wdi = pd.read_csv('WDI_Data.csv')

这将读取文件并将其存储在一个DataFrame对象中，我们将该对象赋值给变量wdi。文件的第一行默认假设包含列标签。我们可以通过运行以下命令来查看表格的开头：

wdi.loc[:5]

请注意，DataFrame类默认是按整数索引的。可以通过向read_csv()方法传递index_col参数来选择数据文件中的某一列作为索引。索引列可以通过其位置或标签在文件中指定。read_csv()的许多选项在pandas.pydata.org/pandas-docs/stable/io.html#io-read-csv-table中有详细讨论。

对文件进行检查后发现，需要对其进行一些处理，才能将其转换为易于使用的格式。文件的每一行包含一个国家和经济指标的年度数据时间序列。一个初步步骤是获取文件中包含的所有国家和经济指标。要获取唯一国家名称的列表，可以使用以下命令：

countries = wdi.loc[:,'Country Name'].unique()

若要查看代表了多少个国家，请运行以下命令：

len(countries)

文件中的一些条目实际上对应的是国家组，如撒哈拉以南非洲。

现在，对于指标，我们可以运行以下命令行：

indicators = wdi.loc[:,'Indicator Code'].unique()

文件中有超过 1300 个不同的经济指标。可以通过运行以下命令来验证：

len(indicators)

为了展示可能感兴趣的不同计算方法，我们以一个单一国家为例，比如巴西。假设我们只关心国内生产总值（GDP）信息。现在，我们将看到如何从表格中选择我们感兴趣的数据。为了简化示例，我们将分两步进行选择。首先，我们使用以下命令选择所有Brazil的行：

wdi_br = wdi.loc[wdi.loc[:,'Country Name']=='Brazil',:]

在前面的命令行代码中，考虑以下表达式：

wdi.loc[:,'Country Name']=='Brazil'

这会选择所有国家名称字符串等于Brazil的行。对于这些行，我们想选择表格的所有列，正如切片操作中的第一个术语中的冒号所示。

现在，让我们选择所有与 GDP 数据相关的行。我们首先定义一个函数，给定一个字符串，判断它是否包含子字符串GDP（忽略大小写）：

select_fcn = lambda string: string.upper().find('GDP') >= 0

我们现在想要在wdi_br中选择那些当select_fcn应用于Indicator Code列时返回True的行。这可以通过以下命令行完成：

criterion = wdi_br.loc[:,'Indicator Code'].map(select_fcn)
wdi_br_gdp = wdi_br.loc[criterion,:]

Series对象的map()方法正是我们需要的功能：它会将函数应用于序列的所有元素。我们将此调用的结果赋值给变量criterion。然后，我们在定义wdi_br_gdp的切片操作中使用criterion。要查看选中了多少行，可以运行以下命令：

len(wdi_br_gdp)

在编写本书时使用的数据集中，前面的命令返回32。这意味着名为Brazil的国家有 32 个与 GDP 相关的指标。由于数据量现在变得可控，我们可以使用以下命令显示一个包含指标代码及其含义的表格：

wdi_br_gdp.loc[:,['Indicator Code', 'Indicator Name']]

上述命令行生成了一个格式良好的表格，包含指标及其对应名称，如下表所示：

假设我们只对四个指标感兴趣：GDP、年 GDP 增长率、人均 GDP 和人均 GDP 增长率。我们可以通过以下命令进一步筛选数据：

wdi_br_gdp = wdi_br_gdp.loc[[37858, 37860, 37864, 37865], :]

这将生成一个相对易于处理的表格，包含 4 行 58 列。每行都包含从 1960 年开始的相应 GDP 数据的时间序列。

当前表格布局的问题在于，它是“转置”的，和 pandas 中常见的布局惯例相反：时间序列横跨表格的行，而不是纵向排列在列中。因此，我们还需要对表格做一些调整。我们希望表格的索引是年份。我们还希望为每个经济指标设置一列，并使用经济指标名称（而不是代码）作为列的标签。以下是如何实现这一点：

idx = wdi_br_gdp.loc[:,'1960':].columns
cols = wdi_br_gdp.loc[:,'Indicator Name']
data = wdi_br_gdp.loc[:,'1960':].as_matrix()
br_data = DataFrame(data.transpose(), columns=cols, index=idx)

以下是对上述命令行功能的解释：

1. 我们首先使用 DataFrame 对象的 columns 字段定义一个与表格中年份对应的 Index 对象。该对象存储在变量 idx 中。

2. 然后，我们创建一个包含列名的对象。这是一个 Series 对象，存储在变量 cols 中。

3. 接下来，我们提取感兴趣的数据，即表格中对应 1960 年之后年份的部分数据。我们使用 DataFrame 对象的 as_matrix() 方法将数据转换为 NumPy 数组，并将其存储在变量 data 中。

4. 最后，我们调用 DataFrame 构造函数来创建新表格。

现在我们已经获得了格式良好的数据，是时候保存它了：

br_data.to_csv('WDI_Brazil_GDP.csv')

此时，我们可以在电子表格程序中打开 WDI_Brazil_GDP.csv 文件进行查看。

现在，让我们开始通过创建一些图表来操作数据。首先绘制自 1980 年起的 GDP 和 GDP 增长图表。由于数据是以美元表示的，我们进行缩放以使其以十亿美元为单位。

pdata = br_data.ix['1970':, 0] / 1E9
pdata.plot(color='DarkRed', lw=2, 
 title='Brazil GDP, billions of current US$')

上述命令行将生成以下图表：

作为最后一个示例，让我们绘制一张比较 BRIC（巴西、俄罗斯、印度和中国）国家在 2000 到 2010 年间人均 GDP 增长百分比的图表。由于我们已经探索了数据的结构，这项任务相对简单：

bric_countries = ['Brazil', 'China', 'India', 'Russian Federation']
gdp_code = 'NY.GDP.PCAP.KD.ZG'
selection_fct = lambda s: s in bric_countries
criterion = wdi.loc[:,'Country Name'].map(selection_fct)
wdi_bric = wdi.loc[criterion & (wdi.loc[:,'Indicator Code'] == gdp_code),:]

我们首先定义一个包含金砖国家名称的列表，并定义一个表示人均 GDP 增长百分比的指标代码字符串。接着，我们定义一个选择函数：如果字符串是金砖国家名称之一，则选择该字符串。然后，使用 map() 方法将选择函数应用于 Country Name 列的所有条目。最后一行命令执行实际的选择操作。请注意使用布尔运算符 & 来结合行选择中使用的两个条件。

接下来，我们对数据进行重新格式化，将相关数据系列排列在表格的列中。命令行与前一个示例中的命令相似：

df_temp = wdi_bric.loc[:, '2000':'2010']
idx = df_temp.columns
cols = wdi_bric.loc[:, 'Country Name']
data = df_temp.as_matrix()
bric_gdp = DataFrame(data.transpose(), columns=cols, index=idx)

完成后，绘制数据就变得非常简单：

bric_gdp.plot(lw=2.5,
 title='Annual per capita GDP growth (%)')

上述命令行将生成以下图表：

概述

在本章中，我们介绍了 pandas 的对象 Series 和 DataFrame，它们是用于数据导向计算的专用容器。我们讨论了如何创建、访问和修改这些对象，包括高级索引和切片操作。我们还探讨了 pandas 提供的计算和图形能力。然后，我们讨论了如何利用这些功能处理真实数据集。

在下一章中，我们将学习如何使用 SciPy 解决建模、科学和工程中的高级数学问题。

第五章。使用 SciPy、Numba 和 NumbaPro 的高级计算

本章中，用户将学习如何使用SciPy执行科学计算。然后，将介绍Numba包作为加速计算的方法。最后，将介绍NumbaPro的 GPU 并行执行能力。

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

• SciPy概述

• 使用SciPy的高级数学算法

• 使用Numba和NumbaPro加速计算

在运行本章中的示例之前，请通过在计算单元中运行以下命令来加载pylab：

%pylab inline

SciPy概述

SciPy是一个广泛应用于数学和科学计算的库。以下是库中所有可用模块的完整列表：

模块	功能
cluster	聚类算法
constants	物理和数学常数
fftpack	快速傅里叶变换
integrate	积分与常微分方程
interpolate	插值与样条
io	输入与输出
linalg	线性代数
ndimage	图像处理
odr	正交距离回归
optimize	优化与根求解
signal	信号处理
sparse	稀疏矩阵
spatial	空间数据结构
special	特殊函数
stats	统计分布
weave	C/C++ 集成

在脚本中导入SciPy模块的标准方式是使用以下命令行：

from scipy import signal

然后，可以使用通常的模块引用语法来调用各个函数，如下所示：

signal.correlate(…)

然而，许多最常用的函数可以在SciPy层次结构的顶层找到。此外，我们使用 IPython 的交互模式，并使用（本书中假设的）魔法命令，如下所示：

%pylab inline

许多函数将无需显式引用模块即可使用。

在下一节中，我们将展示SciPy中可用的一些函数示例。读者不需要深入了解示例中将使用的数学技术和算法。

使用 SciPy 的高级数学算法

在本节中，我们将涵盖SciPy中一些可用的算法。以下每个子节都包含一个来自应用科学重要领域的代表性示例。这些示例的选择不需要广泛的领域知识，但仍然具有现实性。我们呈现的主题和示例如下：

• 解方程和寻找最优值：我们将研究一个市场模型，它需要求解一个非线性系统，以及一个需要非标准优化的设施选址问题。

• 微积分和常微分方程：我们将展示一个利用积分微积分进行的体积计算，并介绍牛顿的炮，这是由艾萨克·牛顿提出的一个思想实验，我们将用常微分方程系统对其建模。最后，我们将介绍一个三维系统——著名的洛伦兹方程，这是早期表现出混沌行为的例子。

求解方程和寻找最优值

为了说明这个主题，我们使用了经济学中一个标准的供需模型。在这个模型中，供需与价格之间有函数关系，市场均衡通过确定供需曲线的交点来找到。我们在示例中使用的数学公式有些是任意的（因此可能不现实），但这些公式超出了教科书中通常假设供需关系为线性关系的范围。

定义供需曲线的公式如下：

我们将使用函数工厂模式。在一个单元格中运行以下代码：

def make_supply(A, B, C):
 def supply_func(q):
 return A * q / (C  - B * q)
 return supply_func
def make_demand(K, L):
 def demand_func(q):
 return K / (1 + L * q)
 return demand_func

上述代码并没有直接定义供需曲线，而是指定了函数工厂。这样做可以更方便地处理参数，这在应用问题中非常常见，因为我们希望相同的模型能适用于多种情况。

接下来，我们设置参数值，并调用函数工厂来定义实际评估供需曲线的函数，如下所示：

A, B, C = 23.3, 9.2, 82.4
K, L = 1.2, 0.54
supply = make_supply(A, B, C)
demand = make_demand(K, L)

以下代码行生成了曲线的图像：

q = linspace(0.01,5,200)
plot(q, supply(q), lw = 2)
plot(q, demand(q), lw = 2)
title('Supply and demand curves')
xlabel('Quantity (thousands of units)')
ylabel('Price ($)')
legend(['Supply', 'Demand'], loc='upper left')

以下是上述代码行输出的图像：

选择的供需曲线反映了合理的假设：随着价格上涨，供应增加而需求减少。即使价格为零，需求也是有限的（反映了市场中对产品感兴趣的人口有限）。另一方面，供给曲线具有垂直渐近线（在图中未显示），表明存在生产限制（即使价格趋近于无限，也只能提供有限的市场供应量）。

市场的均衡点是供需曲线的交点。为了找到均衡点，我们使用optimize模块，该模块除了提供优化函数外，还具有求解数值方程的功能。推荐用来求解一元函数的函数是brentq()，如下代码所示：

from scipy import optimize
def opt_func(q):
 return supply(q) - demand(q)
q_eq = optimize.brentq(opt_func, 1.0, 2.0)
print q_eq, supply(q_eq), demand(q_eq)

brentq()函数假设我们要求解的方程的右边是0。因此，我们首先定义opt_func()函数，该函数计算供给和需求之间的差值。这个函数是brentq()的第一个参数。接下来的两个数值参数给出了包含解的区间。选择一个确切包含方程解的区间非常重要。在我们的例子中，通过查看图形很容易做到这一点，因为可以明显看到曲线在 1 和 2 之间交叉。

运行前面的代码会产生以下输出：

1.75322153719 0.616415252177 0.616415252177

第一个值是均衡点，即在最优价格下能够销售的单位数（以千为单位）。最优价格是通过供给曲线和需求曲线共同计算得出的（以检查这些值是否确实相同）。

为了说明一个有两个变量的优化问题，我们考虑一个最优设施位置问题。假设一座工厂有若干个制造站点，需要从一个单一的供给站点分配材料。工厂的车间是矩形的，分配轨道必须与工厂的墙壁平行。这最后一个要求使问题变得有趣。需要最小化的函数与所谓的出租车****距离相关，具体如以下图所示：

第一步是定义给定制造站点的位置，如下所示：

points = array([[10.3,15.4],[6.5,8.8],[15.6,10.3],[4.7,12.8]])

这些位置存储为一个 4 x 2 的NumPy数组，名为points，每行代表一个点。以下命令会绘制之前命令行中提到的这些点的图：

plot(points[:,0],points[:,1], 'o', ms=12, mfc='LightSkyBlue')

这些点使用圆形标记显示，标记由参数o指定，该参数还关闭了连接这些点的线段。ms和mfc选项分别指定标记的大小（以像素为单位）和颜色。生成的输出为以下图像：

下一步是定义要最小化的函数。我们再次倾向于定义一个函数工厂的方法，如下所示：

def make_txmin(points):
 def txmin_func(p):
 return sum(abs(points - p))
 return txmin_func

这段代码的关键在于计算出租车距离的方式，充分利用了NumPy数组操作的灵活性。此操作在以下代码行中完成：

return sum(abs(points - p))

这段代码首先计算了向量差值points-p。请注意，points是一个 4 x 2 的数组，而p是一个 1 x 2 的数组。NumPy意识到数组的维度不同，并利用其广播规则。结果是，数组p被从points数组的每一行中减去，这正是我们想要的。接着，abs()函数计算结果数组中每个条目的绝对值，最后sum()函数将所有条目相加。这一切都在一行代码中完成，工作量非常大！

然后我们需要使用函数工厂来定义实际计算出租车距离的函数。

txmin = make_txmin(points)

函数工厂只是简单地用包含实际位置的数组作为参数进行调用。此时，问题已经完全设置好，我们准备计算最优解，以下代码可以完成这一过程：

from scipy import optimize
x0 = array([0.,0.])
res = optimize.minimize(
 txmin, x0,
 method='nelder-mead',
 options={'xtol':1e-5, 'disp':True})

最小化通过调用minimize()函数来计算。该函数的前两个参数是之前单元格中定义的目标函数txmin()和初始猜测值x0。我们选择原点作为初始猜测值，但在实际问题中，我们会利用所收集的信息来选择一个接近实际最小值的猜测值。该函数提供了几种优化方法，适用于不同类型的目标函数。我们使用Nelder-Mead 方法，这是一种启发式算法，不要求目标函数具有光滑性，非常适合当前问题。最后，我们为该方法指定了两个选项：所需的容差和用于打印诊断信息的显示选项。这将生成以下输出：

Optimization terminated successfully.
 Current function value: 23.800000
 Iterations: 87
 Function evaluations: 194

上述输出表示已成功找到最小值，并给出了其值。请注意，与任何数值优化方法一样，通常只能保证找到局部最小值。在这种情况下，由于目标函数是凸的，因此最小值必定是全局最小值。该函数的结果存储在SciPy数据结构中，这种数据结构是optimize模块中定义的OptimizeResult类型。要获取设施的最优位置，我们可以使用以下命令：

print res.x

上述命令的输出如下：

[  8.37782286  11.36247412]

为了完成这个示例，我们展示了显示最优解的代码：

plot(points[:,0],points[:,1], 'o', ms=12, mfc='LightSkyBlue')
plot(res.x[0], res.x[1],'o', ms=12, mfc='LightYellow')
locstr = 'x={:5.2f}, y={:5.2f}'.format(res.x[0], res.x[1])
title('Optimal facility location: {}'.format(locstr))

对plot()函数的调用类似于之前示例中的调用。为了给图表添加格式化良好的标题，我们首先定义locstr字符串，用于显示最优位置坐标。这个字符串是一个 Python 格式化字符串，格式规范为{:5.2f}，即宽度为5且精度为2位的小数。结果如下图所示：

微积分和微分方程

作为微积分计算的一个例子，我们将展示如何计算旋转体的体积。该旋转体是通过将下图中显示的曲线绕y轴旋转得到的：

这条曲线是通过以下代码绘制的：

def make_gen(a,b):
 def gen_func(y):
 return a/pi * arccos(2.0 * y / b - 1.0)
 return gen_func
a = 5.
b = 4.
gen = make_gen(a,b)
x = linspace(0,b,200)
y = gen(x)
subplot(111, aspect='equal')
plot(x,y,lw=2)

该曲线本质上是一个拉伸和转置的反余弦函数，定义在make_gen()函数中。它依赖于两个参数，a和b，分别指定其高度和长度。make_gen()函数是一个函数工厂，返回一个实际计算曲线值的函数。实际定义曲线的函数称为gen()（代表生成器），因此这是要绘制的函数。

当这条曲线围绕垂直轴旋转时，我们得到了如下绘制的实体：

当然，前面的图是使用 IPython 生成的，使用了以下代码：

from mpl_toolkits.mplot3d import Axes3D
na = 50
nr = 50
avalues = linspace(0, 2*pi, na, endpoint=False)
rvalues = linspace(b/nr, b, nr)
avalues_r = repeat(avalues[...,newaxis], nr, axis=1)
xvalues = append(0, (rvalues*cos(avalues_r)).flatten())
yvalues = append(0, (rvalues*sin(avalues_r)).flatten())
zvalues = gen(sqrt(xvalues*xvalues+yvalues*yvalues))
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(xvalues, yvalues, zvalues, 
 color='Cyan',alpha=0.65,linewidth=0.)

该代码中的关键函数是最后一行中对plot_trisurf()的调用。该函数接受三个NumPy数组，xvalues、yvalues和zvalues，指定了表面上点的坐标。数组xvalues和yvalues定义了一系列同心圆上的点，如下图所示：

z坐标的值是通过在每个点计算gen(sqrt(x*x+y*y))得到的，这样就将同心圆上的所有点在 3D 图中分配了相同的高度。

要计算实体的体积，我们使用圆柱壳法。该方法的工作原理超出了本书的范围，但归根结底是计算一个积分，如下公式所示：

在这个公式中，f(x)函数代表围绕y轴旋转的曲线。为了计算这个积分，我们使用scipy.integrate包。我们使用quad()函数，适用于没有奇点的函数的通用积分。以下是这个公式的代码：

import scipy.integrate as integrate
int_func  = lambda x: 2 * pi * x * gen(x)
integrate.quad(int_func, 0, b)

导入integrate模块后，我们定义要进行积分的函数。请注意，由于这是一行计算，我们使用了lambda语法。最后，我们调用quad()来执行积分。调用的参数是被积函数和积分的边界（在这种情况下是从0到b）。以下是前面几行代码的输出：

(94.24777961000055, 1.440860870616234e-07)

第一个数字是积分值，第二个是误差估计。

在下一个示例中，我们考虑牛顿的大炮，这是现代物理和微积分根源的一个思维实验。以下图像展示了这种情况，这是艾萨克·牛顿爵士的著作《世界体系论》中的一幅雕刻：

牛顿让我们想象有一门大炮坐落在一座非常高的山顶上。如果大炮发射一颗抛射体，它会飞行一段时间，最终落到地面。抛射体的初速度越大，它落地的距离越远。让我们假设我们可以任意快地发射抛射体，并且没有空气阻力。这样，随着初速度的增加，最终抛射体将绕地球飞行，如果大炮足够迅速地移除，抛射体将永远继续绕地球轨道运行。牛顿用这个例子来解释月球是如何绕地球运行的，而不单单依靠重力作用下掉落。

要建模这个情况，我们需要使用牛顿的引力定律作为一个微分方程组：

我们不会尝试解释这些公式是如何得出的，唯一对我们重要的是有四个状态变量，前两个表示抛射体的位置，后两个表示其速度向量。由于运动发生在通过地球中心的平面内，因此只需要两个位置变量。G和M分别是表示牛顿万有引力常数和地球质量的常量。抛射体的质量不出现，因为引力质量与惯性质量恰好抵消。

使用SciPy来解决这个问题的第一步是定义这一组微分方程，具体实现如下代码：

M = 5.9726E24
G = 6.67384E-11
C = G * M
def ode_func(xvec, t):
 x1, x2, v1, v2 = xvec
 d = (x1 * x1 + x2 * x2) ** 1.5
 return array([v1, v2, -C * x1 / d, -C * x2 / d ])

我们只需要做的是定义一个函数，计算微分方程组的右侧。我们首先定义常量M和G（使用国际单位制），以及辅助常量C，因为G和M仅通过它们的乘积出现在方程中。系统由ode_func()函数表示。这个函数必须至少接受两个参数：一个NumPy数组xvec和一个浮动值t。在我们的例子中，xvec是一个四维向量，因为我们的系统有四个状态变量。t变量在系统中没有被使用，因为没有外力（如果我们发射的是火箭而不是抛射体，那么就会有外力）。然而，它仍然必须作为输入参数列出。

在ode_func()内部，我们首先通过赋值提取xvec向量的各个元素，如下所示：

x1, x2, v1, v2 = xvec

这不是严格必要的，但有助于提高可读性。然后我们计算辅助量d（这是最后两个方程的分母）。最后，输出数组根据系统中的公式进行计算。请注意，没有计算任何导数，因为求解器所需的所有信息都包含在方程的右侧。

我们现在准备使用以下代码行来求解微分方程组：

import scipy.integrate as integrate
earth_radius = 6.371E6
v0 = 8500
h0 = 5E5
ic = array([0.0, earth_radius + h0, v0, 0.0])
tmax = 8700.0
dt = 10.0
tvalues = arange(0.0, tmax, dt)
xsol = integrate.odeint(ode_func, ic, tvalues)

上述代码的第一行导入了integrate模块，该模块用于求解微分方程。然后，我们需要指定抛射物的初始位置和速度。我们假设火炮位于北极，位于 50,000 米高的塔楼上（虽然这显然不现实，我们选择如此大的值是为了增强轨道的可见性）。由于地球不是完美的球体，我们使用半径的平均值。初始速度设置为8500米/秒。

初始条件存储在一个NumPy数组中，通过以下赋值完成：

ic = array([0.0, earth_radius + h0, v0, 0.0])

下一步是定义初始时间（在我们的例子中为零）和解所需的时间数组。这通过以下三行代码完成：

tmax = 8650.0
dt = 60.0
tvalues = arange(0.0, tmax, dt)

我们首先定义tmax为模拟的持续时间（以秒为单位）。变量dt存储我们希望记录解的时间间隔。在上述代码中，解将在 8,650 秒内每 60 秒记录一次。最终时间通过反复试验选择，以大致对应抛射物的一个轨道。

我们现在准备计算解，这通过调用odeint()函数来完成。解存储在向量xsol中，该向量为每个计算解的时间提供一行。要查看该向量的前几行，我们可以运行以下命令：

xsol[:10]

上述命令产生了以下输出：

array([[  0.00000000e+00,   6.87100000e+06,   8.50000000e+03,
 0.00000000e+00],
 [  5.09624253e+05,   6.85581217e+06,   8.48122162e+03,
 -5.05935282e+02],
 [  1.01700042e+06,   6.81036510e+06,   8.42515172e+03,
 -1.00800330e+03],
 [  1.51991202e+06,   6.73500470e+06,   8.33257580e+03,
 -1.50243025e+03],
 [  2.01620463e+06,   6.63029830e+06,   8.20477026e+03,
 -1.98562415e+03],
 [  2.50381401e+06,   6.49702131e+06,   8.04345585e+03,
 -2.45425372e+03],
 [  2.98079103e+06,   6.33613950e+06,   7.85073707e+03,
 -2.90531389e+03],
 [  3.44532256e+06,   6.14878788e+06,   7.62903174e+03,
 -3.33617549e+03],
 [  3.89574810e+06,   5.93624710e+06,   7.38099497e+03,
 -3.74461812e+03],
 [  4.33057158e+06,   5.69991830e+06,   7.10944180e+03,
 -4.12884566e+03]])

这些值是从0秒到360秒的抛射物的位置和速度向量，时间间隔为60秒。

我们确实希望绘制轨道图。这可以通过在一个单元格中运行以下代码来完成：

subplot(111, aspect='equal')
axis(earth_radius * array([-1.5, 1.5, -1.8, 1.2]))
earth = Circle((0.,0.), 
 earth_radius, 
 ec='Black', fc='Brown', lw=3)
gca().add_artist(earth)
plot(xsol[:,0], xsol[:,1], lw=2, color='DarkBlue')
title('Newton\'s Canon, $v_0={}$ m/s'.format(v0))

我们希望两个轴使用相同的尺度，因为两个轴都表示以米为单位的空间坐标。这在第一行代码中完成。第二行设置轴的限制，使轨道的图表能够舒适地适应图像。

然后，我们使用以下代码绘制一个圆形来表示地球：

earth = Circle((0.,0.), 
 earth_radius, 
 ec='Black', fc='Brown', lw=3)
gca().add_artist(earth)

我们在图表中并未强调使用Artist对象，因为这些对象的级别较低，通常在科学绘图中不需要。这里，我们通过指定圆心、半径和外观选项来构建一个Circle对象：黑色边缘颜色、棕色填充颜色和线宽为 3。第二行代码展示了如何将Circle对象添加到图表中。

绘制地球后，我们使用以下代码绘制轨道：

plot(xsol[:,0], xsol[:,1], lw=2, color='DarkBlue')

这是对plot()函数的标准调用。请注意，我们只绘制xsol数组的前两列，因为这表示抛射物的位置（回想一下，其他两列表示速度）。以下图像是我们得到的输出：

微分方程的数值解是一个复杂的话题，完整的处理超出了本书的范围，但我们将呈现odeint()函数的完整形式并评论其中的一些选项。odeint()函数是 Python 对ODEPACK的lsoda求解器的封装，ODEPACK是一个 Fortran 库。有关求解器的详细信息，请访问people.sc.fsu.edu/~jburkardt/f77_src/odepack/odepack.html

以下代码行是odeint()的完整签名：

odeint(ode_func, x0, tvalues, args=(), Dfun=None, col_deriv=0,
 full_output=0, ml=None, mu=None, rtol=None, atol=None,
 tcrit=None, h0=0.0, hmax=0.0, hmin=0.0,ixpr=0, mxstep=0,
 mxhnil=0, mxordn=12, mxords=5, printmessg=0)

参数ode_func、x0和tvalues已经讨论过。参数args允许我们向所求解的方程传递额外的参数。这是一种非常常见的情况，下面的示例将对此进行说明。在这种情况下，定义系统的函数必须具有以下签名：

ode_func(x, t, p1, p2, … , pn)

这里，p1、p2和pn是额外的参数。这些参数对于单个解是固定的，但可以在不同的解之间发生变化（通常用于表示环境）。传递给args的元组长度必须与ode_func()所需的参数数量完全相等。

以下是常见选项意义的部分列表：

• Dfun是一个计算系统雅可比矩阵的函数。这可能提高解的精度。

• col_deriv指定了雅可比矩阵的导数是沿着列（True，更快）还是行（False）排列。

• 如果full_output设置为True，输出将包含关于解算过程的诊断信息。如果错误积累并且解算过程未成功完成，这可能会很有用。

在本节的最后一个示例中，我们展示了洛伦兹振荡器，这是大气对流的简化模型，也是一个在某些参数值下表现出混沌行为的著名方程。我们还将利用这个示例演示如何在三维空间中绘制解。

洛伦兹系统由以下方程定义：

我们通过定义一个表示该系统的 Python 函数开始，如下所示：

def ode_func(xvec, t, sigma, rho, beta):
 x, y, z = xvec
 return array([sigma * (y - x),
 x * (rho - z) - y,
 x * y - beta * z ])

这个系统与前一个系统的唯一区别在于有了sigma、rho和beta参数。注意，它们只是作为额外的参数添加到ode_func()中。求解这个方程几乎与求解前一个示例一样：

tmax = 50
tdelta = 0.005
tvalues = arange(0, tmax, tdelta) 
ic = array([0.0, 1.0, 0.0])
sol = integrate.odeint(ode_func, ic, tvalues, 
 args=(10., 28., 8./3.))

我们定义了时间数组和初始条件，就像在之前的示例中做的那样。请注意，由于这是一个三维问题，因此初始条件是一个包含三个分量的数组。接下来是对odeint()的调用。此时调用增加了一个额外的参数：

args=(10., 28., 8./3.)

这将sigma、rho和beta分别设置为10、28和8/3。这些值已知对应于混沌解。

然后可以通过以下代码绘制解决方案：

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(8,8))
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')
x, y, z = sol.transpose() 
ax.plot(x, y, z, lw=0.5, color='DarkBlue')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')

前三行代码设置了三维绘图的坐标轴。下一行提取了适合绘图的数据格式：

x, y, z = sol.transpose()

这段代码展示了一个常见的模式。数组sol包含了解决方案的坐标（以列的形式），因此我们将数组转置，使得数据沿数组的行排列，然后将每一行分配给变量x、y和z。

其他几行代码也很简单：我们调用了plot()函数，然后为坐标轴添加标签。以下是我们得到的输出图形：

上面的图像被称为经典的洛伦兹蝴蝶，是一个奇异吸引子的 striking 示例。

使用 Numba 和 NumbaPro 加速计算

本节我们将讨论Numba和NumbaPro，这两个非常激动人心的库，用于加速NumPy代码。Numba和NumbaPro由Continuum Analytics公司创建，后者也生产 Anaconda 发行版。Numba是标准的 Anaconda 发行版的一部分，但NumbaPro是一个商业产品，必须单独购买作为Accelerate包的一部分。不过，NumbaPro可以免费试用一段时间。

这些库的独特之处在于它们允许通过增加几行代码来加速代码执行。作为第一个示例，我们来看看以下几行代码来乘法两个矩阵：

def matrix_multiply(A, B):
 m, n = A.shape
 n, r = B.shape
 C = zeros((m, r), float64)
 for i in range(m):
 for j in range(r):
 acc = 0
 for k in range(n):
 acc += A[i, k] * B[k, j]
 C[i, j] = acc
 return C

上面的代码使用了矩阵乘法的简单定义，看起来就像我们在用 C 语言实现算法时写的代码。这不是 Python 风格的代码，显然也没有经过优化。（在实际情况中，通常会直接使用NumPy的内置矩阵乘法。）特别需要注意的是，矩阵的维度并未进行检查：假设A的列数等于B的行数。

首先让我们尝试使用小矩阵进行计算，如下所示：

A = array([[1,2,0],[1,-3,4],[0,-2,1],[3,7,-4]], dtype=float64)
B = array([[3,4],[-2,0],[2,4]], dtype = float64)
C = matrix_multiply(A, B)
print A
print B
print C

我们首先定义矩阵A和B（注意它们的维度是兼容相乘的）。与本节中的所有示例一样，我们小心地指定了数据类型（这可能有助于优化）。然后，我们简单地调用matrix_multiply，将结果存储在数组C中，并打印这三个矩阵。结果如下所示：

[[  1\.   2\.   0.]
 [  1\.  -3\.   4.]
 [  0\.  -2\.   1.]
 [  3\.   7\.  -4.]]
[[  3    4]
 [ -2    0]
 [  2    4]]
[[ -1\.   4.]
 [  17\. 20.]
 [  6\.   4.]
 [ -13\. -4.]]

你可以通过手动检查几个条目来验证算法的正确性。或者，我们也可以检查结果是否与内置的矩阵乘法一致，如下所示：

C - A.dot(B)
array([[ 0\. ,  0.],
 [ 0\. ,  0.],
 [ 0\. ,  0.],
 [ 0\. ,  0.]])

一切看起来都很好。现在，我们希望定义一些更大的随机矩阵，如下所示：

n = 100
A = rand(n, n)
B = rand(n, n)

在 64 位架构下，前面的代码行会自动生成 64 位浮点数矩阵。接下来，我们将矩阵相乘并计时，如下所示：

%%timeit
C = matrix_multiply(A, B)

上述计算的输出结果如下：

1 loops, best of 3: 472 ms per loop

计时结果当然会因运行代码的机器而异。此示例在运行 Microsoft Windows 7 64 位操作系统的 Intel Core i7 处理器，主频为 3.5 GHz，内存为 16 GB 的计算机上运行。

现在让我们看看如何快速优化这个函数。首先，从Numba模块加载jit函数，如下所示：

from numba import jit

然后，定义带有@jit装饰器的函数，如下所示：

@jit
def matrix_multiply_jit(A, B):
 m, n = A.shape
 n, r = B.shape
 C = zeros((m, r), float64)
 for i in range(m):
 for j in range(r):
 acc = 0
 for k in range(n):
 acc += A[i, k] * B[k, j]
 C[i, j] = acc
 return C

注意，代码中唯一的更改是添加装饰器。（我们还更改了函数的名称以避免混淆，但这并非必要。）装饰器是一个高级的 Python 主题，但我们不需要深入讨论它们的工作原理。有关装饰器的更多信息，请参考 Simeon Franklin 在simeonfranklin.com/blog/2012/jul/1/python-decorators-in-12-steps/上的优秀博客文章。

现在让我们计时我们的代码，如下所示：

%%timeit
C = matrix_multiply_jit(A, B)

以下是结果输出：

1000 loops, best of 3: 1.8 ms per loop

这是一行代码实现的 260 倍改进！在这里应该保持适度的观点，因为这种加速不能期望适用于通用代码。请记住，我们故意编写代码，不使用NumPy中已经优化过的函数。为了比较和充分披露，让我们将其与内置的dot()方法进行比较：

%%timeit
C = A.dot(B)

结果输出如下：

10000 loops, best of 3: 28.6 µs per loop

因此，即使加速了，我们的函数也无法与内置的NumPy竞争。我们再次强调，本节的目标是介绍加速技术的概述，而不是深入探讨复杂的优化方法。

了解@jit装饰器的工作原理是值得的。当调用由@jit装饰的函数时，库会尝试推断参数和返回值的数据类型，并动态生成函数的编译版本，然后调用它。结果是一个类似于用 C 编写的代码的函数调用。

可以指定数据类型，而不是让参数和返回值的类型被推断，这可能会提高性能。以下表格列出了Numba支持的数据类型及其缩写：

数据类型	缩写
boolean	b1
bool_	b1
byte	u1
uint8	u1
uint16	u2
uint32	u4
uint64	u8
char	i1
int8	i1
int16	i2
int32	i4
int64	i8
float_	f4
float32	f4
double	f8
float64	f8
complex64	c8
complex128	c16

所有这些名称都在Numba模块中定义。例如，要定义一个将两个浮点值相加的函数，我们使用以下代码：

from numba import jit, f8
@jit (f8(f8,f8))
def my_sum(a, b):
 return a + b

注意装饰器语法，如下所示：

@jit (f8(f8,f8))

这指定了一个接受两个 float64 参数并返回一个 float64 值的函数。然后，函数按如下方式调用：

my_sum(3.5, 6.9)

这将产生预期的结果。然而，如果我们尝试类似以下代码，我们会遇到错误：

a = array([1.4, 2.0])
b = array([2.3, 5,2])
my_sum(a,b)

然而，使用 @jit 装饰器可以使用数组。要定义一个将两个一维数组相加的函数，可以使用以下代码行：

@jit (f8:)
def vector_sum(a, b):
 return a + b

注意如何指定向量。二维数组表示为 f8[:,:]，三维数组表示为 f8[:,:,:]，以此类推。

NumbaPro 是 Numba 的商业版，并增加了多个增强功能。我们将重点讨论使用图形处理单元（GPU）作为示例，进行并行处理，这是一个激动人心的新技术，现已轻松应用于笔记本中。

要运行以下示例，读者必须拥有 NumbaPro、CUDA 兼容的 GPU（以下简称“设备”）以及最新的 CUDA 兼容驱动程序。

CUDA 兼容设备的列表可以在developer.nvidia.com/cuda-gpus找到。验证您的设备兼容性后，请从developer.nvidia.com/cuda-downloads下载并安装适合平台的最新版本 CUDA SDK。CUDA 工具包附带了一些示例，您可以使用这些示例来测试安装情况。

NumbaPro 下载可以通过store.continuum.io/cshop/accelerate/获得。下载并安装 Accelerate 库。

要测试设置，启动一个 IPython 笔记本，并在单元格中运行以下内容：

import numbapro
numbapro.check_cuda()

如果一切正常，系统将打印出由 Anaconda 安装的 CUDA 库的列表以及系统中可用的 CUDA 兼容设备的列表。您还会看到显示末尾的 PASSED 消息。

尽管 CUDA 编程是通向大规模并行计算的相对简单路径，但在我们开始第一个 CUDA 程序之前，仍然需要掌握一些概念。我们将在此简要概述架构的基础内容，只讨论足够运行后续示例的内容。有关完整规格，请参见CUDA 编程指南，该指南可通过docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#programming-model获得。

GPU 最初是为比计算机的 CPU 更快速地处理渲染操作而设计的。这种处理加速在很大程度上是通过大规模并行化渲染管道所需的图形操作来实现的。

一个CUDA兼容的 GPU 由一组流多处理器（SMs）组成。每个 SM 本身无法与当前的 CPU 在速度上竞争。然而，多个 SM 能够合作解决问题，这一点弥补了这一不足。SM 还可以访问 GPU 中的内存，这部分内存称为设备内存。

在CUDA中，运行在 CPU 中的代码和运行在设备（GPU）中的代码有严格的分离。一个特别的限制是，CPU 代码只能访问常规计算机内存，而设备代码只能访问设备内存。运行在设备中的代码被指定在一个名为内核的函数中。内核被编译成 SM 可以理解的低级语言，并异步地运行在每个 SM 中（意味着每个 SM 按照自己的节奏运行，除非遇到特殊的同步指令）。因此，CUDA中的简单计算通常需要以下三个步骤：

1. 将输入数据从计算机内存传输到设备内存。

2. 在设备中启动内核。

3. 将设备内存中的输出数据传输到计算机内存，以便 CPU 能够再次访问它。

正如你将看到的，Python CUDA使内存传输变得透明。（如果需要，仍然可以通过编程方式控制它。）

内核同时在一组 SM 中启动，每个线程独立进行计算。每个 SM 可以并行运行多个线程，并且可以访问设备的所有内存。（架构更为复杂，还有其他类型的内存可用，这里不作讨论。）在最简单的情况下，每个线程只访问几个内存区域，每个区域包含一个 64 位浮动值，而且线程访问的内存不会被其他线程访问。因此，无需同步。在更复杂的问题中，同步可能成为一个主要问题。

正在运行的线程集合具有二级数组结构：

• 线程按块组织。每个块是一个最多具有3个维度的线程数组。块的维度存储在一个名为blockDim的变量中。块中的线程由变量threadIdx标识。这是一个具有三个整数字段的结构：threadIdx.x、threadIdx.y和threadIdx.z。这些字段唯一标识块中的每个线程。

• 块按网格组织。网格是一个最多具有3个维度的块数组。网格的维度存储在一个名为gridDim的变量中。网格中的块由变量gridIdx标识。这是一个具有三个整数字段的结构：gridIdx.x、gridIdx.y和gridIdx.z。这些字段唯一标识网格中的每个块。

下面的图给出了这种组织结构的一个示例：

在前面的示例中，gridDim是(2, 3, 1)，因为有两行三列的块（以及一个单独的空间维度）。网格中的所有块都是一维的，所以blockDim是(4, 1, 1)。例如，底部行中第一个块中的第三个线程由以下代码行标识：

blockIdx.x=0, blockIdx.y=1, blockIdx.z=1
threadIdx.x=2, threadIdx.y=1, threadIdx.z=1

在运行时，每个线程都可以访问这些标识信息。

CUDA架构的一个关键点如下：

• 同一块中的所有线程始终在一个 SM 中并行执行，直到块中的所有线程都执行完毕

• 不同的块可以根据 SM 的可用性并行或串行执行计算

我们现在准备使用 Python CUDA定义内核。我们将编写一个函数，在 GPU 中计算两个向量的和。在一个单元格中运行以下代码：

from numbapro import cuda
@cuda.jit('void(float64[:], float64[:], float64[:])')
def sum(a, b, result):
 i = cuda.threadIdx.x 
 result[i] = a[i] + b[i]

我们假设只有一个线程块，每个线程负责在单个位置添加数组的元素。线程负责的数组位置由threadIdx.x的值确定。请注意，内核没有返回值。我们需要指定一个数组result，用于存储计算的返回值。

现在让我们看看这个函数是如何被调用的。请注意，网格和块的几何形状并没有在内核中定义。（如果需要，内核可以获取几何信息；稍后会详细讲解。）这在启动内核时完成：

a = array([1,2,3,4,5], dtype=float64)
b = array([5,4,3,2,1], dtype=float64)
result = array([0,0,0,0,0], dtype=float64)
sum1,5
print result

上面的代码行给出了以下输出：

[ 6\.  6\.  6\.  6\.  6.]

这段代码的关键点在于以下这一行：

sum1,5

上面的代码行在一个包含1个块、每个块包含5个线程的网格中启动内核。网格和块都是一维的。现在让我们添加更大的向量：

n = 64
a = arange(0,n,dtype=float64)
b = arange(n,0,-1, dtype=float64)
result = zeros(n,dtype=float64)
sum1,n
print result[:5]

上面的代码行本质上与之前的代码相同，只是稍微通用了一些，因为数组的大小可以变化。我们想要做的是增加n的大小。如果你尝试使用n=10000这样的值，将会出现CUDA_ERROR_INVALID_VALUE类型的错误。问题在于单个 SM 可以运行的线程数是有限的，也就是说，单个块中可以执行的线程数有限。为了处理更大的向量，我们需要修改代码，使其能够处理多个块。为此，请以以下方式修改sum()函数的定义：

from numbapro import cuda
@cuda.jit('void(float64[:], float64[:], float64[:], int32)')
def sum(a, b, result, n):
 tx = cuda.threadIdx.x
 bx = cuda.blockIdx.x
 bsz = cuda.blockDim.x
 i = tx + bx * bsz
 if i < n:
 result[i] = a[i] + b[i]

首先需要注意的是，我们包含了一个类型为int32的参数，用于保存正在加法的数组的大小。现在的重点是，位于不同块中的线程必须访问不同的内存区域，因此计算与线程相关联的索引i变得更加复杂。本质上，我们必须知道当前块之前有多少个块，然后将其乘以块的维度，并加上当前线程的索引。然后，在访问相关的内存位置之前，我们会检查索引是否有效。这可以防止线程访问不属于输入/输出数组的区域，并且是更复杂代码中一个至关重要的检查。要测试代码，请运行以下内容：

n = 100000
a = arange(0,n,dtype=float64)
b = arange(n,0,-1, dtype=float64)
result = zeros(n,dtype=float64)
sum1000,64
print result[:5]

上述代码应该可以顺利运行。请注意，我们指定了一个包含1000个块且每个块有64个线程的网格。网格中的块数是无限制的，设备负责以最佳方式分配 SM。请注意，块数必须足够大，以覆盖输入/输出数组。在我们的例子中，这意味着blockDim.x * gridDim.x >= n。

现在我们可以开始使用大向量进行计算。试试下面的代码：

n = 100000
a = rand(n)
b = rand(n)
result = zeros(n, dtype=float64)
bd = 10000
gd = 64
if bd * gd < n:
 print 'Block/grid dimensions too small'
else:
 sumbd,gd
print result[:10]

读者应尝试不同的n、bd和gd值。请记住，gd的最大值取决于你电脑中的设备。一个有趣的实验是检查计算在更大n值下的扩展性。

概述

本章我们介绍了在SciPy中使用高级数学算法，包括解方程、寻找最优值、积分和微分方程。本章最后讨论了如何利用 GPU 并行化来加速计算。

附录 A. IPython 笔记本参考卡

启动笔记本

要启动笔记本，请打开终端窗口并运行以下命令：

ipython notebook

在运行上述命令时，确保你处于包含笔记本的目录中。

键盘快捷键

一些重要的键盘快捷键如下：

• 要进入编辑模式，按下Enter键或点击单元格

• 要进入命令模式，按Esc

编辑模式中的快捷键

编辑模式中使用的一些重要快捷键如下：

• 要运行一个单元格，使用以下快捷键：

  • 要运行一个单元格并跳转到下一个单元格，按Shift + Enter

  • 要运行一个单元格但保持在同一单元格中，按Ctrl + Enter

  • 要运行一个单元格并在其下方插入一个新单元格，按Alt + Enter

  • 要在当前单元格中创建一个新行，按Enter

• 要缩进内容，按Tab

• 要启动代码补全，开始在单元格中输入然后按Tab

• 要全选，按Ctrl + A

• 要撤销一个操作，按Ctrl + Z

• 要重做一个操作，按Ctrl + Y或Ctrl + Shift + Z

• 要跳转到单元格的开头，按Ctrl + Home

• 要跳转到单元格的末尾，按Ctrl + End

• 要拆分一个单元格，按Ctrl + Shift + -

命令模式中的快捷键

• 要列出键盘快捷键，按H

• 要将单元格模式更改为以下之一，快捷键如下：

  • 代码：按Y

  • Markdown：按M

  • 标题：按 1 到 6 之间的数字，选择标题的级别

  • 原始 NBConvert：按R

• 要选择当前单元格上方的单元格，按上箭头或K

• 要选择当前单元格下方的单元格，按下箭头或J

• 要将一个单元格向上移动一个位置，按Ctrl + K

• 要将一个单元格向下移动一个位置，按Ctrl + J

• 要在当前单元格上方插入一个新单元格，按A

• 要在当前单元格下方插入一个新单元格，按B

• 要剪切一个单元格，按X

• 要复制一个单元格，按C

• 要将一个单元格粘贴到当前单元格下方，按V

• 要将一个单元格粘贴到当前单元格上方，按Shift + V

• 要删除一个单元格，按D

• 要撤销删除操作，按Z

• 要将当前单元格与下方的单元格合并，按Shift + M

• 要切换行号显示，按L

导入模块

加载一些重要模块的步骤如下：

• 要加载NumPy和 matplotlib 以进行交互式工作，并支持内联图形，请运行以下命令：

  pylab inline
  
  

• 要加载NumPy和 matplotlib 而不将名称导入当前命名空间，并且支持内联图形，请运行以下命令行：

  pylab -–no-import-all inline
  
  

• 要加载SciPy模块，使用以下任何标准的 Python 导入命令：

  import scipy.<module>
  import scipy.<module> as <local-name>
  from scipy.<module> import <function>
  
  

如果使用了–no-import-all选项，函数和对象必须加上适当的模块名称前缀，如下所示：

• 对于NumPy函数和对象，使用numpy或np。

• 对于交互式图形，使用pyplot或plt。

系统中安装的库模块以及用户创建的.py扩展名的模块，可以通过标准的 Python 机制进行导入。

获取帮助

获取帮助的方式有很多种：

• 要启动交互式帮助，可以运行以下命令：

  help()
  
  

• 要获取函数或对象的帮助，可以运行以下命令：

  help(object)
  help(function)
  <object>?
  <function>?
  
  

• 对于自动补全，开始输入函数/对象/方法的名称，然后按Tab。

• 要获取工具提示，开始输入函数/对象/方法的名称，然后按Shift + Tab。

附录 B. Python 简要回顾

介绍

本附录将简要介绍 Python 语法。这并不是一个 Python 编程课程，而是供不熟悉该语言的读者作为快速入门使用。以下主题将在本附录中涵盖：

• 基本类型、表达式、变量及其赋值

• 序列类型

• 字典

• 控制结构

• 函数、对象和方法

基本类型、表达式、变量及其赋值

任何可以在 Python 代码中引用的数据都被视为对象。对象用于表示从原子数据（如数字）到非常复杂的数据结构（如多维数组、数据库连接和各种格式的文档）。

在对象层次结构的根部是数值数据类型。这些包括以下内容：

• 整数：Python 中有三种类型的整数。

  • 普通整数：它们在本地架构中表示，大多数系统中通常是 32 位或 64 位的有符号值。

  • 长整型：它们是具有无限范围的整数，取决于可用内存。大多数情况下，程序员不需要关心普通整数和长整型之间的区别。Python 会透明地处理这两种类型之间的转换。

  • 布尔值：它们表示False和True两个值。在大多数情况下，它们分别等价于0和1。

• 浮点数：它们表示本地的双精度浮点数。

• 复数：它们表示复数，表示为一对双精度浮点数。

下表展示了每种数据类型的字面量（即常量）示例：

数据类型	字面量
整数	0, 2, 4, …, 43882838388``5L, 5l（长整型）0xFE4（十六进制）03241（八进制）
实数（float）	5.34, 1.2, 3., 0``1.4e-32（科学记数法）
复数	1.0+3.4j, 1+2j, 1j, 0j, complex(4.3, 2.5)

虚数单位由j表示，但只有在它跟随数字字面量时（否则它表示名为j的变量）。因此，要表示虚数单位，必须使用1j，复数零为0j。复数的实部和虚部总是以双精度浮点值存储。

注意

请注意，NumPy大大扩展了数值类型的集合，以便进行高效的数值计算。

赋值语句用于将值存储到变量中，如下所示：

a = 3
b = 2.28
c = 12
d = 1+2j

Python 支持多重同时赋值，因此前四行代码可以等效地写成一行，如下所示：

a, b, c, d = 3, 2.28, 12, 1+2j

在多重赋值中，右侧的所有表达式都会在赋值之前被求值。例如，交换两个变量的值的常见惯用法如下所示：

v, w = w, v

作为练习，读者可以尝试预测以下语句的结果，前提是已知前面的变量赋值：

a, b, c = a + b, c + d, a * b * c * d
print a, b, c, d

以下示例展示了如何计算二次方程的两个解：

a, b, c = 2., -1., -4.
x1, x2 = .5 * (-b - (b ** 2 - 4 * a * c) ** 0.5), .5 * (-b + (b ** 2 - 4 * a * c) ** 0.5)
print x1, x2

请注意，我们通过使用小数点强制将变量 a、b 和 c 转换为浮点值。在进行数值计算时，这是一种良好的实践。下表包含了 Python 运算符的部分列表：

运算符	Python 运算符
算术	+（加法）-（减法，取负）*（乘法）/（除法，见表格下方的注释）//（整数除法）%（余数）
比较	==（等于）>（大于）<（小于）>=（大于或等于）<=（小于或等于）!=（不等于）
布尔	and or not
位运算布尔	&（与）&#124;（或）^（异或）~（非）
位移	<<（左移）>>（右移）

注释

使用除法运算符（/）时需要小心。如果操作数是整数，则此操作的结果是整数商。例如，34/12 结果为 2。要获得浮点结果，我们必须输入浮点操作数，例如 34./12.，或者添加以下语句：

from __future__ import division

// 运算符始终表示整数除法。

算术运算符遵循运算顺序规则，括号可以改变该顺序。比较运算符的优先级低于算术运算符，而 or、and 和 not 运算符的优先级更低。因此，像以下的表达式会得到预期的结果：

2 + 3 < 5 ** 2 and 4 * 3 != 13

换句话说，前面的命令行被解析如下：

(((2 + 3) < (5 ** 2)) and ((4 * 3) != 13))

逻辑运算符 and 和 or 采用短路求值，因此，例如，命令中的第二个比较不会被评估：

2 < 3 or 4 > 5

位运算符和移位运算符的优先级规则可能不像直观那样清晰，因此建议始终使用括号来指定运算顺序，这也有助于代码的清晰度。

Python 还支持增强赋值。例如，以下命令行首先将值 5 赋给 a，然后将 a 的值加一：

a = 5
a += 1

注释

Python 不支持递增/递减运算符，如 C 语言中的 a++ 和 ++a。

所有 Python 运算符都有对应的增强赋值语句。任何运算符 $ 的一般语义如下所示：

v $= <expression>

前面的语句等价于以下：

v = v $ (<expression>)

注释

请注意，$ 不是一个有效的 Python 运算符，它仅作为一个占位符表示通用运算符。

序列类型

Python 序列类型用于表示有序的对象集合。它们分为可变序列类型和不可变序列类型。在这里，我们只讨论lists（可变）和tuples与strings（均为不可变）。其他序列类型将在本节末尾提到。

列表

以下示例演示如何在 Python 中构造一个列表并将其赋值给变量：

numbers = [0, 1.2, 234259399992, 4+3j]

列表中的个别条目可以通过索引表示法访问，如下所示：

numbers[2]

请注意，索引总是从0开始。允许使用负索引，它们表示从列表末尾开始的位置。例如，numbers[-1]是最后一个元素，numbers[-2]是倒数第二个元素，以此类推。

由于列表是可变的序列类型，我们可以就地修改其中的条目：

numbers[0] = -3
numbers[2] += numbers[0]
print numbers

另一种重要的引用 Python 序列类型元素的方法是切片，它允许从列表中提取子列表。由于这个话题对NumPy数组非常重要，我们将讨论延迟到附录 C，NumPy 数组。

Python 列表具有一组非常实用的功能，其中一些在以下代码示例中得到了体现：

• 要查找列表的长度，请使用以下命令：

  len(numbers)
  
  

• 要就地反转列表，请使用以下命令：

  numbers.reverse()
  print numbers
  
  

• 要添加一个新元素，请使用以下命令：

  numbers.append(35)
  print numbers
  
  

• 要就地对列表进行排序，请使用以下命令：

  values = [1.2, 0.5, -3.4, 12.6, 3.5]
  values.sort()
  print values
  values.sort(reverse=True)
  print values
  
  

• 要在指定位置插入一个值，请使用以下命令：

  values.insert(3, 6.8)
  print values
  
  

• 要扩展列表，请使用以下命令：

  values.extend([7,8,9])
  print values
  
  

Python 有一些便捷的方式来构造常用的列表。range()函数返回一个等间隔的整数列表。最简单的形式返回一个从0开始的连续整数列表：

range(10)

上述命令返回以下列表：

[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

请注意，最后一个元素比函数调用中给出的参数少 1。通常规则是，range(n)返回一个包含n个元素的列表，从零开始，因此最后一个元素是n-1。要从非零值开始，请使用以下带有两个参数的版本：

range(3, 17)

第三个参数指定一个增量。以下命令行会生成一个包含所有小于 100 的正 6 的倍数的列表：

range(6,100,6)

也可以使用负增量：

range(20, 2, -3)

列表支持连接，表示为+运算符：

l1 = range(1, 10)
l2 = range(10, 0, -1)
l3 = l1 + l2
print l3

注意

请注意，对于NumPy数组，+运算符被重新定义为表示向量/矩阵加法。

乘法运算符（*）可用于通过重复给定列表的元素来构造一个新列表，如下所示：

l4 = 3*[4,-1,5]
print l4

在 Python 中构造列表的最灵活方法是使用列表推导式。完整的讨论超出了本附录的范围，但以下示例说明了其中的一些可能性：

• 要显示从0到10（包括10）的整数的平方列表，请使用以下命令行：

  [n ** 2 for n in range(11)]
  
  

• 要显示一个整数的所有除数列表，请使用以下命令行：

  k = 60
  [d for d in range(1, k+1) if k % d == 0]
  
  

• 要显示小于 100 的所有素数列表，请使用以下命令行（效率较低）：

  [k for k in range(2,101) if len([d for d in range(1, k+1) if k % d == 0])==2]
  
  

• 要显示具有整数坐标的点的元组列表及其与原点的距离，请使用以下命令行：

  [(i,j,(i*i+j*j)**0.5) for i in range(5) for j in range(6)]
  
  

元组

元组与列表类似，但它们是不可变的——一旦创建，其元素就不能修改。

会发生变化。以下命令行将导致错误信息：

t1 = (2,3,5,7)
t1[2] = -4

元组在 Python 中有一些特殊用途。它们可以作为字典中的索引（因为它们是不可变的）。它们还构成了 Python 用来从函数返回多个值的机制。例如，内建函数 divmod() 返回一个元组，其中包含整数商和余数：

divmod(213, 43)

元组支持与列表相同的序列接口，除了会修改元组的那些方法。例如，没有名为 sort() 的方法可以就地对元组进行排序。

字符串

Python 字符串表示一个不可变的字符序列。有两种字符串类型：str，表示 ASCII 字符串，和 unicode，表示 Unicode 字符串。

字符串文字是由单引号或双引号括起来的字符序列，如下所示：

s1 = 'I am a string'
s2 = "I am a string"
print s1
print s2

单引号和双引号之间没有语义上的区别，唯一的不同是单引号字符串可以包含双引号，而双引号字符串可以包含单引号。例如，以下命令行是正确的：

s3 = "I'm a string"
print s3

字符串有两个主要用途：作为字典的索引和打印信息。当打印信息时，字符串有一个 format() 方法，可以轻松地显示信息。我们经常使用这个功能为图形添加注释。这里有一个例子：

n = 3
message = 'The square root of {:d} is approximately {:8.5f}.'.format(n, n ** 0.5)
print message

在前面的例子中，有两个格式说明符：

• {:d}：这指定了整数值的十进制格式

• {:8.5f}：这指定了一个宽度为 8，小数点后有 5 位的浮点数格式

格式说明符会按顺序与参数匹配，在这个例子中是 n 和 n ** 0.5。

字符串有丰富的接口。如果你需要编写与字符串相关的代码，很可能有一个内建函数可以实现，且只需要很少的修改。所有可用字符串方法的列表以及格式化功能，可以参考 docs.python.org/2/library/stdtypes.html#string-methods。

字典

Python 字典是一种包含键值对的数据结构。键必须是不可变类型，通常是字符串或元组。以下是一个构建字典的示例：

grades = {'Pete':87, 'Annie':92, 'Jodi':78}

要访问某个项，我们提供键作为索引，如下所示：

print grades['Annie']

字典是可变的，所以我们可以使用它们来更改条目的值。如果 Jodi 做了额外的工作来提高她的成绩，我们可以这样更改：

grades['Jodi'] += 10
print grades['Jodi']

要向字典添加条目，只需为新键赋值：

grades['Ivan']=94

然而，尝试访问一个不存在的键会导致错误。

一个需要注意的重要点是，字典是无序的。以下代码是迭代字典的标准惯用法：

for key, item in grades.iteritems():
 print "{:s}'s grade in the test is {:d}".format(key, item)

这里的关键点是，输出与字典中条目的添加顺序完全无关。

如需了解更多字典接口的细节，可以参考docs.python.org/2/library/stdtypes.html#mapping-types-dict。

控制结构

控制结构允许改变代码执行的流程。我们关注的有两种结构：分支和循环。

分支结构根据测试结果执行不同的代码。以下示例展示了一个改进版的代码，用于求解二次方程。使用了if-then-else结构来处理实数解和虚数解的情况，如下所示：

a, b, c = 2., -4., 5.
discr = b ** 2 - 4 * a * c
if discr >= 0:
 sqroot = discr ** 0.5
 x1 = 0.5 * (-b + sqroot)
 x2 = 0.5 * (-b - sqroot)
else:
 sqroot = (-discr) ** 0.5
 x1 = 0.5 * (-b + sqroot * 1j)
 x2 = 0.5 * (-b - sqroot * 1j)
print x1, x2

上面的代码首先计算二次方程的判别式。然后，使用if-then-else语句来判断根是实数还是虚数，这取决于判别式的符号。注意代码的缩进，Python 通过缩进来定义语句块的边界。if-then-else结构的一般形式如下：

if <condition>:
 <statement block T>
else:
 <statement block F>

首先，条件<condition>被评估。如果为True，则执行语句<statement block T>；否则，执行语句<statement block F>。else:子句可以省略。

Python 中最常见的循环结构是for语句。下面是一个示例：

numbers = [2, 4, 3, 6, 2, 1, 5, 10]
for n in numbers:
 r = n % 2
 if r == 0:
 print 'The integer {:d} is even'.format(n)
 else:
 print 'The integer {:d} is odd'.format(n)

我们首先定义一个整数列表。for语句使得变量n依次取列表中的每个值，并对每个值执行缩进块的代码。注意，for循环内部有一个if-then-else结构。此外，print语句是双重缩进的。

for循环常用于执行简单的搜索。一个常见的场景是，当满足某个条件时需要跳出循环。以下代码查找整数范围内的第一个完全平方数：

for n in range(30, 90):
 if int(n ** 0.5) ** 2 == n:
 print n
 break

对于给定范围内的每个n值，我们先计算n的平方根，取其整数部分，再计算平方。如果结果等于n，则进入if块，打印n，然后跳出循环。

如果在范围内没有完全平方数怎么办？将前面的函数range(30, 60)改为range(125, 140)。运行命令行时不会输出任何内容，因为在125到140之间没有完全平方数。现在，将命令行改为如下：

for n in range(125, 140):
 if int(n ** 0.5) ** 2 == n:
 print n
 break
else:
 print 'There are no perfect squares in the range'

else子句只有在执行没有跳出循环时才会被执行，此时会打印出相应的消息。

另一个常见的情况是某些值在迭代过程中需要跳过。以下示例中，我们打印一系列-1到1之间随机数的平方根，但只有当这些数为正时才打印：

import random
numbers = [-1 + 2 * rand() for _ in range(20)]
for n in numbers:
 if n < 0:
 continue
 print 'The square root of {:8.6} is {:8.6}'.format(n, n ** 0.5)

当 Python 遇到continue语句时，它会跳过当前执行块的剩余部分，并继续执行控制变量的下一个值。

另一个经常使用的控制结构是while循环。该结构在条件为真时执行一组命令。例如，假设我们想计算一组随机生成的值的累加和，但仅当和超过某个值时才停止。可以使用以下代码来实现：

import random
bound = 10.
acc = 0.
n = 0
while acc < bound:
 v = random.random()
 acc += v
 print 'v={:5.4}, acc={:6.4}'.format(v, acc)

另一种比预期更常见的情况需要一种叫做永远循环的模式。这种情况发生在循环开始时需要检查的条件不可用时。例如，下面的代码实现了著名的3n+1游戏：

n = 7
while True:
 if n % 2 == 0:
 n /= 2
 else:
 n = 3 * n + 1
 print n
 if n == 1:
 break

游戏从一个任意整数开始，这里是7。然后，在每次迭代中，我们测试n是否为偶数。如果是，我们将其除以2；否则，将其乘以3并加上1。然后，我们检查是否达到了1。如果是，我们退出循环。由于我们不知道是否需要在循环结束前退出，因此我们使用一个永远循环，如下所示：

while True:
 <statements>
 if <condition>:
 break
 <possibly more statements>

一些程序员避免使用这种结构，因为如果不小心，很容易导致无限循环。然而，它在某些情况下确实非常有用。顺便说一下，3n+1问题中的循环是否对所有初始值都会停止还是一个开放性问题！读者可以尝试使用初始值n=27来感受一下。

函数、对象和方法

现在我们来介绍使 Python 如此灵活和强大的构造，它的面向对象特性。我们在前面的章节中已经看到了一些面向对象的代码示例（面向对象的范式是 Python 的核心，很难写出不使用它的代码），但现在我们将对这些特性进行更具体的讨论。

函数

我们已经看到许多使用函数的例子。例如，len()函数用于计算列表的长度：

lst = range(1000)
print len(lst)

调用函数的最基本语法如下：

function_name(arg1, arg2, …, argn)

在这种情况下，arg1、arg2、…、argn被称为位置参数，因为它们是根据出现的位置进行匹配的。例如，我们考虑内置函数pow()。这个函数最多接受三个参数：

pow(b, n, m)

在这种形式下，前面的函数使用优化算法来计算b的n次方模m。（如果你想知道，这是公共密钥密码学中的一个重要操作。）参数b、n和m是根据它们的位置来关联的。例如，要计算12的十次方模15，我们使用以下命令：

pow(12, 10, 15)

Python 还支持任意大小的参数序列。例如，max()函数计算任意序列中的最大值：

max(2,6,8,-3,3,4)

前面的命令返回值8。

传递参数给函数的第三种方式是使用关键字参数。这非常有用，因为通常很难记住位置参数的准确顺序。（例如，我不太愿意编写超过三个或四个位置参数的函数。）

例如，内置的int()函数可以用于将字符串转换为整数。可选的关键字参数base让我们可以指定转换的进制。例如，以下命令行将给定的2进制整数赋值给n：

n = int('10111010100001', base=2)
print n

关键字参数总是有默认值。在我们的示例中，如果没有指定基数，默认假设为10。

我们经常需要编写自己的函数。这是通过关键字def来实现的。作为示例，假设我们要编写代码来实现著名的bisection方法以数值求解方程。一个可能的解决方案如下：

def bisection(f, a, b, tol=1e-5, itermax=1000):
 fa = f(a)
 fb = f(b)
 if fa * fb > 0:
 raise ValueError('f(a) and f(b) must have opposite signs')
 niter = 0
 while abs(a-b) > tol and niter < itermax:
 m = 0.5 * (a + b)
 fm = f(m)
 if fm * fa < 0:
 b, fb = m, fm
 else:
 a, fa = m, fm
 return min(a, b), max(a, b)

上述函数需要三个重要且必要的参数：

• f函数接受一个浮点值作为输入，并返回一个浮点值作为输出。

• 浮点值a和b，它们指定一个包含函数零点的区间。

另外两个参数是可选的。参数tol指定结果所需的容差，itermax指定最大迭代次数。要使用bisection()函数，必须先定义函数f。我们将借此机会展示另一种定义 Python 函数的方式，如下：

from math import cos, pi
f = lambda x: cos(x) - x

我们现在准备通过以下命令来调用函数：

bisection(f, 0, pi/2)

上述函数返回以下输出：

(0.7390851262506977, 0.7390911183631504)

请注意，我们设计了该函数返回一个包含零的区间。该区间的长度小于tol，除非达到了最大迭代次数。如果我们希望更小的容差，可以使用以下函数：

bisection(f, 0, pi/2, tol=1E-10)

现在，假设我们关注计算所花费的时间。我们可以通过以下方式限制最大次数：

bisection(f, 0, pi/2, itermax=10, tol=1E-20)

请注意，关键字参数的顺序是无关紧要的，并且在前面的示例中没有达到所需的容差。

对象和方法

对象是 Python 中最通用的数据抽象。实际上，从程序员的角度来看，Python 中的一切都是对象。

对象不过是结构化数据的集合，并且有一个操作这些数据的接口。对象是通过class构造定义的，但我们这里的目标并不是展示如何定义类。尽管设计新类是一个高级话题，但使用现有类却相当简单。

作为示例，让我们探索内置类型str。首先，我们定义一个可以操作的str对象，如下所示：

message = 'Mathematics is the queen of science'

首先，让我们将消息转换为大写，如下所示：

message.upper()

我们说前述语句调用了message对象的upper()方法。方法仅仅是与对象关联的函数。以下是str对象的其他一些方法：

• 要查找子字符串的第一次出现（如果未找到该字符串，返回-1），请使用以下命令行：

  message.find('queen')
  
  

• 要将字符串拆分为单词，请使用以下命令行：

  words = message.split()
  print words
  
  

• 要计算s子字符串出现的次数，请使用以下命令行：

  message.count('e')
  
  

• 要用其他内容替换子字符串，请使用以下命令行：

  message.replace('Mathematics', 'Mme. Curie')
  
  

注意

请注意，前述方法不会改变原始字符串对象，而是返回新的修改过的字符串。字符串是不可变的。对于可变对象，方法可以自由地改变对象中的数据。

总结

在本附录中，我们概述了 Python 语法和特性，涵盖了基本类型、表达式、变量和赋值、基本数据结构、函数、对象和方法。

附录 C. NumPy 数组

介绍

数组是NumPy引入的基本数据结构，它们是我们在本书中讨论的所有科学计算和数据分析库的基础。本附录将简要概述以下数组特性：

• 数组创建和成员访问

• 索引和切片

数组创建和成员访问

NumPy数组是ndarray类的对象，该类表示一个固定大小的多维同质数据集合。

在这里，我们假设已经通过以下命令行导入了NumPy库：

import numpy as np

一旦完成这些步骤，我们可以从一个列表的列表创建ndarray（以后简称为数组对象或简单称为数组），如以下命令行所示：

a = np.array([[-2,3,-4,0],[2,-7,0,0],[3,-4,2,1]],dtype=np.float64)
print a

与 Python 列表和元组不同，数组对象的所有元素必须是相同类型的。这些类型由NumPy对象表示，并称为数组的dtype（数据类型）。在前面的示例中，我们明确指定dtype为float64，它表示一个 64 位浮动小数值。

数组具有几个属性，用于提供关于数据布局的信息。最常用的属性如下：

• 数组的形状可以通过以下命令计算：

  a.shape
  
  

  前面的命令返回元组(3, 4)，因为这是一个具有三行四列的二维数组。有些人可能会感到惊讶，shape属性并非只读，我们可以利用它来重塑数组：

  a.shape = (6,2)
  print a
  
  

  运行前面的示例后，运行a.shape(3,4)可以返回原始维度。

• 数组的维度数量可以通过以下命令获取：

  a.ndim
  
  

  当然，这将返回2。在NumPy中，一个重要的概念是数组的轴。二维数组有两个轴，编号为 0 和 1。如果我们把数组看作是一个数学矩阵，那么轴 0 是垂直的，指向下方，轴 1 是水平的，指向右侧。某些数组方法有一个可选的axis关键字参数，允许用户指定在哪个轴上执行操作。

• 要获取数组中元素的个数，可以使用以下命令：

  a.size
  
  

  在前面的示例中，返回的输出是12，正如预期的那样。

• 数组的一个最终属性是计算数组的转置。这可以通过以下命令完成：

  b = a.T
  print b
  
  

  这创建的一个重要内容是数组a的视图。NumPy包被设计成能高效处理非常大的数组，并且在大多数情况下，除非绝对必要，或明确指示，否则避免复制数据。

• 运行以下代码行：

  print a
  b[1,2] = 11
  print a
  
  

  请注意，数组a的条目2, 1已被更改，这表明变量a和b都指向内存中的同一位置。

• 可以通过empty()函数如下创建一个包含未初始化数据的数组：

  c = np.empty(shape=(3,2), dtype=np.float64)
  print c
  
  

• 使用未初始化的数据是不推荐的，因此最好使用 zeros() 或 ones() 函数，如下所示：

  • 要使用 zeros() 函数，执行以下命令：

    d = np.zeros(shape=(3,2), dtype=np.float64)
    print d
    
    

  • 要使用 ones() 函数，执行以下命令：

    e = np.ones(shape=(3,2), dtype=np.float64)
    print e
    
    

    a_like = np.zeros_like(a)
    print a_like
    
    

  也有一些函数可以创建具有与现有数组相同形状和数据类型的新数组：

• ones_like() 和 empty_like() 函数生成与给定数组相同形状的全 1 数组和未初始化的数据数组。

• NumPy 还有一个 eye() 函数，可以返回给定维度和 dtype 的单位矩阵：

  f = np.eye(5, dtype=np.float64)
  print f
  
  

  行和列的数量不必相同。在这种情况下，生成的矩阵将仅是左单位矩阵或右单位矩阵，具体取决于情况：

  g = np.eye(5, 3, dtype=np.float64)
  print g
  
  

• 数组也可以从现有数据创建。copy() 函数可以如下面这样克隆数组：

  aa = np.copy(a)
  print a
  print aa
  
  

• frombuffer() 函数从暴露（单维）缓冲区接口的对象创建数组。以下是一个示例：

  ar = np.arange(0.0, 1.0, 0.1, dtype=np.float64)
  v = np.frombuffer(ar)
  v.shape = (2, 5)
  print v
  
  

  arange() 函数是 NumPy 对 Python range 的扩展。它的语法类似，但允许浮动范围值。

• loadtxt() 函数从文本文件中读取数组。假设文本文件 matrix.txt 包含以下数据：

   1.3  4.6  7.8
  -3.6  0.4  3.54
   2.4  1.7  4.5
  
  

  然后，我们可以使用以下命令读取数据：

  h = np.loadtxt('matrix.txt', dtype=np.float64)
  print h
  
  

• 数组也可以以 .npy 格式保存和加载：

  np.save('matrix.npy',h)
  hh = np.load('matrix.npy')
  print hh
  
  

索引和切片

为了说明索引，我们先使用以下命令创建一个包含随机数据的数组：

import numpy.random
a = np.random.rand(6,5)
print a

这会创建一个形状为 (6,5) 的数组，包含随机数据。数组的各个元素可以通过常规的索引表示法来访问，例如，a[2,4]。

操作 NumPy 数据的一个重要技巧是使用 切片。切片可以被认为是数组的一个子数组。例如，假设我们想要提取数组 a 中的中间两行和前两列的子数组。考虑以下命令：

b = a[2:4,0:2]
print b

现在，让我们做一个非常重要的观察。切片仅仅是数组的一种视图，并没有实际复制数据。可以通过运行以下命令来验证：

b[0,0]=0
print a

所以，b 的更改会影响数组 a！如果我们确实需要一个副本，我们需要明确表示我们想要一个副本。可以使用以下命令行来实现：

c = np.copy(a[2:4,0:2])
c[0,0] = -1
print a

在切片表示法 i:j 中，我们可以省略 i 或 j，在这种情况下，切片将表示对应轴的开始或结束：

print a[:4,3:]

省略 i 和 j 表示整个轴：

print a[:,2:4]

最后，我们可以使用 i:j:k 的表示法来指定切片中的步幅 k。在以下示例中，我们首先创建一个更大的随机数组来说明这一点：

a = np.random.rand(10,6)
print a
print
print a[1:7:2,5:0:-3]

现在，让我们考虑更高维度数组的切片。我们将通过创建一个非常大的三维数组来开始：

d1, d2, d3 = 4, 5, 3
a = np.random.rand(d1, d2, d3)
print a

假设我们想提取最后一个轴中索引为 1 的所有元素。这可以通过使用省略号对象轻松实现，示例如下：

print a[...,1]

上述命令行等同于以下命令：

print a[:,:,1]

在切片时，也可以沿着某一轴扩展矩阵，如下所示：

print a[0, :, np.newaxis, 0]

将前一个命令行的输出与以下输出进行比较：

print a[0, :, 0]