NumPy-基础知识-全-

NumPy 基础知识（全）

原文：Numpy Essentials

协议：CC BY-NC-SA 4.0

零、前言

无论您是科学/分析编程的新手还是经验丰富的专家，这本书都将为您提供成功创建，优化和分发 Python / NumPy 分析模块所需的技能。

从一开始，这本书将涵盖 NumPy 数组的关键功能以及调整数据格式以使其最适合您的分析需求的详细信息。 然后，您将获得各种多维，数据类型的分析所共有的核心和子模块的演练。 接下来，您将继续进行关键技术实现，例如线性代数和傅立叶分析。 最后，您将学习如何使用 Cython 和 NumPy C API 扩展 NumPy 的功能和性能。 本书的最后一章还提供了高级材料，可帮助您自己进一步学习。

如果您打算在分析项目中使用 NumPy，则本指南是非常宝贵的教程。

这本书涵盖的内容

第 1 章，“NumPy 简介”是本书的入门章节，它提供了帮助您设置环境的说明。 首先介绍科学 Python 模块系列（SciPy 栈），并说明 NumPy 在 Python 科学计算中所起的关键作用。

第 2 章，“NumPy ndarray对象”涵盖了 NumPy ndarray对象的基本用法，包括初始化，基本属性，数据类型和内存布局。 它还涵盖了操作下的理论，使您可以清晰地了解ndarray。

第 3 章，“使用 Numpy 数组”是有关 NumPy ndarray使用的高级章节，该章延续了第 2 章 NumPy ndarray对象。 它涵盖了 NumPy 中的通用函数，并向您展示了加快代码速度的技巧。 它还显示形状控制和广播规则。

第 4 章，“Numpy 核心和子模块”包括两个部分。 第一部分详细说明了 NumPy ndarray分配内存的方式与 CPU 缓存的交互之间的关系。 本章的第二部分介绍了特殊的 NumPy 数组，其中包含多种数据类型（结构/记录数组）。 此外，本章还将探讨 NumPy 中的实验性datetime64模块。

第 5 章，“NumPy 中的线性代数”从利用线性代数模块的矩阵和数学计算开始。 它向您展示了解决数学问题的多种方法：使用矩阵，向量分解和多项式。 它还提供了曲线拟合和回归的具体实践。

第 6 章，“NumPy 中的傅立叶分析”涵盖了使用 NumPy FFT 模块进行的信号处理以及傅立叶在放大信号/放大图像时的失真应用。 它还提供了 Python 中 matplotlib 包的基本用法。

第 7 章，“构建和分发 NumPy 代码”涵盖了有关使用 Python 打包和发布代码的基本细节。 它基本介绍了 NumPy 特定的安装文件以及如何构建扩展模块。

第 8 章，“使用 Cython 加速 NumPy” ，向用户介绍 Cython 编程语言，并向读者介绍可用于加速现有 Python 代码的技术 。

第 9 章，“NumPy C-API 简介”提供了对 NumPy C API 的基本介绍，以及一般而言，如何编写现有包装程序 C/C++ 库。 本章旨在提供一个简短的介绍，并为读者提供有关如何创建新包装程序和了解现有程序的基本知识。

第 10 章，“扩展阅读”是本书的最后一章。 它总结了我们在书中所学的内容，并探索了 4 个依赖于 NumPy 数组的 SciPy 栈的 Python 模块，这些模块为您提供了进一步科学 Python 编程的想法。

这本书需要什么

对于本书，您将需要以下设置：

• Python 2.x 或 3.x
• NumPy 1.9（或更高版本）
• IPython 笔记本
• Matplotlib 1.3（或更高版本）
• Windows 中的 gnu gcc 编译器或等效版本
• 设置工具

这本书适合谁

如果您了解 Python，但对科学编程还是陌生的，并且想进入科学计算的世界，或者您是具有分析经验的 Python 开发人员，但想获得洞察力来增强您的分析技能。 无论哪种情况，NumPy 或这本书都是您的理想选择。 学习 NumPy 以及如何将 NumPy 应用于 Python 程序非常适合作为构建专业分析应用的下一步。 熟悉基本的编程概念和数学会有所帮助，但是不需要任何先验经验。 后面的章节涵盖了诸如包分发，加速代码和 C/C++ 集成之类的概念，这些概念需要一定数量的编程和调试知识。 假定读者能够在其首选的 OS 中构建 C/C++ 程序（在 Linux 和 cygwin / migw 中使用 gcc，在 Windows 中使用 gcc）。

约定

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

文本，数据库表名称，文件夹名称，文件名，文件扩展名，路径名，虚拟 URL，用户输入和 Twitter 句柄中的代码字如下所示：“请注意，SciPy 可以表示很多东西，例如名为scipy。”

代码块设置如下：

In [42]: print("Hello, World!")

任何命令行输入或输出的编写方式如下：

In [6]: x
Out[6]:
array([[1, 2, 3],
       [2, 3, 4]])
In [7]: x[0,0]
Out[7]: 1
In [8]: x[1,2]
Out[8]: 4

新术语和重要词用粗体显示。

注意

警告或重要提示会出现在这样的框中。

提示

提示和技巧如下所示。

一、NumPy 简介

“我宁愿使用通用语言进行数学运算，也不愿尝试使用数学语言进行通用编程。”

-- John D Cook

在过去的十年中，Python 已成为科学计算中最受欢迎的编程语言之一。 其成功的原因很多，随着您着手本书，这些原因将逐渐变得明显。 与许多其他数学语言（例如 MATLAB，R 和 Mathematica）不同，Python 是一种通用编程语言。 因此，它为构建科学应用并将其进一步扩展到任何商业或学术领域提供了合适的框架。 例如，考虑一个（某种）简单的应用，该应用要求您编写软件并预测博客文章的受欢迎程度。 通常，这些是您要执行此操作的步骤：

1. 生成博客文章及其相应评级的语料库（假设此处的评级可适当量化）。
2. 制定一个模型，该模型根据与博客文章相关的内容和其他数据生成评分。
3. 根据您在步骤 1 中找到的数据来训练模型。请继续进行此操作，直到您对模型的可靠性充满信心为止。
4. 将模型部署为 Web 服务。

通常，在执行这些步骤时，您会发现自己在不同的软件栈之间跳转。 第 1 步需要进行大量的网页抓取。 Web 抓取是一个非常普遍的问题，几乎每种编程语言都有一些工具可以抓取 Web（如果您已经在使用 Python，则可能会选择 BeautifulSoup 或 Scrapy）。 第 2 步和第 3 步涉及解决机器学习问题，并且需要使用复杂的数学语言或框架，例如 Weka 或 MATLAB，这只是提供机器学习功能的众多工具中的少数几种。 同样，步骤 4 可以使用许多不同的工具以多种方式实现。 没有一个正确的答案。 由于这个问题已经被许多科学家和软件开发人员充分研究和解决（在一定程度上），因此，找到可行的解决方案并不困难。 但是，诸如稳定性和可伸缩性之类的问题可能会严重限制您在问题的每个步骤中对编程语言，Web 框架或机器学习算法的选择。 这就是 Python 胜过大多数其他编程语言的地方。 前面的所有步骤（以及更多步骤）都只能使用 Python 和一些第三方 Python 库来完成。 这种使用 Python 开发软件的灵活性和便捷性正是使其成为科学计算生态系统的理想宿主。 在 Ivan Idris 所写的《Python 数据分析》中可以找到关于 Python 作为成熟的应用开发语言的非常有趣的解释。精确地讲，Python 是一种用于快速原型制作的语言，并且由于其随着时间的推移而获得了广泛的科学生态系统，它也被用于构建生产质量的软件。 这个生态系统的基础是 NumPy。

数字 Python（NumPy）是 Numeric 包的后续产品。 它最初由 Travis Oliphant 编写，是 Python 科学计算环境的基础。 它在 2005 年初从功能更广泛的 SciPy 模块分支出来，并于 2006 年中首次稳定发布。 从那以后，它在从事数学，科学和工程领域的 Python 爱好者中越来越受欢迎。 本书的目的是使您对 NumPy 足够熟悉，以便您能够使用它并使用它构建复杂的科学应用。

Python 科学栈

让我们首先简要浏览一下 Python 科学计算（SciPy）栈。

注意

请注意，SciPy 可能意味着很多事情：名为 scipy 的 Python 模块，整个 SciPy 栈，或在世界各地举行的有关科学 Python 的三个会议中的任何一个。

Figure 1: The SciPy stack, standard, and extended libraries

IPython 的主要作者 Fernando Perez 在 2012 年加拿大 PyCon 的主题演讲中说：

“科学计算的发展不仅仅是因为软件的发展，而且还因为我们作为科学家所做的不仅仅是浮点运算。”

这正是 SciPy 栈拥有如此丰富的功能的原因。 大多数 SciPy 栈的演进是由试图以通用编程语言解决科学和工程问题的科学家和工程师团队推动的。 对 NumPy 为什么如此重要的一个单方面的解释是，它提供了科学计算中大多数任务所必需的核心多维数组对象。 这就是为什么它是 SciPy 栈的根本原因。 NumPy 使用久经考验的科学库，提供了一种简单的方法来与遗留的 Fortran 和 C/C++ 数字代码对接，我们知道该库已经运行了数十年。 全世界的公司和实验室都使用 Python 将已经存在很长时间的遗留代码粘合在一起。 简而言之，这意味着 NumPy 允许我们站在巨人的肩膀上。 我们不必重新发明轮子。 这是每个其他 SciPy 包的依赖项。 NumPy ndarray对象实际上是下一章的主题，它是 Pythonic 接口，用于用 Fortran，C 和 C++ 编写的库所使用的数据结构。 实际上，NumPy ndarray对象使用的内部内存布局实现 C 和 Fortran 布局。 这将在以后的章节中详细讨论。

栈的下一层包括 SciPy，matplotlib，IPython（Python 的交互式外壳；我们将在整本书中将其用作示例，其安装和使用的详细信息将在后面的部分中提供）以及 SymPy 模块。 SciPy 提供了生态系统主要部分所依赖的大部分科学和数值功能。 Matplotlib 是 Python 中的事实绘图和数据可视化库。 IPython 是用于 Python 中科学计算的日益流行的交互式环境。 实际上，该项目已经进行了如此积极的开发并享有很高的知名度，以至于它不再局限于 Python，而且将其功能扩展到其他科学语言，尤其是 R 和 Julia。 栈中的这一层可以看作是 NumPy 的面向核心数组的功能与栈较高层提供的特定于域的抽象之间的桥梁。 这些特定于领域的工具通常称为 SciKits，受欢迎的工具包括 scikit-image（图像处理），scikit-learn（机器学习），statsmodels（统计信息），pandas（高级数据分析）等等。 由于科学的 Python 社区非常活跃，因此几乎不可能在 Python 中列出每个科学包，并且针对大量的科学问题总是有很多发展。 跟踪项目的最佳方法是参与社区。 加入邮件列表，编写代码，将软件用于日常计算需求并报告错误，这非常有用。 本书的目标之一是使您足够感兴趣，以积极地参与科学的 Python 社区。

NumPy 数组的必要性

初学者提出的一个基本问题是。 为什么数组对于科学计算完全必要？ 当然，可以对任何抽象数据类型（如列表）执行复杂的数学运算。 答案在于数组的众多属性，这些属性使它们明显更有用。 在本节中，让我们看一下其中的一些属性，以强调为什么诸如 NumPy ndarray对象之类的东西根本不存在。

表示矩阵和向量

矩阵和向量的抽象数学概念是许多科学问题的核心。 数组为这些概念提供了直接的语义链接。 确实，每当一本数学文献提到矩阵时，就可以安全地将数组视为代表矩阵的软件抽象。 在科学文献中，A[ij]等表达式通常用于表示数组A的第i行和j列的元素。 NumPy 中的相应表达将简单地是A[i, j]。 对于矩阵运算，NumPy 数组还支持向量化（有关详细信息，请参见第 3 章，“使用 Numpy 数组”），这大大加快了执行速度。 向量化使代码更简洁，更易于阅读，并且更类似于数学符号。 像矩阵一样，数组也可以是多维的。 数组的每个元素都可以通过一组称为索引的整数来寻址，而使用整数集访问数组的元素的过程称为索引。 确实可以在不使用数组的情况下实现此功能，但这将很麻烦并且非常不必要。

效率

效率在软件中可能意味着很多事情。 该术语可用于指代程序的执行速度，其数据检索和存储性能，其内存开销（程序执行时消耗的内存）或其整体吞吐量。 就几乎所有这些特性而言，NumPy 数组都比大多数其他数据结构要好（只有少数例外，例如 pandas，DataFrame或 SciPy 的稀疏矩阵，我们将在后面的章节中介绍）。 由于 NumPy 数组是静态类型且同质的，因此可以用编译语言实现快速数学运算（默认实现使用 C 和 Fortran）。 效率（在同类数组上运行快速算法的可用性）使 NumPy 变得流行且重要。

易于开发

NumPy 模块是用于数学任务的现成功能的强大平台。 它极大地增加了 Python 的开发难度。 以下是该模块包含的内容的简要概述，我们将在本书中探讨其中的大部分内容。 有关 NumPy 模块的详细介绍，请参见 Travis Oliphat 所写的权威的《NumPy 指南》。 NumPy API 非常灵活，以至于科学 Python 社区已广泛采用它作为构建科学应用的标准 API。 可以参考 Van Der Walt 等人所写的《NumPy 数组：有效数值计算的结构》：

子模块	内容
numpy.core	基本对象
lib	其他工具
linalg	基本线性代数
fft	离散傅立叶变换
random	随机数生成器
distutils	增强的构建和发行版
testing	单元测试
f2py	Fortran 代码的自动包装

学术界和工业界的 NumPy

有人说，如果您在时代广场站足够长的时间，就会遇到世界上每个人。 现在，您必须已经确信 NumPy 是 SciPy 的时代广场。 如果您使用 Python 编写科学应用，那么不需深入研究 NumPy，您将无能为力。 图 2 显示了不同抽象级别的科学计算中 SciPy 的范围。 红色箭头表示科学软件应具有的各种低级功能，蓝色箭头表示利用这些功能的不同应用领域。 配备 SciPy 栈的 Python 处于提供这些功能的语言的最前沿。

Google 学术搜索 NumPy 会返回近 6,280 个结果。 其中一些是有关 NumPy 和 SciPy 栈本身的论文和文章，还有许多是有关 NumPy 在各种研究问题中的应用的。 学者们喜欢 Python，SciPy 栈作为世界上无数大学和研究实验室中科学编程的主要语言越来越受欢迎，这说明了 Python。 许多科学家和软件专业人员的经验已发布在 Python 网站上：

Figure 2: Python versus other languages

书中使用的代码约定

现在，Python 和 NumPy 的信誉已经确立，让我们动手吧。

本书中所有 Python 代码使用的默认环境是 IPython。 下一节将介绍如何安装 IPython 和其他工具。 在整本书中，您只需在命令窗口或 IPython 提示符下输入输入即可。 除非另有说明，否则code将引用 Python 代码，command将引用 bash 或 DOS 命令。

所有 Python 输入代码都将按照以下代码段格式进行格式化：

 In [42]: print("Hello, World!")

上一小段中的In [42]:表示这是第 42 个 IPython 会话输入。 同样，所有命令行输入的格式将如下：

 $ python hello_world.py 

在 Windows 系统上，相同的命令如下所示：

C:\Users\JohnDoe> python hello_world.py

为了保持一致，无论操作系统如何，$符号都将用于表示命令行提示符。 诸如C:\Users\JohnDoe>之类的提示将不会出现在书中。 通常，$符号表示 Unix 系统上的 bash 提示，但是相同的命令（无需键入实际的美元符号或任何其他字符）也可以在 Windows 上使用。 但是，如果您使用的是 Cygwin 或 Git Bash，则也应该能够在 Windows 上使用 Bash 命令。

请注意，如果您在 Windows 上安装 Git，则默认情况下 Git Bash 可用。

安装要求

让我们看一下在继续之前需要设置的各种要求。

使用 Python 发行版

本书所需的三个最重要的 Python 模块是 NumPy，IPython 和 matplotlib。 在本书中，代码基于 Python 3.4 / 2.7 兼容版本，NumPy 1.9 版本和 matplotlib 1.4.3。 安装这些要求（甚至更多）的最简单方法是安装完整的 Python 发行版，例如 Enthought Canopy，EPD，Anaconda 或 PythonXY。 一旦安装了其中任何一个，就可以安全地跳过本节的其余部分，并且应该可以开始了。

注意

面向 Canopy 用户的注意事项：可以使用 Canopy GUI，该 GUI 包括嵌入式 IPython 控制台，文本编辑器和 IPython Notebook 编辑器。 使用命令行时，为了获得最佳效果，请使用 Canopy 的“工具”菜单中的 Canopy 终端。

Windows OS 用户注意事项：除了 Python 发行版，您还可以从 Ghristoph Gohlke 网站安装预构建的 Windows python 扩展包。

使用 Python 包管理器

您还可以使用以下命令之一使用 Python 包管理器（例如 enpkg，Conda，PIP 或 EasyInstall）来安装需求。 将numpy替换为您要安装的其他任何包名称，例如ipython，matplotlib等：

$ pip install numpy
$ easy_install numpy
$ enpkg numpy # for Canopy users
$ conda install numpy # for Anaconda users

使用本机包管理器

如果要使用的 Python 解释器随操作系统一起提供，而不是第三方安装，则您可能更喜欢使用特定于操作系统的包管理器，例如 aptitude，yum 或 Homebrew。 下表说明了包管理器和用于安装 NumPy 的相应命令：

包管理器	命令
apt	$ sudo apt-get install python-numpy
yum	$ yum install python-numpy
brew	$ brew install numpy

请注意，在带有 Homebrew 的 OS X 系统上安装 NumPy（或任何其他 Python 模块）时，Python 最初应与 Homebrew 一起安装。

有关详细的安装说明，请参见 NumPy，IPython 和 matplotlib 的相应网站。 作为预防措施，要检查 NumPy 是否已正确安装，请打开 IPython 终端并键入以下命令：

 In [1]: import numpy as np 
 In [2]: np.test()

如果第一个语句看起来什么也没做，则表明这是一个好兆头。 如果执行时没有任何输出，则表示已安装 NumPy 并将其正确导入到 Python 会话中。 第二条语句运行 NumPy 测试套件。 这不是绝对必要的，但永远不要太谨慎。 理想情况下，它应运行几分钟并产生测试结果。 它可能会生成一些警告，但是这些都不会引起警报。 如果您愿意，也可以运行 IPython 和 matplotlib 的测试套件。

注意

请注意，只有从源代码安装了 matplotlib 时，matplotlib 测试套件才能可靠运行。 但是，测试 matplotlib 并不是非常必要。 如果您可以导入 matplotlib 而没有任何错误，则表明它可以使用了。

恭喜你！ 我们现在准备开始。

总结

在本章中，我们向自己介绍了 NumPy 模块。 我们了解了 NumPy 是如何为从事科学计算工作的人们提供的有用软件工具。 我们安装了完成本书其余部分所需的软件。

在下一章中，我们将介绍功能强大的 NumPy ndarray对象，向您展示如何有效地使用它。

十、扩展阅读

NumPy 是 Python 中功能强大的科学模块； 希望在前九章中，我们已经向您展示了足以向您证明这一点的内容。 ndarray是所有其他 Python 科学模块的核心。 使用 NumPy 的最佳方法是使用numpy.ndarray作为基本数据格式，并将其与其他科学模块组合以进行预处理，分析，计算，导出等。 在本章中，我们的重点是向您介绍可以与 NumPy 一起使用的两个模块，并使您的工作/研究效率更高。

在本章中，我们将讨论以下主题：

• Pandas
• scikit-learn
• netCDF4
• SciPy

Pandas

到目前为止，pandas 是 Python 中最可取的数据预处理模块。 它处理数据的方式与 R 非常相似。它的数据框不仅为您提供视觉上吸引人的表打印输出，而且还允许您以更直观的方式访问数据。 如果您不熟悉 R，请尝试考虑以编程方式使用电子表格软件（例如 Microsoft Excel 或 SQL 表）。 这涵盖了 Pandas 所做的很多事情。

您可以从 Pandas 官方网站下载并安装 Pandas。 一种更可取的方法是使用 pip 或安装 Python 科学发行版，例如 Anaconda。

还记得我们如何使用numpy.genfromtxt()读取第 4 章， “Numpy 核心和子模块”中的csv数据吗？ 实际上，使用 Pandas 来读取表格并将经过预处理的数据传递给ndarray（简单地执行np.array(data_frame)会将数据帧转换为多维ndarray）对于分析来说是更可取的工作流程。 在本节中，我们将向您展示 Pandas 的两个基本数据结构：Series（用于一维）和DataFrame（用于二维或多维）。然后，我们将向您展示如何使用 Pandas 来读取表并将数据传递给它。

然后，我们将向您展示如何使用 Pandas 读取表并将数据传递给ndarray进行进一步分析。 让我们从pandas.Series开始：

In [1]: import pandas as pd 
In [2]: py_list = [3, 8, 15, 25, 11] 
In [3]: series = pd.Series(py_list) 
In [4]: series 
Out[4]: 
0     3 
1     8
2    15
3    25 
4    11 
dtype: int64 

在前面的示例中，您可以看到我们已经将 Python 列表转换为 Pandasseries，并且在打印series时，这些值完美对齐并具有与之关联的索引号（0 到 4）。 当然，我们可以指定自己的索引（以1开头或以字母形式）。 看下面的代码示例：

In [5]: indices = ['A', 'B', 'C', 'D', 'E'] 
In [6]: series = pd.Series(py_list, index = indices) 
In [7]: series 
Out[7]: 
A     3 
B     8 
C    15 
D    25 
E    11 
dtype: int64 

我们将索引从数字更改为从A-E的字母。 更方便的是，当我们将 Python 词典转换为 Pandas Series时，执行此操作所需的键将自动成为索引。 尝试练习转换字典。 接下来，我们将探索DataFrame，这是 Pandas 中最常用的数据结构：

In [8]: data = {'Name': ['Brian', 'George', 'Kate', 'Amy', 'Joe'], 
   ...:         'Age': [23, 41, 26, 19, 35]} 
In [9]: data_frame = pd.DataFrame(data) 
In [10]: data_frame 
Out[10]: 
   Age    Name 
0   23   Brian 
1   41  George
2   26    Kate 
3   19     Amy 
4   35     Joe 

在前面的示例中，我们创建了DataFrame，其中包含两列：第一个是Name，第二个是Age。 您可以从打印输出中看到它看起来像一张表格，因为它的格式正确。 当然，您也可以更改数据帧的索引。 但是，数据帧的优势远不止于此。 我们可以访问每列中的数据或对其进行排序（按列名，访问data_frame.column_name或data_frame[column_name]需要两个符号）； 我们甚至可以分析汇总统计信息。 为此，请看以下代码示例：

In [11]: data_frame = pd.DataFrame(data) 
In [12]: data_frame['Age'] 
Out[12]: 
0    23 
1    41 
2    26 
3    19 
4    35 
Name: Age, dtype: int64 
In [13]: data_frame.sort(columns = 'Age') 
Out[13]: 
   Age    Name 
3   19     Amy 
0   23   Brian 
2   26    Kate 
4   35     Joe 
41  George 
In [14]: data_frame.describe() 
Out[14]: 
             Age 
count   5.000000 
mean   28.800000 
std     9.011104 
min    19.000000 
25%    23.000000 
50%    26.000000 
75%    35.000000 
max    41.000000 

在前面的示例中，我们仅获得Age列，并按Age对DataFrame进行排序。 当我们使用describe()时，它将计算所有数字字段的摘要统计信息（包括计数，均值，标准差，最小值，最大值和百分位数）。在本节的最后部分，我们将使用 Pandas 读取

在本节的最后部分，我们将使用 Pandas 读取csv文件并将一个字段值传递给ndarray以进行进一步的计算。 example.csv文件来自国家统计局（ONS）。 请访问这里了解更多详细信息。 我们将在 ONS 网站上使用Sale counts by dwelling type and local authority, England and Wales（房屋类型和地方当局（英格兰和威尔士）的销售计数）。 您可以按主题名称搜索它，以访问下载页面或选择您感兴趣的任何数据集。在以下示例中，我们将示例数据集重命名为sales.csv：

In [15]: sales = pd.read_csv('sales.csv') 
In [16]: sales.shape 
Out[16]: (348, 97) 
In [17]: sales.columns[:3] 
Out[17]: Index([u'LA_Code', u'LA_Name', u'1995_COUNT_ALL_TYPES'], dtype='object') 
In [18]: sales['1995_COUNT_ALL_TYPES'].head() 
Out[18]: 
0    1,188 
1    1,652 
2    1,684 
3    2,314 
4    1,558 
Name: 1995_COUNT_ALL_TYPES, dtype: object 

首先，我们将sale.csv读入名为sales的DataFrame对象; 当我们打印销售的shapes时，我们发现数据框中有 384 条记录和 97 列。 DataFrame column属性的返回列表是一个普通的 Python 列表，我们在数据中打印了前三列：LA_Code，LA_Name和1995_COUNT_ALL_TYPES。 然后，我们使用head()函数在1995_COUNT_ALL_TYPES中打印了前五个记录（tail()函数将打印后五个记录）。

同样，pandas 是 Python 中一个功能强大的预处理模块（通常，其数据处理能力超过其预处理功能，但在前面的示例中，我们仅介绍了预处理部分），并且它具有许多方便的功能来帮助您清理数据并准备数据。 您的分析。 本节仅作介绍； 由于空间限制，我们无法涵盖很多内容，例如数据透视，datetime等。 希望您能理解并开始提高脚本的效率。

scikit-learn

Scikit 是 SciPy 工具包的缩写，它是 SciPy 的附加包。 它提供了广泛的分析模块，而 scikit-learn 是其中之一。 这是迄今为止最全面的 Python 机器学习模块。 scikit-learn 提供了一种简单有效的方法来执行数据挖掘和数据分析，并且它拥有非常活跃的用户社区。

您可以从 scikit-learn 的官方网站下载并安装它。 如果您使用的是 Python 科学发行版（例如 Anaconda），则也包含在其中。

现在，是时候使用 scikit-learn 进行一些机器学习了。 scikit-learn 的优点之一是它提供了一些用于实践的样本数据集（演示数据集）。 让我们首先加载糖尿病数据集。

In [1]: from sklearn.datasets import load_diabetes 
In [2]: diabetes = load_diabetes() 
In [3]: diabetes.data 
Out[3]: 
array([[ 0.03807591,  0.05068012,  0.06169621, ..., -0.00259226, 
         0.01990842, -0.01764613], 
       [-0.00188202, -0.04464164, -0.05147406, ..., -0.03949338, 
        -0.06832974, -0.09220405], 
       [ 0.08529891,  0.05068012,  0.04445121, ..., -0.00259226, 
         0.00286377, -0.02593034], 
       ..., 
       [ 0.04170844,  0.05068012, -0.01590626, ..., -0.01107952, 
        -0.04687948,  0.01549073], 
       [-0.04547248, -0.04464164,  0.03906215, ...,  0.02655962, 
         0.04452837, -0.02593034], 
       [-0.04547248, -0.04464164, -0.0730303 , ..., -0.03949338, 
        -0.00421986,  0.00306441]]) 
In [4]: diabetes.data.shape 
Out[4]: (442, 10) 

我们从sklearn.datasets中加载了一个名为diabetes的样本数据集； 它包含 442 个观测值，10 个维度，范围从 -2 到 2。 Toy数据集还提供了标记数据用于监督学习（如果您不熟悉机器学习，请尝试将标记数据视为类别）。 在我们的示例中，可以从diabetes.target调用diabetes数据集中的标记数据，范围为 25 到 346。

还记得我们在第 5 章， “numpy 中的线性代数”中如何进行线性回归吗？ 我们将使用 scikit-learn 再次执行它。 同样，我建议您在开发脚本以帮助您进行研究或分析时，请使用 NumPy ndarray作为常规数据格式； 但是，对于计算，使用 scipy，scikit-learn 或其他科学模块会更好。 机器学习的优势之一是模型评估（您可以在其中训练和测试结果）。 使用此方法，我们将数据分为两个数据集：训练数据集和测试数据集，然后将这两个数据集传递给线性回归：

In [5]: from sklearn.cross_validation import train_test_split 
In [6]: X_train, X_test, y_train, y_test  = 
        train_test_split(diabetes.data, 
diabetes.target, 
random_state = 50) 

在前面的示例中，我们使用train_test_split()函数将糖尿病数据集分为训练和测试数据集（针对数据及其类别）。 前两个参数是我们要拆分的数组。 random_state参数是可选的，这意味着伪随机数生成器状态用于随机采样。 默认拆分率是 0.25，这意味着 75% 的数据拆分为训练集，而 25% 的数据拆分为测试集。 您可以尝试打印出我们刚刚创建的训练/测试数据集，以查看其分布情况（在前面的代码示例中，X_train代表糖尿病数据的训练数据集，X_test代表糖尿病测试数据，y_train代表分类的糖尿病训练数据，y_test代表分类的糖尿病测试数据。

接下来，我们将数据集拟合为线性回归模型：

In [7]: from sklearn.linear_model import LinearRegression 
In [8]: lr = LinearRegression() 
In [9]: lr.fit(X_train, y_train) 
Out[9]: LinearRegression(copy_X = True, fit_intercept = True, 
     Normalize = False) 
In [10]: lr.coef_ 
Out[10]: 
array([  80.73490856, -195.84197988,  474.68083473,  371.06688824, 
       -952.26675602,  611.63783483,  174.40777144,  159.78382579, 
        832.01569658,   12.04749505]) 

首先，我们从sklearn.linear_model创建了一个LinearRegression对象，并使用fit()函数拟合了X_train和y_train数据集。 我们可以通过调用coef_属性来检查线性回归的估计系数。 此外，我们可以使用拟合的线性回归进行预测。 看下面的例子：

In [11]: lr.predict(X_test)[:10] 
Out[11]: 
array([  71.96974998,   82.55916305,  265.71560021,   79.37396336, 
         72.48674613,   47.01580194,  149.11263906,  185.36563936, 
         94.88688296,  132.08984366]) 

predict()函数用于基于我们拟合训练数据集的线性回归来预测测试数据集； 在前面的示例中，我们打印了前 10 个预测值。 这是y的预测值和测试值的图：

In [12]: lr.score(X_test, y_test) 
Out[12]: 0.48699089712593369 

然后，我们可以使用score()函数检查预测的确定 R 平方。

这几乎是 scikit-learn 中的标准拟合和预测过程，并且非常直观且易于使用。 当然，除了回归之外，scikit-learn 还可以执行许多分析，例如分类，聚类和建模。 希望本节对您的脚本有所帮助。

netCDF4

netCDF4 是 netCDF 库的第四个版本，该库是在 HDF5（分层数据格式，旨在存储和组织大量数据）的基础上实现的，从而可以管理非常大和复杂的多维数据。 netCDF4 的最大优点是，它是一种完全可移植的文件格式，对集合中数据对象的数量或大小没有限制，并且在可归档的同时也可以追加。 许多科研组织将其用于数据存储。 Python 还具有访问和创建此类数据格式的接口。

您可以从官方文档页面，或从这里下载并安装该模块。 它不包含在标准的 Python 科学发行版中，但已内置在 NumPy 中，可以与 Cython 一起构建（建议但并非必需）。

对于以下示例，我们将使用 Unidata 网站上的示例netCDF4文件，该文件位于这里，并且我们将以气候系统模型为例：sresa1b_ncar_ccsm3-example.nc

首先，我们将使用netCDF4模块稍微探索一下数据集，并提取我们需要进行进一步分析的值：

In [1]: import netCDF4 as nc 
In [2]: dataset = nc.Dataset('sresa1b_ncar_ccsm3-example.nc', 'r') 
In [3]: variables = [var for var in dataset.variables] 
In [4]: variables 
Out[4]: 
['area', 'lat', 'lat_bnds', 'lon', 'lon_bnds', 'msk_rgn', 'plev', 
'pr', 'tas', 'time', 'time_bnds', 'ua'] 

我们导入了 python netCDF4模块，并使用Dataset()函数读取了示例netCDF4文件。 r参数表示文件处于只读模式，因此当我们要附加文件或w创建新文件时，我们也可以指定a。 然后，我们获得了存储在数据集中的所有变量，并将它们保存到名为变量的列表中（请注意，variables属性将返回变量对象的 Python 字典）。 最后，我们使用以下命令在数据集中打印出变量：

In [5]: precipitation = dataset.variables['pr'] 
In [6]: precipitation.standard_name 
Out[6]: 'precipitation_flux' 
In [7]: precipitation.missing_value 
Out[7]: 1e+20 
In [8]: precipitation.ndim 
Out[8]: 3 
In [9]: precipitation.shape 
Out[9]: (1, 128, 256) 
In [10]: precipitation[:, 1, :10] 
Out[10]: 
array([[  8.50919207e-07,   8.01471970e-07,   7.74396426e-07, 
          7.74230614e-07,   7.47181844e-07,   7.21426375e-07, 
          7.19294349e-07,   6.99790974e-07,   6.83397502e-07, 
          6.74683179e-07]], dtype=float32) 

在前面的示例中，我们选择了一个名为pr的变量并将其保存到precipitation中。 众所周知netCDF4是一种自我描述的文件格式; 您可以创建和访问存储在变量中的任何用户定义属性，尽管最常见的是standard_name，它告诉我们该变量代表降水通量。 我们检查了另一个常用属性missing_value，该属性表示存储在netCDF4文件中的无数据值。 然后，我们通过ndim来打印降水量的维数，并通过shape属性来打印形状。 最后，我们要获取第 1 行的值，即netCDF4文件中的前 10 列; 为此，只需像往常一样使用索引。

接下来，我们将介绍创建netCDF4文件并将三维 NumPy ndarray作为变量存储的基础知识：

In [11]: import numpy as np 
In [12]: time = np.arange(10) 
In [13]: lat = 54 + np.random.randn(8) 
In [14]: lon = np.random.randn(6) 
In [15]: data = np.random.randn(480).reshape(10, 8, 6) 

首先，我们准备了一个三维ndarray（数据）以存储在netCDF4文件中； 数据建立在三个维度中，分别是时间（time，大小为 10），纬度（lat，大小为 8）和经度（lon，大小为 6）。 在netCDF4中，时间不是datetime对象，而是从定义的开始时间（在unit属性中指定）开始的时间单位数（可以是秒，小时，天等）。 稍后再向您解释）。 现在，我们拥有了要存储在文件中的所有数据，因此让我们构建 netCDF 结构：

In [16]: output = nc.Dataset('test_output.nc', 'w') 
In [17]: output.createDimension('time', 10) 
In [18]: output.createDimension('lat', 8) 
In [19]: output.createDimension('lon', 6) 
In [20]: time_var = output.createVariable('time', 'f4', ('time',)) 
In [21]: time_var[:] = time 
In [22]: lat_var = output.createVariable('lat', 'f4', ('lat',)) 
In [23]: lat_var[:] = lat 
In [24]: lon_var = output.createVariable('lon', 'f4', ('lon',)) 
In [25]: lon_var[:] = lon 

我们通过指定文件路径并使用w写入模式来初始化netCDF4文件。 然后，我们使用createDimension()构建结构以指定尺寸：time，lat和lon。 每个尺寸都有一个变量来表示其值，就像轴的刻度一样。 接下来，我们将把三维数据保存到文件中：

In [26]: var = output.createVariable('test', 'f8', ('time', 'lat', 'lon')) 
In [27]: var[:] = data 

变量的创建始终从createVariable()函数开始，并指定变量名称，变量数据类型以及与其关联的维。 第二步是将ndarray 的相同形状传递到声明的变量中。 现在我们已经将整个数据存储在文件中，我们可以指定属性以帮助描述数据集。 以下示例使用time变量说明如何指定属性：

In [28]: time_var.standard_name = 'Time' 
In [29]: time_var.units = 'days since 2015-01-01 00:00:00' 
In [30]: time_var.calendar = 'gregorian' 

因此，现在时间变量已与单位和日历相关联，因此ndarray时间将根据我们指定的单位和日历转换为日期； 这类似于所有变量。 完成netCDF4文件的创建后，最后一步是关闭文件连接：

In [31]: output.close() 

上面的代码向您展示了 Python netCDF4 API 的用法，以便读取和创建netCDF4文件。 该模块不包含任何科学计算（因此它不包含在任何 Python 科学发行版中），但是目标位于文件 I/O 的接口中，该接口可以是研究和分析的第一阶段或最后阶段。

SciPy

SciPy 是一个着名的 Python 库，专注于科学计算（它包含用于优化，线性代数，积分，内插以及诸如 FFT，信号和图像处理等特殊功能的模块）。 它建立在 NumPy 数组对象的基础上，并且 NumPy 是整个 SciPy 栈的一部分（请记住，我们在第 1 章， “NumPy 简介”）。 但是，SciPy 模块包含各种主题，而我们不能仅在一个部分中进行介绍。 让我们看一个图像处理（降噪）示例，以帮助您了解 SciPy 可以做什么：

In [1]: from scipy.misc import imread, imsave, ascent 
In [2]: import matplotlib.pyplot as plt 
In [3]: image_data = ascent() 

首先，我们从 SciPy 的其他例程中导入三个函数：imread，imsave和ascent。 在下面的示例中，我们使用内置图像ascent，它是512 x 512灰度图像。 当然，您可以使用自己的图像。 只需调用imread('your_image_name')，它将作为ndarray加载。

我们在此处导入的matplotlib模块的pyplot结果仅用于显示图像； 我们在第 6 章，“在 NumPy 中进行傅里叶分析”这样做了。 这是内置图像ascent：

现在，我们可以为图像添加一些噪点，并使用pyplot模块显示噪点图像：

In [4]: import numpy as np 
In [5]:noise_img = image_data +  image_data.std() * np.random.random(image_data.shape) 
In [6]: imsave('noise_img.png', noise_img) 
In [7]: plt.imshow(noise_img) 
Out[7]: <matplotlib.image.AxesImage at 0x20066572898> 
In [8]: plt.show() 

在此代码段中，我们导入numpy以根据图像形状生成一些随机噪声。 然后，将噪点图像保存到noise_img.png。 噪点图像如下所示：

接下来，我们将使用 SciPy scipy.ndimage中的多维图像处理模块对噪波图像应用滤镜以使其平滑。 ndimage模块提供各种过滤器； 有关详细列表，请参考这里，但是在以下示例中，我们将仅使用高斯和均匀滤波器：

In [9]: from scipy import ndimage 
In [10]: gaussian_denoised = ndimage.gaussian_filter(noise_img, 3) 
In [11]: imsave('gaussian_denoised.png', gaussian_denoised ) 
In [12]: plt.imshow(gaussian_denoised) 
Out[12]: <matplotlib.image.AxesImage at 0x2006ba54860> 
In [13]: plt.show() 
In [14]: uniform_denoised = ndimage.uniform_filter(noise_img) 
In [15]: imsave('uniform_denoised.png', uniform_denoised) 
In [16]: plt.imshow(gaussian_denoised) 
Out[17]: <matplotlib.image.AxesImage at 0x2006ba80320> 
In [18]: plt.show() 

首先，我们从 SciPy 导入ndimage，在noise_image上应用高斯滤波器，将sigma（高斯内核的标准偏差）设置为3，然后将其保存到gaussian_denoised.png。 查看上一张图片的左侧。 通常，sigma越大，图像将越平滑，这意味着细节丢失。 我们应用的第二个过滤器是均匀过滤器，并采用了所有参数的默认值，这将导致上一张图像的右侧。 尽管统一滤波器保留了原始图像的更多细节，但图像仍包含噪点。

上一个示例是使用 SciPy 的简单图像处理示例。 但是，SciPy 不仅可以处理图像处理，还可以执行许多类型的分析/科学计算。 有关详细信息，请参阅《SciPy 数值和科学计算学习手册第二版》，Packt Publishing。

总结

NumPy 当然是使用 Python 进行科学计算的核心：许多模块都基于 NumPy。 尽管有时您可能会发现 NumPy 没有分析模块，但它无疑为您提供了一种接触广泛科学模块的方法。

我们希望本书的最后一章为您提供了一个关于将这些模块与 NumPy 一起使用的好主意，并使您的脚本更加有效（本书中无法涵盖很多便捷的 NumPy 模块；仅在 GitHub 或 PyPI 上度过一个下午，您可能会发现其中的少数几个）。 最后但并非最不重要的一点，感谢您花时间与我们一起完成许多功能。 现在与 NumPy 一起玩吧！

二、NumPy ndarray对象

面向数组的计算是计算科学的核心。 这是大多数 Python 程序员都不习惯的。 尽管列表或字典的理解是相对于数组的，有时与数组的用法类似，但是在性能和​​操作上，列表/字典和数组之间还是存在巨大差异。 本章介绍了 NumPy 中的基本数组对象。 它涵盖了可以从 NumPy 数组的固有特性中收集的信息，而无需对该数组执行任何外部操作。

本章将涉及的主题如下：

• numpy.ndarray以及如何使用它-面向基本数组的计算
• numpy.ndarray内存访问，存储和检索的性能
• 索引，切片，视图和副本
• 数组数据类型

numpy.ndarray入门

在本节中，我们将介绍 numpy ndarray的一些内部结构，包括其结构和行为。 开始吧。 在 IPython 提示符中输入以下语句：

In [1]: import numpy as np 

In [2]: x = np.array([[1,2,3],[2,3,4]]) 

In [3]: print(x)

NumPy 与其他模块（例如 Python 标准库中的math模块）中的函数共享其函数名称。 不建议使用如下所示的导入：

from numpy import * 

因为它可能会覆盖全局名称空间中已经存在的许多函数，所以不建议这样做。 这可能会导致您的代码出现意外行为，并可能在其中引入非常细小的错误。 这也可能会在代码本身中造成冲突（例如 numPy 具有any并会与系统any关键字发生冲突），并可能在检查或调试一段代码时引起混乱。 因此，重要的是建议始终使用带有显式名称的导入 numPy，例如第一行中使用的np约定：-import numpy as np，这是用于导入目的的标准约定，因为它有助于开发人员找出函数的来源。 这可以避免大型程序中的许多混乱。

如我们将看到的，可以用多种方式创建 NumPy 数组。 创建数组的最简单方法之一是使用array函数。 注意，我们向函数传递了一个列表列表，组成列表的长度相等。 每个组成列表成为数组中的一行，并且这些列表的元素填充了结果数组的列。 array函数可以在列表甚至嵌套列表上调用。 由于此处输入的嵌套级别是 2，因此生成的数组是二维的。 这意味着可以使用两个整数集对数组进行索引。 计算数组维数的最简单方法是检查数组的ndim属性：

In [4]: x.ndim 

Out [4]: 2 

这也可以通过检查数组的shape属性以其他（间接）方式来实现。 数组的维数将等于您在shape属性中看到的数字。 （但是请注意，这并不是shape属性的目的。）：

In [5]: x.shape 

Out [5]: (2, 3) 

这意味着该数组具有两行三列。 重要的是要注意，与 MATLAB 和 R 不同，NumPy 数组的索引是从零开始的。 也就是说，NumPy 数组的第一个元素索引为零，而最后一个元素索引为整数n-1，其中n是数组沿相应维度的长度。 因此，对于我们刚刚创建的数组，可以使用一对零来访问数组左上角的元素，并且可以使用索引来访问右下角的元素（1，2）：

In [6]: x 

Out[6]: 
array([[1, 2, 3], 
       [2, 3, 4]]) 

In [7]: x[0,0] 

Out[7]: 1 

In [8]: x[1,2] 

Out[8]: 4 

ndarray对象具有许多有用的方法。 要获取可以在ndarray对象上调用的方法的列表，请在 IPython 提示符下键入array变量（在前面的示例中为x），然后按TAB。 这应该列出该对象可用的所有方法。 作为练习，尝试与其中一些玩耍。

数组索引和切片

NumPy 为数组提供了强大的索引功能。 NumPy 中的索引功能变得如此流行，以至于其中许多功能又重新添加到 Python 中。

在许多方面，为 NumPy 数组建立索引与为列表或元组建立索引非常相似。 存在一些差异，随着我们的进行，这些差异将变得显而易见。 首先，让我们创建一个尺寸为100 x 100的数组：

In [9]: x = np.random.random((100, 100)) 

简单的整数索引的工作方式是在一对方括号内键入索引，然后将其放置在数组变量旁边。 这是一种广泛使用的 Python 构造。 具有__getitem__方法的任何对象都将响应此类索引。 因此，要访问第 42 行和第 87 列中的元素，只需键入：

In [10]: y = x[42, 87] 

像列表和其他 Python 序列一样，也支持使用冒号为一系列值建立索引。 以下语句将打印x矩阵的第k行。

In [11]: print(x[k, :]) 

冒号可以被认为是所有元素的字符。 因此，前面的语句实际上意味着打印第k行的所有字符。 同样，可以使用x[:,k]访问列。 反转数组也类似于反转列表，例如x[::-1]。

数组的索引部分也称为数组的切片，它创建端口或整个数组的副本（我们将在后面的部分中介绍副本和视图） 。 在数组的上下文中，单词“切片”和“索引”通常可以互换使用。

下图显示了不同切片和索引技术的非常简洁的概述：

ndarray的内存布局

ndarray对象的一个​​特别有趣的属性是flags。 输入以下代码：

In [12]: x.flags 

它应该产生如下内容：

Out[12]: 
  C_CONTIGUOUS : True 
  F_CONTIGUOUS : False 
  OWNDATA : True 
  WRITEABLE : True 
  ALIGNED : True 
  UPDATEIFCOPY : False 

flags属性保存有关数组的内存布局的信息。 输出中的C_CONTIGUOUS字段指示该数组是否为 C 样式数组。 这意味着该数组的索引就像 C 数组一样完成。 在 2D 数组的情况下，这也称为行优先索引。 这意味着，当在数组中移动时，行索引将首先增加，然后列索引将增加。 在多维 C 样式数组的情况下，最后一个维度首先递增，然后是最后一个，但最后一个递增，依此类推。

同样，F_CONTIGUOUS属性指示该数组是否为 Fortran 样式的数组。 据说这样的数组具有列主索引（R，Julia 和 MATLAB 使用列主数组）。 这意味着，当在数组中移动时，第一个索引（沿着列）首先增加。

知道索引样式之间的差异很重要，尤其是对于大型数组，因为如果以正确的方式应用索引，则可以大大加快对数组的操作。 让我们通过练习来演示这一点。

声明一个数组，如下所示：

In [13]: c_array = np.random.rand(10000, 10000) 

这将产生一个名为c_array的变量，它是一个二维数组，其随机数为一亿。 （我们使用了 NumPy 中random子模块中的rand函数，我们将在下一部分中对其进行处理）。 接下来，从c_array创建一个 Fortran 样式的数组，如下所示：

In [14]: f_array = np.asfortranarray(c_array) 

您可以通过读取flags属性来分别检查c_array和f_array是否确实为 C 和 Fortran 风格。 接下来，我们定义以下两个函数：

In [15]: def sum_row(x):
         '''
         Given an array `x`, return the sum of its zeroth row.
         '''
         return np.sum(x[0, :])
In [16]: def sum_col(x):
         '''
         Given an array `x`, return the sum of its zeroth column.
         '''
         return np.sum(x[:, 0])

现在，让我们使用 IPython 的%timeit魔术函数在两个数组上测试这两个函数的性能：

注意

IPython 提供了一些魔术函数来帮助我们更好地理解代码。 有关更多详细信息，请参见这里。

In [17]: %timeit sum_row(c_array) 
10000 loops, best of 3: 21.2 µs per loop 

In [18]: %timeit sum_row(f_array) 
10000 loops, best of 3: 157 µs per loop 

In [19]: %timeit sum_col(c_array) 
10000 loops, best of 3: 162 µs per loop 

In [20]: %timeit sum_col(f_array) 
10000 loops, best of 3: 21.4 µs per loop 

如我们所见，对 C 数组的行求和比对其列求和要快得多。 这是因为，在 C 数组中，一行中的元素被放置在连续的内存位置中。 对于 Fortran 数组，情况恰好相反，其中列的元素布置在连续的内存位置中。

提示

请注意，确切的数字可能会因所使用的操作系统，RAM 和 Python 发行版而异，但是执行时间之间的相对顺序应保持不变。

这是一个重要的区别，它使您可以根据要执行的算法或运算的类型将数据适当地排列在一个数组中。 知道这种区别可以帮助您将代码加速几个数量级。

视图和副本

通过切片和索引访问数据主要有两种方法。 它们称为副本和视图：您可以直接从数组访问元素，也可以创建仅包含访问的元素的数组副本。 由于视图是原始数组的引用（在 Python 中，所有变量都是引用），因此修改视图也将修改原始数组。 对于副本，情况并非如此。

NumPy 其他例程中的may_share_memory函数可用于确定两个数组是彼此的副本还是彼此的视图。 尽管此方法在大多数情况下都能胜任，但由于使用启发式方法，因此并不总是可靠的。 它也可能返回错误的结果。 但是，出于介绍性目的，我们将其视为理所当然。

通常，切片数组会创建一个视图，对其进行索引会创建一个副本。 让我们通过一些代码片段研究这些差异。 首先，让我们创建一个随机的100x10数组。

In [21]: x = np.random.rand(100,10) 

现在，让我们提取数组的前五行，并将它们分配给变量y。

In [22]: y = x[:5, :] 

让我们看看y是否是x的视图。

In [23]: np.may_share_memory(x, y) 

Out[23]: True 

现在让我们修改数组y，看看它如何影响x。 将y的所有元素设置为零：

In [24]: y[:] = 0 

In [25]: print(x[:5, :]) 
Out[25]: [[ 0\.  0\.  0\.  0\.  0\.  0\.  0\.  0\.  0\.  0.] 
[ 0\.  0\.  0\.  0\.  0\.  0\.  0\.  0\.  0\.  0.] 
[ 0\.  0\.  0\.  0\.  0\.  0\.  0\.  0\.  0\.  0.] 
[ 0\.  0\.  0\.  0\.  0\.  0\.  0\.  0\.  0\.  0.] 
          [ 0\.  0\.  0\.  0\.  0\.  0\.  0\.  0\.  0\.  0.]] 

该代码段打印出五行零。 这是因为y只是一个视图，是对x的引用。

接下来，让我们创建一个副本以查看区别。 我们使用前面的方法使用随机函数来创建x数组，但是这次我们使用numpy.empty初始化y数组以首先创建一个空数组，然后将值从x复制到y。 因此，现在y不再是x的视图/参考； 它是一个独立的数组，但具有与x相同的值。 让我们再次使用may_share_memory函数来验证y是x的副本：

In [26]: x = np.random.rand(100,10) 

In [27]: y = np.empty([5, 10]) 

In [28]: y[:] = x[:5, :] 

In [29]: np.may_share_memory(x, y) 
Out[29]: False 

让我们更改y中的值，并检查x的值是否也发生变化：

In [30]: y[:] = 0 
In [31]: print(x[:5, :]) 

您应该看到前面的代码段在初始化x时打印了五行随机数，因此将y更改为0不会影响x。

创建数组

数组可以通过多种方式创建，例如从其他数据结构，通过读取磁盘上的文件或从 Web 创建。 就本章而言，其目的是使我们熟悉 NumPy 数组的核心特性，我们将使用列表或各种 NumPy 函数创建数组。

从列表创建数组

创建数组的最简单方法是使用array函数。 若要创建有效的数组对象，数组函数的参数必须至少满足以下条件之一：

• 它必须是有效的可迭代值或序列，可以嵌套
• 它必须具有返回有效的 numpy 数组的__array__方法

考虑以下代码段：

In [32]: x = np.array([1, 2, 3]) 

In [33]: y = np.array(['hello', 'world']) 

第一个条件对于 Python 列表和元组始终为true。 从列表或元组创建数组时，输入可能包含不同的（异构）数据类型。 但是，数组函数通常会将所有输入元素转换为数组所需的最合适的数据类型。 例如，如果列表同时包含浮点数和整数，则结果数组将为float类型。 如果它包含一个整数和一个布尔值，则结果数组将由整数组成。 作为练习，请尝试从包含任意数据类型的列表创建数组。

使用range函数是创建整数列表以及数组的最便捷方法之一：

In [34]: x = range(5) 

In [35]: y = np.array(x) 

NumPy 具有称为arange的便捷函数，该函数结合了range和array函数的功能。 前两行代码与此等效：

In [36]: x = np.arange(5) 

对于多维数组，只需嵌套输入列表，如下所示：

In [37]: x = np.array([[1, 2, 3],[4, 5, 6]]) 

In [38]: x.ndim 
Out[38]: 2 

In [39]: x.shape 
Out[39]: (2, 3) 

前面的示例仅说明如何从现有数组或一系列数字创建 NumPy 数组。 接下来，我们将讨论创建具有随机数的数组。

创建随机数组

NumPy 中的random模块提供了各种函数来创建任何数据类型的随机数组。 在整本书中，我们将非常频繁地使用此模块来演示 NumPy 中函数的工作。 random模块大致由以下函数组成：

• 创建随机数组
• 创建数组的随机排列
• 生成具有特定概率分布的数组

我们将在本书中详细介绍所有这些内容。 在本章中，我们将重点介绍random模块中的两个重要函数-rand和random。 这是一个简单的代码段，展示了这两个函数的使用：

In [40]: x = np.random.rand(2, 2, 2) 

In [41]: print(x.shape) 
Out[41]: (2, 2, 2) 

In [42]: shape_tuple = (2, 3, 4) 

In [43]: y = np.random.random(shape_tuple) 

In [44]: print(y.shape) 
Out[44]: (2, 3, 4) 

注意传递给两个函数的参数之间的细微差别。 随机函数接受元组作为参数，并创建维数等于元组长度的数组。 各个尺寸的长度等于元组的元素。 另一方面，rand函数采用任意数量的整数参数，并返回一个随机数组，使得其维数等于传递给该函数的整数参数的数量 ，并且各个维度的长度等于整数参数的值。 因此，前面的代码段中的x是三维数组（传递给函数的参数数量），并且x的三个维度中的每个维度的长度均为2（每个参数的值）。 rand是random的便捷函数。 这两个函数可以互换使用，只要传递的参数分别对两个函数均有效即可。

但是，这两个函数的主要缺点是-它们只能创建浮点数组。 如果我们想要一个随机整数数组，则必须将这些函数的输出转换为整数。 但是，这也是一个重大问题，因为 NumPy 的int函数将浮点数截断为最接近 0 的整数（这与floor函数等效）。 因此，将rand或random的输出强制转换为整数将始终返回零数组，因为这两个函数都返回[0, 1)区间内的浮点数。 可以使用randint函数解决此问题，如下所示：

In [45]: LOW, HIGH = 1, 11 

In [46]: SIZE = 10 

In [47]: x = np.random.randint(LOW, HIGH, size=SIZE) 

In [48]: print(x) 
Out[48]: [ 6  9 10  7  9  5  8  8  9  3] 

randint函数带有三个参数，其中两个是可选的。 第一个参数表示输出值的期望下限，第二个可选参数表示输出值的（专有）上限。 可选的size参数是一个元组，用于确定输出数组的形状。

还有许多其他函数，例如将随机数生成器植入随机子模块中。 有关详细信息，请参考这里。

其他数组

还有一些其他的数组创建函数，例如zeros()，ones()，eye()和其他一些函数（类似于 MATLAB 中的函数）可用于创建 NumPy 数组。 它们的使用非常简单。 数组也可以从文件或从 Web 填充。 我们将在下一章中处理文件 I/O。

数组的数据类型

数据类型是 NumPy 数组的另一个重要内在方面，它的内存布局和索引也是如此。 只需检查数组的dtype属性即可找到 NumPy 数组的数据类型。 尝试以下示例检查不同数组的数据类型：

In [49]: x = np.random.random((10,10)) 

In [50]: x.dtype 
Out[50]: dtype('float64') 
In [51]: x = np.array(range(10)) 

In [52]: x.dtype 
Out[52]: dtype('int32') 

In [53]: x = np.array(['hello', 'world']) 

In [54]: x.dtype 
Out [54]: dtype('S5') 

许多数组创建函数提供默认的数组数据类型。 例如，np.zeros和np.ones函数默认创建充满浮点数的数组。 但是也可以使它们创建其他数据类型的数组。 考虑以下示例，这些示例演示如何使用 dtype 参数创建任意数据类型的数组。

In [55]: x = np.ones((10, 10), dtype=np.int) 

In [56]: x.dtype 
Out[56]: dtype('int32') 

In [57]: x = np.zeros((10, 10), dtype='|S1') 

In [58]: x.dtype 
Out[58]: dtype('S1') 

有关 NumPy 支持的数据类型的完整列表，请参考这里。

总结

在本章中，我们介绍了 NumPy ndarray对象的一些基础知识。 我们研究了一些创建 NumPy 数组的基本方法。 我们还研究了数组的副本和视图之间的差异，以及它们如何影响使用索引和切片的情况。 我们看到了 NumPy 提供的内存布局之间的细微差别。 现在，我们已经配备了ndarray对象的基本词汇，并且可以开始使用 NumPy 的核心功能。 在下一章中，我们将探索ndarray的更多细节，并使用某些技巧和窍门（通用函数和形状操作）向您展示其中的一些技巧，以使您的 NumPy 脚本加速！

三、使用 NumPy 数组

NumPy 数组的优点在于您可以使用数组索引和切片来快速访问数据或执行计算，同时保持 C 数组的效率。 还支持许多数学运算。 在本章中，我们将深入研究使用 NumPy 数组。 在本章结束之后，您将对使用 NumPy 数组及其大部分功能感到满意。

这是本章将涉及的主题列表：

• NumPy 数组的基本操作和属性
• 通用函数（ufunc）和辅助函数
• 广播规则和形状操作
• 屏蔽 NumPy 数组

向量化运算

所有 NumPy 操作都是向量化的，您可以将操作应用于整个数组，而不是分别应用于每个元素。 与使用循环相比，这不仅整齐方便，而且还提高了计算性能。 在本节中，我们将体验 NumPy 向量化操作的强大功能。 在开始探索此主题之前，一个值得牢记的关键思想是始终考虑整个数组集而不是每个元素。 这将帮助您享受有关 NumPy 数组及其性能的学习。 让我们从标量和 NumPy 数组之间进行一些简单的计算开始：

In [1]: import numpy as np 
In [2]: x = np.array([1, 2, 3, 4]) 
In [3]: x + 1 
Out[3]: array([2, 3, 4, 5]) 

数组中的所有元素通过1同时添加。 这与 Python 或大多数其他编程语言有很大不同。 NumPy 数组中的元素都具有相同的dtype； 在前面的示例中，这是numpy.int（根据计算机的不同是 32 位或 64 位）; 因此，NumPy 可以节省在运行时检查每个元素的类型的时间，这通常是由 Python 完成的。 因此，只需应用以下算术运算：

In [4]: y = np.array([-1, 2, 3, 0]) 
In [5]: x * y 
Out[5]: array([-1,  4,  9,  0]) 

两个 NumPy 数组逐个元素相乘。 在前面的示例中，两个数组的形状相同，因此此处不应用广播（我们将在后面的部分中解释不同的形状，NumPy 数组操作和广播规则。）数组x中的第一个元素乘以数组y中的第一个元素，依此类推。 这里要注意的重要一点是，两个 NumPy 数组之间的算术运算不是矩阵乘法。 结果仍然返回相同形状的 NumPy 数组。 NumPy 中的矩阵乘法将使用numpy.dot()。 看一下这个例子：

In [6]: np.dot(x, y) 
Out[6]: 12 

NumPy 还支持两个数组之间的逻辑比较，并且比较也被向量化。 结果返回一个布尔值，并且 NumPy 数组指示两个数组中的哪个元素相等。 如果比较两个不同形状的数组，结果将仅返回一个False，这表明两个数组不同，并且实际上将比较每个元素：

In [7]: x == y 
Out[7]: array([False,  True,  True, False], dtype=bool) 

从前面的示例中，我们可以深入了解 NumPy 的元素操作，但是使用它们的好处是什么？ 我们怎么知道通过这些 NumPy 操作进行了优化？ 我们将使用上一章中介绍的 IPython 中的%timeit函数，向您展示 NumPy 操作和 Python for循环之间的区别：

In [8]: x = np.arange(10000) 
In [9]: %timeit x + 1 
100000 loops, best of 3: 12.6 µs per loop 
In [10]: y = range(10000) 
In [11]: %timeit [i + 1 for i in y] 
1000 loops, best of 3: 458 µs per loop 

x 和y这两个变量的长度相同，并且执行相同的工作，其中包括向数组中的所有元素添加值。 在 NumPy 操作的帮助下，性能比普通的 Python for循环要快得多（我们在这里使用列表推导来编写整洁的代码，这比普通的 Python for循环要快，但是与普通的 Python for循环相比，NumPy 的性能却更好）。 知道这个巨大的区别可以通过用 NumPy 操作替换循环来帮助您加速代码。

正如我们在前面的示例中提到的，性能的提高归因于 NumPy 数组中一致的dtype。 可以帮助您正确使用 NumPy 数组的技巧是在执行任何操作之前始终考虑dtype ，因为您很可能会在大多数编程语言中进行此操作。 下面的示例将为您展示使用相同操作的巨大不同结果，但这是基于不同的dtype数组：

In [12]: x = np.arange(1,9) 
In [13]: x.dtype 
Out[13]: dtype('int32') 
In [14]: x = x / 10.0 
In [15]: x 
Out[15]: array([ 0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8]) 
In [16]: x.dtype 
Out[16]: dtype('float64') 
In [17]: y = np.arange(1,9) 
In [18]: y /= 10.0 
In [19]: y 
Out[19]: array([0, 0, 0, 0, 0, 0, 0, 0]) 
In [20]: y.dtype 
Out[20]: dtype('int32') 

两个变量x和y完全相同：都是numpy.int32数组，范围从1到8（如果使用 64 位计算机，则可能会得到numpy.int64）并除以float 10.0。 但是，当x除以浮点数时，将使用dtype = numpy.float64创建一个新的 NumPy 数组。 这是一个全新的数组，但是具有相同的变量名x，因此x中的dtype进行了更改。 另一方面，y使用/=符号，该符号始终沿用y数组的dtype值。 因此，当它除以10.0时，不会创建新的数组; 仅更改y元素中的值，但dtype 仍为numpy.int32。 这就是x和y最终具有两个不同数组的原因。 请注意，从 1.10 版本开始，NumPy 不允许将浮点结果强制转换为整数。 因此，必须提高TypeError。

通用函数（ufunc）

NumPy 具有许多通用函数（所谓的ufunc），因此可以利用它们来发挥自己的优势，从而尽可能地减少循环以优化代码。 ufunc在数学，三角函数，汇总统计信息和比较运算方面有很好的覆盖范围。 有关详细的ufunc列表，请参考在线文档。

由于 NumPy 中有大量的ufunc，我们很难在一章中涵盖所有这些函数。 在本节中，我们仅旨在了解如何以及为何应使用 NumPy ufuncs。

基本ufunc入门

大多数ufunc都是一元或二进制的，这意味着它们只能接受一个或两个参数，然后逐元素地或在数学中应用它们。 这称为向量化运算或 NumPy 算术运算，我们已在前面的部分中进行了说明。 以下是一些常见的ufunc：

In [21]: x = np.arange(5,10) 
In [22]: np.square(x) 
Out[22]: array([25, 36, 49, 64, 81]) 

ufuncs 中广泛支持数学运算，其中一些基本与numpy.square()或numpy.log()基本相同，而另一些则是高级三角运算，例如numpy.arcsin()，numpy.rad2deg()等。 在这里np.mod()检索除法的余数：

In [23]: y = np.ones(5) * 10 
In [24]: np.mod(y, x) 
Out[24]: array([ 0.,  4.,  3.,  2.,  1.]) 

一些ufunc具有相似的名称，但它们的功能和行为却大不相同。 首先查看在线文档，以确保获得期望的结果。 这是numpy.minimum() 和numpy.min()的示例：

In [25]: np.minimum(x, 7) 
Out[25]: array([5, 6, 7, 7, 7]) 
In [26]: np.min(x) 
Out[26]: 5 

如您所见，numpy.minimum()比较两个数组并返回两个数组的最小值。 1 是数组值的形状，其值为 7，因此将其广播到[7, 7, 7, 7, 7]。 我们将在下一节中讨论 NumPy 广播规则。 numpy.min()，仅接受一个必需的参数，并返回数组中最小的元素。

使用更高级的函数

大多数ufunc具有可选参数，以在使用它们时提供更大的灵活性。 以下示例将使用numpy.median()。 这是在numpy.repeat()函数创建的二维数组上使用可选的axis 参数完成的，以重复x数组三次并将其分配给z变量：

In [27]: z = np.repeat(x, 3).reshape(5, 3) 
In [28]: z 
Out[28]: 
array([[5, 5, 5], 
       [6, 6, 6], 
       [7, 7, 7], 
       [8, 8, 8], 
       [9, 9, 9]]) 
In [29]: np.median(z) 
Out[29]: 7.0 
In [30]: np.median(z, axis = 0) 
Out[30]: array([ 7.,  7.,  7.]) 
In [31]: np.median(z, axis = 1) 
Out[31]: array([ 5.,  6.,  7.,  8.,  9.]) 

我们可以不使用axis参数就可以看到numpy.median()函数默认情况下会展平数组并返回中值元素，因此仅返回一个值。 使用axis自变量，如果将其应用于 0，则该操作将基于该列； 因此，我们获得了一个新的 NumPy 数组，其长度为3（z变量中总共有3列）。 虽然axis = 1，它基于行执行操作，所以我们有了一个包含五个元素的新数组。

ufuncs 不仅提供可选参数来调整操作，而且其中许多还具有一些内置方法，从而提供了更大的灵活性。 以下示例使用numpy.add()中的accumulate()累积对所有元素应用add()的结果：

In [32]: np.add.accumulate(x) 
Out[32]: array([ 5, 11, 18, 26, 35]) 

第二个示例将numpy.multiply()上的矩阵外部运算应用于来自两个输入数组的所有元素对。 在此示例中，两个数组来自x。 multiply()的外部产品的最终形状为5x5：

In [33]: np.multiply.outer(x, x) 
Out[33]: 
array([[25, 30, 35, 40, 45], 
       [30, 36, 42, 48, 54], 
       [35, 42, 49, 56, 63], 
       [40, 48, 56, 64, 72], 
       [45, 54, 63, 72, 81]]) 

如果您需要更高级的函数，则可以考虑构建自己的ufunc，这可能需要使用 Python-C API，或者您也可以使用 Numba 模块（向量化装饰器）来实现自定义的ufunc。 在本章中，我们的目标是了解 NumPy ufunc，因此我们将不介绍自定义的ufunc。 有关更多详细信息，请参阅 NumPy 的联机文档，名为编写自己的ufunc或 Numba 文档，创建 Numpy 通用函数。

广播和形状操作

NumPy 操作大部分是按元素进行的，这需要一个操作中的两个数组具有相同的形状。 但是，这并不意味着 NumPy 操作不能采用两个形状不同的数组（请参阅我们在标量中看到的第一个示例）。 NumPy 提供了在较大的数组上广播较小尺寸的数组的灵活性。 但是我们不能将数组广播成几乎任何形状。 它需要遵循某些约束； 我们将在本节中介绍它们。 要记住的一个关键思想是广播涉及在两个不同形状的数组上执行有意义的操作。 但是，不当广播可能会导致内存使用效率低下，从而减慢计算速度。

广播规则

广播的一般规则是确定两个数组是否与尺寸兼容。 需要满足两个条件：

• 两个数组的大小应相等
• 其中之一是 1

如果不满足上述条件，将引发ValueError异常，以指示数组具有不兼容的形状。 现在，我们将通过三个示例来研究广播规则的工作原理：

In [35]: x = np.array([[ 0, 0, 0], 
   ....:               [10,10,10], 
   ....:               [20,20,20]]) 
In [36]: y = np.array([1, 2, 3]) 
In [37]: x + y 
Out[37]: 
array([[ 1,  2,  3], 
       [11, 12, 13], 
       [21, 22, 23]]) 

让我们将前面的代码制作成图表，以帮助我们理解广播。 x变量的形状为(3, 3)，而y的形状仅为 3。但是在 NumPy 广播中，y的形状转换为1x3； 因此，该规则的第二个条件已得到满足。 通过重复将y广播到x的相同形状。 +操作可以按元素应用。

Numpy broadcasting on different shapes of arrays, where x(3,3) + y(3)

接下来，我们将向您展示广播两个数组的结果：

In [38]: x = np.array([[0], [10], [20]]) 
In [39]: x 
Out[39]: 
array([[ 0], 
       [10], 
       [20]]) 
In [40]: x + y 
Out[40]: 
array([[ 1,  2,  3], 
       [11, 12, 13], 
       [21, 22, 23]]) 

前面的示例向您展示x和y的广播方式。 x按列广播，而y按行广播，因为它们的形状在形状上均等于1。 满足第二个广播条件，并且新结果数组是3x3。

让我们看一下最后一个示例，其中两个数组不能满足广播规则的要求：

In [41]: x = np.array([[ 0, 0, 0], 
   ....:               [10,10,10], 
   ....:               [20,20,20]]) 
In [42]: y = np.arange(1,5) 
In [43]: x + y 
ValueError: operands could not be broadcast together with shapes (3,3) (4) 

在第三个示例中，由于x和y在行维度上具有不同的形状，并且它们都不等于1，因此无法执行广播。 因此，不能满足任何广播条件。 NumPy 抛出ValueError，告诉您形状不兼容。

重塑 NumPy 数组

了解广播规则之后，这里的另一个重要概念是重塑 NumPy 数组，尤其是在处理多维数组时。 通常只在一个维度上创建一个 NumPy 数组，然后将其重塑为多维，反之亦然。 这里的一个关键思想是，您可以更改数组的形状，但不应更改元素的数量。 例如，您无法将3xe数组整形为10x1数组。 整形前后，元素的总数（或ndarray内部组织中的所谓数据缓冲区）应保持一致。 或者，您可能需要调整大小，但这是另一回事了。 现在，让我们看一些形状操作：

In [44]: x = np.arange(24) 
In [45]: x.shape = 2, 3, -1 
In [46]: x 
Out[46]: 
array([[[ 0,  1,  2,  3], 
        [ 4,  5,  6,  7], 
        [ 8,  9, 10, 11]], 
       [[12, 13, 14, 15], 
        [16, 17, 18, 19], 
        [20, 21, 22, 23]]]) 

基本的重塑技术会更改numpy.shape属性。 在前面的示例中，我们有一个形状为(24,1)的数组，更改了shape属性后，我们获得了一个相同大小的数组，但是形状已更改为2x3x4组成。 注意， -1的形状是指转移数组的剩余形状尺寸。

In [47]: x = np.arange(1000000) 
In [48]: x.shape = 100, 100, 100 
In [49]: %timeit x.flatten() 
1000 loops, best of 3: 1.14 ms per loop 
In [50]: %timeit x.ravel() 
1000000 loops, best of 3: 330 ns per loop 

前面的示例是将100 x 100 x 100数组整形为一个尺寸; 在这里，我们应用numpy.flatten()和numpy.ravel()这两个函数来折叠数组，同时我们还比较了执行时间。 我们注意到numpy.flatten()和numpy.ravel()之间的速度差异很大，但是它们都比三层 Python 循环快得多。 两个函数在性能上的差异是np.flatten()从原始数组创建副本，而np.ravel()只是更改视图（如果您不记得副本和视图之间的区别，请回到第 2 章， “NumPy ndarray对象”）。

这个例子只是向您展示了 NumPy 提供了许多功能，其中一些可以产生相同的结果。 选择满足您目的的功能，同时为您提供优化的性能。

向量堆叠

重塑会更改一个数组的形状，但是如何通过大小相等的行向量构造二维或多维数组呢？ NumPy 为这种称为向量堆叠的解决方案提供了解决方案。 在这里，我们将通过三个示例使用三个不同的栈函数来实现基于不同维度的两个数组的组合：

In [51]: x = np.arange (0, 10, 2) 
In [52]: y = np.arange (0, -5, -1) 
In [53]: np.vstack([x, y]) 
Out[53]: 
array([[ 0,  2,  4,  6,  8], 
          [ 0, -1, -2, -3, -4]]) 

Numpy.vstack()通过垂直堆叠两个输入数组来构造新数组。 新数组是二维的：

In [54]: np.hstack([x, y]) 
Out[54]: array([ 0,  2,  4,  6,  8,  0, -1, -2, -3, -4]) 

numpy.hstack()水平合并两个数组时，新数组仍是一维的：

In [55]: np.dstack([x, y]) 
Out[55]: 
array([[[ 0,  0], 
        [ 2, -1], 
        [ 4, -2], 
        [ 6, -3], 
        [ 8, -4]]]) 

numpy.dstack()有点不同：它沿三维方向在深度方向上按顺序堆叠数组，因此新数组是三维的。

在下面的代码中，如果您使用numpy.resize()更改数组大小，则您正在放大数组，它将重复自身直到达到新大小； 否则，它将把数组截断为新的大小。 这里要注意的一点是ndarray也具有resize()操作，因此在此示例中，您还可以通过键入x.resize(8)来使用它来更改数组的大小； 但是，您会注意到放大部分填充了零，而不是重复数组本身。 另外，如果您已将数组分配给另一个变量，则无法使用ndarray.resize()。 Numpy.resize()创建一个具有指定形状的新数组，该数组的限制比ndarray.resize()少，并且是在需要时用于更改 NumPy 数组大小的更可取的操作：

In [56]: x = np.arange(3) 
In [57]: np.resize(x, (8,)) 
Out[57]: array([0, 1, 2, 0, 1, 2, 0, 1]) 

布尔掩码

在 NumPy 中，索引和切片非常方便且功能强大，但是使用布尔掩码，效果会更好！ 让我们首先创建一个布尔数组。 请注意，NumPy 中有一种特殊的数组，称为掩码数组。 在这里，我们不讨论它，但是我们还将解释如何使用 NumPy 数组扩展索引和切片：

In [58]: x = np.array([1,3,-1, 5, 7, -1]) 
In [59]: mask = (x < 0) 
In [60]: mask 
Out[60]: array([False, False,  True, False, False,  True], dtype=bool) 

从前面的示例中我们可以看到，通过应用<逻辑符号，我们将标量应用于 NumPy 数组，并将新数组命名为mask，它仍被向量化并返回与x形状相同的True/False 布尔值，表示x中的哪个元素符合标准：

In [61]: x [mask] = 0 
In [62]: x 
Out[62]: array([1, 3, 0, 5, 7, 0]) 

使用掩码，我们可以在不知道数组索引的情况下访问或替换数组中的任何元素值。 不用说，无需使用for循环即可完成此操作。

以下示例显示了如何对掩码数组求和，其中True代表 1，False 代表 0。我们创建了 50 个随机值，范围从0到1，其中 20 个大于0.5； 但是，对于随机数组，这是非常期望的：

In [69]: x = np.random.random(50) 
In [70]: (x > .5).sum() 
Out[70]: 20 

辅助函数

除了 Python 和其他在线文档中的help()和dir()函数之外，NumPy 还提供了一个辅助函数numpy.lookfor()来帮助您找到所需的正确函数。 参数是一个字符串，可以采用函数名称或任何与之相关的形式。 让我们尝试查找与resize相关的操作的更多信息，我们在前面的部分中进行了介绍：

In [71]: np.lookfor('resize') 
Search results for 'resize' 
--------------------------- 
numpy.ma.resize 
    Return a new masked array with the specified size and shape. 
numpy.chararray.resize 
    Change shape and size of array in-place. 
numpy.oldnumeric.ma.resize 
    The original array's total size can be any size. 
numpy.resize 
    Return a new array with the specified shape. 

总结

在本章中，我们介绍了 NumPy 及其ufunc的基本操作。 我们看了 NumPy 操作和 Python 循环之间的巨大差异。 我们还研究了广播的工作原理以及应避免的情况。 我们也试图理解掩蔽的概念。

使用 NumPy 数组的最好方法是尽可能地消除循环，并在 NumPy 中使用 ufuncs。 请记住广播规则，并谨慎使用它们。 将切片和索引与掩码一起使用可提高代码效率。 最重要的是，在使用时要玩得开心。

在接下来的几章中，我们将介绍 NumPy 的核心库，包括日期/时间和文件 I/O，以帮助您扩展 NumPy 的使用体验。

四、NumPy 核心和子模块

在上一章介绍了这么多 NumPy 函数之后，我希望您仍然记得 NumPy 的核心，即ndarray对象。 我们将完成ndarray的最后一个重要属性：步幅，它将为您提供完整的内存布局图。 另外，该向您展示 NumPy 数组不仅可以处理数字，还可以处理各种类型的数据； 我们将讨论记录数组和日期时间数组。 最后，我们将展示如何从文件中读取/写入 NumPy 数组，并开始使用 NumPy 进行一些实际的分析。

本章将涉及的主题是：

• NumPy 数组的核心：内存布局
• 结构数组（记录数组）
• NumPy 数组中的日期时间
• NumPy 数组中的文件 I/O

步幅

步幅是 NumPy 数组中的索引方案，它指示要跳转以查找下一个元素的字节数。 我们都知道 NumPy 的性能改进来自具有固定大小项的同构多维数组对象numpy.ndarray对象。 我们已经讨论了ndarray对象的shape（维度），数据类型和顺序（C 风格的行主要索引数组和 Fortran 风格的列主要数组）。现在让我们仔细看看步幅。

让我们首先创建一个 NumPy 数组并更改其形状以查看步幅的差异。

1. 创建一个 NumPy 数组，看一下步幅：

         In [1]: import numpy as np
         In [2]: x = np.arange(8, dtype = np.int8)
         In [3]: x
         Out[3]: array([0, 1, 2, 3, 4, 5, 6, 7])
         In [4]: x.strides
         Out[4]: (1,)
         In [5]: str(x.data)
         Out[5]: '\x00\x01\x02\x03\x04\x05\x06\x07'
   

   创建一维数组x，其数据类型为 NumPy 整数8，这意味着数组中的每个元素都是 8 位整数（每个 1 字节，总共 8 个字节）。 步幅表示遍历数组时在每个维度上步进的字节元组。 在上一个示例中，它是一个维度，因此我们将元组获得为(1, )。 每个元素与其前一个元素相距 1 个字节。 当我们打印出x.data时，我们可以获得指向数据开头的 Python 缓冲区对象，在示例中从x01到x07。

2. 更改形状并查看步幅变化：

         In [6]: x.shape = 2, 4 
         In [7]: x 
         Out[7]: 
         array([[0, 1, 2, 3], 
                [4, 5, 6, 7]], dtype=int8) 
         In [8]: x.strides 
         Out[8]: (4, 1) 
         In [9]: str(x.data) 
         Out[9]: '\x00\x01\x02\x03\x04\x05\x06\x07' 
         In [10]: x.shape = 1,4,2 
         In [11]: x.strides 
         Out[11]: (8, 2, 1) 
         In [12]: str(x.data) 
         Out[12]: '\x00\x01\x02\x03\x04\x05\x06\x07' 
   
   

   现在我们将x的尺寸更改为 2 乘以 4，然后再次检查步幅。 我们可以看到它变为(4, 1)，这意味着第一维中的元素相隔四个字节，并且数组需要跳转四个字节才能找到下一行，但是第二维中的元素仍相隔 1 个字节， 跳一个字节以查找下一列。 让我们再次打印出x.data，我们可以看到数据的内存布局保持不变，但是步幅改变了。 当我们将形状更改为三维时，会发生相同的行为：1 x 4 x 2数组。 （如果我们的数组是按照 Fortran 样式顺序构建的，会怎样？由于形状的变化，步幅将如何变化？尝试创建一个以列为主的数组，并执行相同的操作来检验这一点。）

3. 所以现在我们知道了什么是步幅，以及它与ndarray对象的关系，但是步幅如何改善我们的 NumPy 体验？ 让我们进行一些步幅操作以更好地理解这一点：两个数组的内容相同，但步幅却有所不同：

         In [13]: x = np.ones((10000,)) 
         In [14]: y = np.ones((10000 * 100, ))[::100] 
         In [15]: x.shape, y.shape 
         Out[15]: ((10000,), (10000,)) 
         In [16]: x == y 
         Out[16]: array([ True,  True,  True, ...,  True,  True,
         True], dtype=bool) 
   
   

4. 我们创建两个 NumPy 数组x和y并进行比较； 我们可以看到两个数组相等。 它们具有相同的形状，所有元素都是一个，但是实际上这两个数组在内存布局方面是不同的。 让我们简单地使用您在第 2 章， “NumPy ndarray对象”中了解的flags属性来检查两个数组的内存布局。

         In [17]: x.flags 
         Out[17]: C_CONTIGUOUS : True 
                  F_CONTIGUOUS : True 
                  OWNDATA : True 
                  WRITEABLE : True 
                  ALIGNED : True 
                  UPDATEIFCOPY : False 
   
         In [18]: y.flags 
         Out[18]: C_CONTIGUOUS : False 
                  F_CONTIGUOUS : False 
                  OWNDATA : False 
                  WRITEABLE : True 
                  ALIGNED : True 
                  UPDATEIFCOPY : False 
   
   

5. 我们可以看到x数组在 C 和 Fortran 顺序中都是连续的，而y则不是。 让我们检查一下步幅差异：

         In [19]: x.strides, y.strides 
         Out[19]: ((8,), (800,)) 
   
   

   数组x是连续创建的，因此在同一维中，每个元素相隔八个字节（numpy.ones的默认dtype是 64 位浮点数）； 但是，y是由每 100 个元素10000 * 100的子集创建的，因此内存布局中的索引架构不是连续的。

6. 尽管x和y具有相同的形状，但y中的每个元素彼此相距 800 个字节。 使用 NumPy 数组x和y时，您可能不会注意到索引的差异，但是内存布局确实会影响性能。 让我们使用 IPython 中的%timeit函数进行检查：

         In [18]: %timeit x.sum() 
         100000 loops, best of 3: 13.8 µs per loop 
         In [19]: %timeit y.sum() 
         10000 loops, best of 3: 25.9 µs per loop 
   
   

通常，在固定的高速缓存大小的情况下，当步幅大小变大时，命中率（在高速缓存中找到数据的内存访问量的百分比）相对较低，而未命中率（必须访问的内存访问量的百分比） 到内存）将会更高。 缓存命中时间和未命中时间构成了平均数据访问时间。 让我们尝试从缓存的角度再次来看我们的示例。 步幅较小的数组x快于y步幅较大的数组。 性能差异的原因是 CPU 将数据从主存储器分块地拉到其缓存中，步幅越小意味着需要的传输越少。 有关详细信息，请参见下图，其中红线表示 CPU 缓存的大小，蓝框表示包含数据的内存布局。

显然，如果同时需要x和y，则 100 个蓝色框的数据将减少x所需的缓存时间。

Cache and the x, y array in the memory layout

结构化数组

结构化数组或记录数组在执行计算时很有用，同时您可以将密切相关的数据保持在一起。 例如，当您处理事件数据并且每个事件都包含地理坐标和发生时间时，在计算最终结果时，您可以轻松地找到相关的地理位置和时间点以进行进一步的可视化。 NumPy 还提供了创建记录数组的强大功能，因为一个 NumPy 数组中存在多种数据类型。 但是，在 NumPy 中仍然需要遵守的一个原则是，每个字段中的数据类型（您可以将其视为记录中的列）需要是同质的。 以下是一些简单的示例，向您展示其工作方式：

In [20]: x = np.empty((2,), dtype = ('i4,f4,a10')) 
In [21]: x[:] = [(1,0.5, 'NumPy'), (10,-0.5, 'Essential')] 
In [22]: x 
Out[22]: 
array([(1, 0.5, 'NumPy'), (10, -0.5, 'Essential')], 
      dtype=[('f0', '<i4'), ('f1', '<f4'), ('f2', 'S10')]) 

在上一个示例中，我们使用numpy.empty()创建了一个一维记录数组，并指定了元素的数据类型-第一个元素是i4（32 位整数，其中i代表有符号整数， 4表示 4 个字节，例如np.int32），第二个元素是 32 位浮点数（f代表float也是 4 个字节），第三个元素是长度小于或等于 10 的字符串。我们按照指定的数据类型顺序将值分配给定义的数组。

您可以看到x的打印输出，该输出现在包含三种不同类型的记录，并且我们还在dtype中获得了默认字段名称：f0，f1和f2。 当然，您可以指定字段名称，如以下示例所示。

这里要注意的一件事是，我们使用了打印输出数据类型-i4和f4前面有一个<，而<代表字节顺序大端（指示内存地址增加顺序）：

In [23]: x[0] 
Out[23]: (1, 0.5, 'NumPy') 
In [24]: x['f2'] 
Out[24]: 
array(['NumPy', 'Essential'], dtype='|S10') 

检索数据的方式保持不变，我们使用索引来获取记录，但是此外，我们可以使用字段名称来获取某些字段的值，因此在上一个示例中，我们使用f2来获取字符串字段。 在下面的示例中，我们将创建名为y的x视图，并查看其如何与原始记录数组交互：

In [25]: y = x['f0'] 
In [26]: y 
Out[26]: array([ 1, 10]) 
In [27]: y[:] = y * 10 
In [28]: y 
Out[28]: array([ 10, 100]) 
In [29]: y[:] = y + 0.5 
In [30]: y 
Out[30]: array([ 10, 100]) 
In [31]: x 
Out[31]: 
array([(10, 0.5, 'NumPy'), (100, -0.5, 'Essential')], 
      dtype=[('f0', '<i4'), ('f1', '<f4'), ('f2', 'S10')]) 

这里y是x中字段f0的视图。 在记录数组中，NumPy 数组的特征仍然保留。 当您将标量 10 乘以时，它仍然适用于y的整个数组（广播规则），并且始终采用数据类型。 您可以在乘法之后看到，我们在y上添加了0.5，但是由于字段f0的数据类型是 32 位整数，结果仍然是[10, 100]。 另外，y是x中f0的视图，因此它们共享相同的内存块。 当我们在y中进行计算后打印出x时，我们发现x中的值也已更改。

在进一步介绍记录数组之前，让我们先整理一下如何定义记录数组。 最简单的方法如上一个示例所示，在该示例中，我们初始化 NumPy 数组，并使用字符串参数指定字段的数据类型。

NumPy 可以接受多种形式的字符串参数（有关详细信息，请参见这里； 最优选的可以选自以下之一：

数据类型	表示形式
b1	字节
i1，i2，i4，i8	1、2、4 和 8 个字节的带符号整数
u1，u2，u4，u8	1、2、4 和 8 个字节的无符号整数
f2，f4，f8	2、4 和 8 个字节的浮点数
c8，c16	8 和 16 个字节的复数
a<n>	长度为n的定长字符串

您也可以在字符串参数前面加上重复的数字或形状，以定义字段的维，但是在记录数组中，它仍仅被视为一个字段。 在下面的示例中，让我们尝试使用形状作为字符串参数的前缀：

In [32]: z = np.ones((2,), dtype = ('3i4, (2,3)f4')) 
In [32]: z 
Out[32]: 
array([([1, 1, 1], [[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]), 
       ([1, 1, 1], [[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]])], 
      dtype=[('f0', '<i4', (3,)), ('f1', '<f4', (2, 3))]) 

在前面的示例中，字段f0是尺寸为3的一维数组，f1是形状为(2, 3)的二维数组。 现在，我们很清楚记录数组的结构以及如何定义它。 您可能想知道默认字段名称是否可以更改为对您的分析有意义的名称？ 当然可以！ 这是这样的：

In [33]: x.dtype.names 
Out[33]: ('f0', 'f1', 'f2') 
In [34]: x.dtype.names = ('id', 'value', 'note') 
In [35]: x 
Out[35]: 
array([(10, 0.5, 'NumPy'), (100, -0.5, 'Essential')], 
      dtype=[('id', '<i4'), ('value', '<f4'), ('note', 'S10')]) 

通过将新的字段名称分配回dtype对象中的名称属性，我们可以获得自定义的字段名称。 或者，您可以通过使用带有元组的列表或字典来初始化记录数组时执行此操作。 在以下示例中，我们将使用列表和字典来创建两个具有自定义字段名称的相同记录数组：

In [36]: list_ex = np.zeros((2,), dtype = [('id', 'i4'), ('value', 'f4', (2,))]) 
In [37]: list_ex 
Out[37]: 
array([(0, [0.0, 0.0]), (0, [0.0, 0.0])], 
      dtype=[('id', '<i4'), ('value', '<f4', (2,))]) 
In [38]: dict_ex = np.zeros((2,), dtype = {'names':['id', 'value'], 'formats':['i4', '2f4']}) 
In [39]: dict_ex 
Out[39]: 
array([(0, [0.0, 0.0]), (0, [0.0, 0.0])], 
      dtype=[('id', '<i4'), ('value', '<f4', (2,))]) 

在列表示例中，我们为每个字段创建一个元组（字段名称，数据类型和形状）。 shape参数是可选的； 您也可以使用数据类型参数指定形状。 使用字典定义字段时，有两个必需的键（names和formats），每个键都有大小相等的值列表。

在继续下一节之前，我们将向您展示如何一次访问记录数组中的多个字段。 以下示例仍使用我们在自定义字段开头创建的数组x：id，value和note：

In [40]: x[['id', 'note']] 
Out[40]: 
array([(10, 'NumPy'), (100, 'Essential')], 
      dtype=[('id', '<i4'), ('note', 'S10')]) 

您可能会发现此示例过于简单； 如果是这样，您可以尝试使用来自维基百科的数据从包含国家名称，人口和排名的真实示例创建一个新的记录数组。 这样会更有趣！

NumPy 中的日期和时间

当您进行时间序列分析时，日期和时间很重要，从简单的事情（如在博物馆中累积每天的访客）到复杂的事情（如对犯罪预测的趋势回归）。 从 NumPy 1.7 开始，NumPy 核心支持日期时间类型（尽管它仍处于试验阶段，并且可能会发生变化）。 为了区别于 Python 中的datetime对象，数据类型称为datetime64。

本节将介绍numpy.datetime64的创建，时间增量算法以及单位与 Python 之间的转换datetime。 让我们使用 ISO 字符串创建一个numpy.datetime64对象：

In [41]: x = np.datetime64('2015-04-01') 
In [42]: y = np.datetime64('2015-04') 
In [43]: x.dtype, y.dtype 
Out[43]: (dtype('<M8[D]'), dtype('<M8[M]')) 

x和y都是numpy.datetime64对象，并由 ISO 8601 字符串构造（通用日期格式-有关详细信息，请参见这里。 但是x的输入字符串包含天单位，而y的字符串则没有。 创建 NumPy datetime64时，它将自动从输入字符串的形式中进行选择，因此当我们为x和y都打印出dtype时，我们可以看到x的单位为D。 ]代表几天，而y和单位M代表几个月。 <也是字节序，这里是大端，M8是datetime64的缩写（从np.int64实现）。 numpy.datetime64支持的默认日期单位是年（Y），月（M），周（W）和天（D），而时间单位是小时（h）。 ，分钟（m），秒（s）和毫秒（ms）。

当然，我们在创建数组时可以指定单位，也可以使用numpy.arange()方法创建数组的序列。 请参阅以下示例：

In [44]: y = np.datetime64('2015-04', 'D') 
In [45]: y, y.dtype 
Out[45]: (numpy.datetime64('2015-04-01'), dtype('<M8[D]')) 
In [46]: x = np.arange('2015-01', '2015-04', dtype = 'datetime64[M]') 
In [47]: x 
Out[47]: array(['2015-01', '2015-02', '2015-03'], dtype='datetime64[M]') 

但是，当 ISO 字符串仅包含日期单位时，不允许指定时间单位。 将触发TypeError，因为日期单位和时间单位之间的转换需要选择时区和给定日期的特定时间：

In [48]: y = np.datetime64('2015-04-01', 's') 
TypeError: Cannot parse "2015-04-01" as unit 's' using casting rule 'same_kind' 

接下来，我们将对两个numpy.datetime64数组进行减法运算，您将看到只要两个数组之间的日期/时间单位是可转换的，广播规则仍然有效。 我们使用先前创建的相同数组x，并为以下示例创建一个新的y：

In [49]: x 
Out[49]: array(['2015-01', '2015-02', '2015-03'], dtype='datetime64[M]') 
In [50]: y = np.datetime64('2015-01-01') 
In [51]: x - y 
Out[51]: array([ 0, 31, 59], dtype='timedelta64[D]') 

有趣的是，x减去y的结果数组是[0, 31, 59]，而不是日期，并且dtype更改为timedelta64[D]。 由于 NumPy 的核心没有物理量系统，因此创建了timedelta64数据类型以补充datetime64。 在上一个示例中，[0, 31, 59]是x中每个元素中从2015-01-01开始的单位，单位是天（D）。 您也可以在datetime64和timedelta64之间进行算术运算，如以下示例所示：

In [52]: np.datetime64('2015') + np.timedelta64(12, 'M') 
Out[52]: numpy.datetime64('2016-01') 
In [53]: np.timedelta64(1, 'W') / np.timedelta64(1, 'D') 
Out[53]: 7.0 

在本节的最后部分，我们将讨论numpy.datetime64和 Python datetime之间的转换。 尽管datetime64对象从 NumPy 数组继承了许多特征，但是使用 Python datetime对象（例如date和year属性，isoformat等）仍然有一些好处，反之亦然 。 例如，您可能具有datetime对象的列表，并且可能希望将其转换为用于算术或其他 NumPy 函数的numpy.datetime64。 在下面的示例中，我们将以两种方式将现有的datetime64数组x转换为 Python datetime列表：

In [54]: x 
Out[54]: array(['2015-01', '2015-02', '2015-03'], dtype='datetime64[M]') 
In [55]: x.tolist() 
Out[55]: 
[datetime.date(2015, 1, 1), 
 datetime.date(2015, 2, 1), 
 datetime.date(2015, 3, 1)] 
In [56]: [element.item() for element in x] 
Out[56]: 
[datetime.date(2015, 1, 1), 
 datetime.date(2015, 2, 1), 
 datetime.date(2015, 3, 1)] 

我们可以看到带有for循环的numpy.datetime64.tolist()和numpy.datetime64.item()可以实现相同的目标，即将数组转换为 Python datetime对象的列表。 但是不用说，我们都知道哪种方法更适合进行转换（如果您不知道答案，请快速浏览第 3 章，“使用 Numpy 数组”。）另一方面，如果您已经有了 Python datetime的列表，并想将其转换为 NumPy datetime64数组，则只需使用numpy.array()函数。

NumPy 文件 I/O

现在我们可以执行 NumPy 数组计算和操作，并且知道如何构造记录数组，现在是时候进行一些实际的分析了，方法是将文件读入 NumPy 数组并将结果数组输出到文件中以进行进一步的分析。

我们应该谈论先读取文件然后导出文件。 但是现在，我们将逆转此过程，先创建一个记录数组，然后将其输出到 CSV 文件。 我们将导出的 CSV 文件读入 NumPy 记录数组，并将其与原始记录数组进行比较。 我们将要创建的样本数组将包含一个带有连续整数的id字段，一个包含随机浮点数的value字段和一个带有numpy.datetime64['D']的date字段。 本练习将使用您从前面的章节中获得的所有知识。 让我们开始创建记录数组：

In [57]: id = np.arange(1000) 
In [58]: value = np.random.random(1000) 
In [59]: day = np.random.random_integers(0, 365, 1000) * np.timedelta64(1,'D') 
In [60]: date = np.datetime64('2014-01-01') + day 
In [61]: rec_array = np.core.records.fromarrays([id, value, date], names='id, value, date', formats='i4, f4, a10') 
In [62]: rec_array[:5] 
Out[62]: 
rec.array([(0, 0.07019801437854767, '2014-07-10'), 
       (1, 0.4863224923610687, '2014-12-03'), 
       (2, 0.9525277614593506, '2014-03-11'), 
       (3, 0.39706873893737793, '2014-01-02'), 
       (4, 0.8536589741706848, '2014-09-14')], 
      dtype=[('id', '<i4'), ('value', '<f4'), ('date', 'S10')]) 

我们首先创建代表所需字段的三个 NumPy 数组：id，value和date。 创建date字段时，我们将numpy.datetime64与大小为1000的随机 NumPy 数组结合使用，以模拟从2014-01-01到2014-12-31范围内的随机日期（365 天）。

然后我们使用numpy.core.records.fromarrays()函数将三个数组合并为一个记录数组，并分配names（字段名称）和formats（数据类型）。 这里要注意的一件事是，记录数组不支持numpy.datetime64对象，因此我们将其作为长度为 10 的日期/时间字符串存储在数组中。

如果您使用的是 Python 3，则会在记录数组（例如b'2014-09-25'）的日期/时间字符串的前面找到前缀b。 b在这里代表“字节字面值”，这意味着它仅包含 ASCII 字符（Python 3 中的所有字符串类型均为 Unicode，这是 Python 2 和 3 之间的一大变化）。 因此，在 Python 3 中，将对象（datetime64）转换为字符串将添加前缀以区分普通字符串类型。 但是，这不会影响下一步将记录数组导出到 CSV 文件中的操作：

In [63]: np.savetxt('./record.csv', rec_array, fmt='%i,%.4f,%s') 

我们使用numpy.savetxt()函数来处理导出，并使用fmt参数将第一个参数指定为导出文件的位置，数组名称和格式。 我们有三个具有三种不同数据类型的字段，我们想在 CSV 文件的每个字段之间添加,。 如果您更喜欢其他定界符，请在fmt参数中替换逗号。 我们还消除了value字段中的冗余数字，因此我们使用%.4f仅在文件的小数点后指定四位数字。 现在，您可以转到我们在第一个参数中指定的文件位置，以检查 CSV 文件。 在电子表格软件程序中将其打开，您将看到以下内容：

接下来，我们将 CSV 文件读取到记录数组中，并使用value字段生成一个名为mask的掩码字段，该掩码字段表示一个大于或等于 0.75 的值。 然后，我们将新的mask字段追加到记录数组。 让我们先阅读 CSV 文件：

In [64]: read_array = np.genfromtxt('./record.csv', dtype='i4,f4,a10', delimiter=',', skip_header=0) 
In [65]: read_array[:5] 
Out[65]: 
array([(0, 0.07020000368356705, '2014-07-10'), 
       (1, 0.486299991607666, '2014-12-03'), 
       (2, 0.9524999856948853, '2014-03-11'), 
       (3, 0.3971000015735626, '2014-01-02'), 
       (4, 0.8536999821662903, '2014-09-14')], 
      dtype=[('f0', '<i4'), ('f1', '<f4'), ('f2', 'S10')]) 

我们使用numpy.genfromtxt()将文件读入 NumPy 记录数组。 第一个参数仍然是我们要访问的文件位置，dtype参数是可选的。 如果未指定，NumPy 将使用各列的内容分别确定dtype参数。 由于我们清楚地了解数据，因此建议您在每次读取文件时指定一次。

delimiter参数也是可选的，默认情况下，任何连续的空格都用作分隔符。 但是，对于 CSV 文件，我们使用了,。 我们在方法中使用的最后一个可选参数是skip_header。 尽管我们没有在文件中的记录顶部添加字段名称，但是 NumPy 提供了跳过文件开头多行的功能。

除了skip_header之外，numpy.genfromtext()函数还支持 22 种以上的操作参数以微调数组，例如定义缺失值和填充值。 有关更多详细信息，请参阅这里。

现在将数据读到记录数组中，您将发现第二个字段是小数点后四位数以上，这是我们在导出 CSV 时指定的。 这样做的原因是因为我们在读取时使用f4作为其数据类型。NumPy 将填充空数字，但有效的四位数字与文件中的数字相同。 您可能还会注意到我们丢失了字段名，因此让我们指定它：

In [66]: read_array.dtype.names = ('id', 'value', 'date') 

本练习的最后一部分是创建一个基于value字段的值大于或等于0.75的掩码变量。 我们将新的掩码数组作为新列附加到read_array中：

In [68]: mask = read_array['value'] >= 0.75 
In [69]: from numpy.lib.recfunctions import append_fields 
In [70]: read_array = append_fields(read_array, 'mask', data=mask, dtypes='i1') 
In [71]: read_array[:5] 
Out[71]: 
masked_array(data = [(0, 0.07020000368356705, '2014-07-10', 0) 
 (1, 0.486299991607666, '2014-12-03', 0)

 (2, 0.9524999856948853, '2014-03-11', 1) 
 (3, 0.3971000015735626, '2014-01-02', 0) 
dtype = [('id', '<i4'), ('value', '<f4'), ('date', 'S10'), ('mask','i1')]) 

仅当直接导入numpy.lib.recfunctions且模块中具有append_field()函数时，才能访问它。 追加一个记录数组就像追加一个 NumPy 数组一样简单：第一个参数是基本数组；第二个参数是基本数组。 第二个参数是新字段名称mask以及与之关联的数据; 最后一个参数是数据类型。 由于掩码是布尔数组，因此 NumPy 会自动将掩码应用于记录数组，但是我们仍然可以看到在read_array中添加了一个新字段，掩码的值反映了阈值（>= 0.75） value字段。 这只是向您展示如何将 NumPy 数组与数据文件连接的开始。 现在是时候对您的数据进行一些真实的分析了！

总结

在本章中，我们介绍了ndarray对象的最后一个重要组成部分：步幅。 当您使用不同的方式初始化 NumPy 数组时，我们看到了内存布局和性能上的巨大差异。 我们还了解了记录数组（结构化数组）以及如何在 NumPy 中操纵日期/时间。 最重要的是，我们看到了如何使用 NumPy 读写数据。

NumPy 的强大功能不仅在于其性能或功能，还在于它使分析变得如此容易。 尽可能多地将 NumPy 与您的数据一起使用！

接下来，我们将研究使用 NumPy 进行线性代数和矩阵计算。

五、NumPy 中的线性代数

NumPy 专为数值计算而设计； 在引擎盖下，它仍然是功能强大的ndarray对象，但同时 NumPy 提供了不同类型的对象来解决数学问题。 在本章中，我们将介绍矩阵对象和多项式对象，以帮助您使用非 ndarray 方法解决问题。 同样，NumPy 提供了许多标准的数学算法并支持多维数据。 虽然矩阵无法执行三维数据，但更可取的是使用ndarray对象以及线性代数和多项式的 NumPy 函数（更广泛的 SciPy 库是线性代数的另一个不错的选择，但是 NumPy 是我们关注的重点） 书）。 现在让我们使用 NumPy 进行一些数学运算！

本章将涉及的主题是：

• 矩阵和向量运算
• 分解
• 多项式运算
• 回归和曲线拟合

矩阵类型

对于线性代数，使用矩阵可能更直接。 NumPy 中的矩阵对象继承了ndarray的所有属性和方法，但严格来说是二维的，而ndarray可以是多维的。 使用 NumPy 矩阵的众所周知的优点是它们提供矩阵乘法作为*表示法； 例如，如果x和y是矩阵，则x * y是它们的矩阵乘积。 但是，从 Python 3.5 / NumPy 1.10 开始，新的运算符“

但是，从 Python 3.5 / NumPy 1.10 开始，新的运算符@支持本机矩阵乘法。 因此这是使用ndarray的另一个很好的理由。

但是，矩阵对象仍提供方便的转换，例如逆和共轭转置，而ndarray不提供。 让我们从创建 NumPy 矩阵开始：

In [1]: import numpy as np 
In [2]: ndArray = np.arange(9).reshape(3,3) 
In [3]: x = np.matrix(ndArray) 
In [4]: y = np.mat(np.identity(3)) 
In [5]: x 
Out[5]: 
matrix([[0, 1, 2], 
        [3, 4, 5], 
        [6, 7, 8]]) 
In [6]: y 
Out[6]: 
matrix([[1., 0., 0.], 
        [0., 1., 0.], 
        [0., 0., 1.]]) 

有两种创建或转换为 NumPy 矩阵对象的方法，更优选的方法是使用numpy.mat()或numpy.matrix()。 两种方法都创建矩阵，但是numpy.matrix()创建一个副本，而numpy.mat()仅更改视图; 等效于numpy.matrix(data, copy = False)。 在前面的示例中，我们创建了两个矩阵，两个矩阵均来自ndarray对象（np.identity(3)返回3 x 3的单位数组）。 当然，您可以使用字符串或列表来创建矩阵，例如：np.matrix('0 1 2; 3 4 5; 6 7 8')和np.matrix([[0,1,2],[3,4,5],[6,7,8]])将创建与x相同的矩阵。 在以下示例中，我们将执行一些基本的矩阵运算：

In [7]: x + y 
Out[7]: 
matrix([[ 1.,  1.,  2.], 
        [ 3.,  5.,  5.], 
        [ 6.,  7.,  9.]]) 
In [8]: x * x 
Out[8]: 
matrix([[ 15,  18,  21], 
        [ 42,  54,  66], 
        [ 69,  90, 111]]) 
In [9]: np.dot(ndArray, ndArray) 
Out[9]: 
array([[ 15,  18,  21], 
       [ 42,  54,  66], 
       [ 69,  90, 111]]) 
In [10]: x**3 
Out[10]: 
matrix([[ 180,  234,  288], 
        [ 558,  720,  882], 
        [ 936, 1206, 1476]])
In [11]: z = np.matrix(np.random.random_integers(1, 50, 9).reshape(3,3)) 
In [12]: z 
Out[12]: 
matrix([[32, 21, 28], 
        [ 2, 24, 22], 
        [32, 20, 22]]) 
In [13]: z.I 
Out[13]: 
matrix( [[-0.0237 -0.0264  0.0566] 
         [-0.178   0.0518  0.1748] 
         [ 0.1963 -0.0086 -0.1958]]) 

In [14]: z.H 
Out[14]: 
matrix([[32  2 32] 
        [21 24 20] 
        [28 22 22]]) 

您可以从前面的示例中看到，当我们使用*表示法时，当您为ndarray使用numpy.dot()时，它会应用矩阵乘法（我们将在下一节中讨论）。 同样，**功率符号以矩阵方式完成。 我们还根据随机函数创建了一个矩阵z，以显示该矩阵何时可逆（非奇异）。 您可以使用numpy.matrix.I获得逆矩阵。 我们还可以使用numpy.matrix.H进行共轭（Hermitian）转置。

现在我们知道了如何创建矩阵对象并执行一些基本操作，是时候进行一些练习了。 让我们尝试求解一个简单的线性方程。 假设我们有一个线性方程A x = b，我们想知道x的值。 可能的解决方案如下：

A-1A x = A-1 b 
I x = A-1 b 
x = A-1 b 

我们通过将A和b的逆数相乘来获得x，所以我们用numpy.matrix来做到这一点：

In [15]: A = np.mat('3 1 4; 1 5 9; 2 6 5') 
In [16]: b = np.mat([[1],[2],[3]]) 
In [17]: x = A.I * b 
In [18]: x 
Out[18]: 
matrix([[ 0.2667], 
        [ 0.4667], 
        [-0.0667]]) 
In [21]: np.allclose(A * x, b) 
Out[21]: True 

我们获得了x，并使用numpy.allclose()在公差范围内比较了 LHS 和 RHS。 默认的绝对公差为1e-8。 结果返回True，这意味着 LHS 和 RHS 在公差范围内相等，这验证了我们的解决方案。 尽管numpy.matrix()采用普通矩阵形式，但在大多数情况下ndarray足以满足您进行线性代数的需要。 现在我们将简单比较ndarray和matrix的性能：

In [20]: x = np.arange(25000000).reshape(5000,5000) 

In [21]: y = np.mat(x) 

In [22]: %timeit x.T 
10000000 loops, best of 3: 176 ns per loop 

In [23]: %timeit y.T 
1000000 loops, best of 3: 1.36 µs per loop 

此示例显示了转置时ndarray和matrix之间的巨大性能差异。 x和y都具有5,000 x 5,000元素，但是x是二维ndarray，而y将其转换为相同的形状matrix。 即使计算已通过 NumPy 优化，NumPy 矩阵也将始终以矩阵方式进行运算。

虽然默认情况下ndarray会反转尺寸而不是对轴进行置换（矩阵始终对轴进行置换），但这是ndarray中完成的一项巨大的性能改进技巧。 因此，考虑到线性代数的性能，ndarray 特别适用于大型数据集。 仅在必要时使用matrix。 在继续下一节之前，让我们看一下另外两个matrix对象属性，这些属性将matrix转换​​为基本的ndarray：

In [24]: A.A 
Out[24]: 
array([[3, 1, 4], 
       [1, 5, 9], 
       [2, 6, 5]]) 
In [25]: A.A1 
Out[25]: array([3, 1, 4, 1, 5, 9, 2, 6, 5]) 

前面的示例使用了我们在线性方程实践中创建的矩阵A。 numpy.matrix.A返回基本的ndarray，numpy.matrix.A1返回一维的ndarray。

NumPy 中的线性代数

在进入 NumPy 的线性代数类之前，我们将在本节的开头介绍五种向量积。 让我们从numpy.dot()产品开始逐一回顾它们：

In [26]: x = np.array([[1, 2], [3, 4]]) 
In [27]: y = np.array([[10, 20], [30, 40]]) 
In [28]: np.dot(x, y) 
Out[28]: 
array([[ 70, 100], 
       [150, 220]]) 

numpy.dot()函数执行矩阵乘法，其详细计算如下所示：

numpy.vdot()处理多维数组的方式与numpy.dot()不同。 它不执行矩阵乘积，而是首先将输入参数展平到一维向量：

In [29]: np.vdot(x, y) 
Out[29]: 300 

numpy.vdot()的详细计算如下：

numpy.outer()函数是两个向量的外积。 如果输入数组不是一维的，则它将变平。 假设扁平化的输入向量A的形状为(M, )，扁平化的输入向量B的形状为(N, )。 那么结果形状将是(M, N)：

In [100]: np.outer(x,y) 
Out[100]: 
array([[ 10,  20,  30,  40], 
       [ 20,  40,  60,  80], 
       [ 30,  60,  90, 120], 
       [ 40,  80, 120, 160]]) 

numpy.outer()的详细计算如下：

最后一个是numpy.cross()乘积，它是三维空间中两个向量的二进制运算（并且仅适用于向量），其结果是垂直于两个输入数据的向量（a，b）。 如果您不熟悉外部产品，请参考这里。 以下示例显示a和b是向量数组，以及（a，b）和（b，a）的叉积：

In [31]: a = np.array([1,0,0]) 
In [32]: b = np.array([0,1,0]) 
In [33]: np.cross(a,b) 
Out[33]: array([0, 0, 1]) 
In [34]: np.cross(b,a) 
Out[34]: array([ 0,  0, -1]) 

下图显示了详细的计算，并且两个向量a与b的叉积由a x b表示：

NumPy 为标准向量例程提供了前面的功能。 现在，我们将讨论本节的关键主题：用于线性代数的numpy.linalg子模块。 将 NumPy ndarray与numpy.linalg结合使用会比numpy.matrix()更好。

注意

如果您希望将 scipy 作为程序的依赖项，则scipy.linalg具有比numpy.linalg更高级的功能，例如矩阵中的三角函数以及更多的分解选择。 但是，NumPy 包含所有基本操作。

在以下示例中，我们将介绍numpy.linalg的其余基本操作，并使用它们来求解矩阵部分中的线性方程：

In [35]: x = np.array([[4,8],[7,9]]) 
In [36]: np.linalg.det(x) 
Out[36]: -20.000000000000007 

前面的示例计算方阵的行列式。 当然，我们可以使用numpy.linalg.inv()来计算数组的逆数，就像我们使用numpy.matrix.I一样：

In [37]: np.linalg.inv(x) 
Out[37]: 
array([[-0.45,  0.4 ], 
       [ 0.35, -0.2 ]]) 
In [38]: np.mat(x).I 
Out[38]: 
matrix([[-0.45,  0.4 ], 
        [ 0.35, -0.2 ]]) 

从前面的示例中，我们可以看到numpy.linalg.inv()提供与numpy.matrix.I相同的结果。 唯一的区别是numpy.linalg返回ndarray。 接下来，我们将再次回到线性方程式A x = b，以了解如何使用numpy.linalg.solve()获得与使用矩阵对象相同的结果：

In [39]: x = np.linalg.solve(A,b) 
In [40]: x 
Out[40]: 
matrix([[ 0.2667], 
        [ 0.4667], 
        [-0.0667]]) 

numpy.linalg.solve(A,b)计算x的解，其中第一个输入参数（A）代表系数数组，第二个参数（b）代表坐标或因变量值。 numpy.linalg.solve()函数支持输入数据类型。 在示例中，我们使用矩阵作为输入，因此输出还返回一个矩阵x。 我们也可以使用ndarray 作为输入。

使用 NumPy 进行线性代数运算时，最好仅使用一种数据类型，即ndarray或matrix。 不建议在计算中使用混合类型。 原因之一是减少了不同数据类型之间的转换。 另一个原因是要避免两种类型的计算中的意外错误。 由于ndarray对数据尺寸的限制较少，并且可以执行所有类似矩阵的运算，因此与matrix相比，ndarray与numpy.linalg结合使用是首选的。

分解

numpy.linalg提供了分解，在本节中，我们将介绍两种最常用的分解：奇异值分解（svd）和 QR 因式分解。 让我们首先计算特征值和特征向量。 在我们开始之前，如果您不熟悉特征值和特征向量，可以在这里进行检查。 开始吧：

In [41]: x = np.random.randint(0, 10, 9).reshape(3,3) 
In [42]: x 
Out[42]: 
array([[ 1,  5,  0] 
       [ 7,  4,  0] 
       [ 2,  9,  8]]) 
In [42]: w, v = np.linalg.eig(x) 
In [43]: w 
Out[43]: array([ 8.,  8.6033,  -3.6033]) 
In [44]: v 
Out[44]: 
array([[ 0.,  0.0384,  0.6834] 
       [ 0.,  0.0583, -0.6292] 
       [ 1.,  0.9976,  0.3702]] 
) 

在前面的示例中，首先我们使用numpy.random.randint ()创建了一个3 x 3的ndarray，然后使用np.linalg.eig()计算了特征值和特征向量。 该函数返回两个元组：第一个元组是特征值，每个元组根据其多重性重复;第二个元组是规范化的特征向量，其中v[: , i]列是与特征值w[i]相对应的特征向量。 在此示例中，我们将元组解压缩为w和v。 如果输入ndarray是复数值，则计算出的特征向量也将是复数类型，如下面的示例所示：

In [45]: y = np.array([[1, 2j],[-3j, 4]]) 
In [46]: np.linalg.eig(y) 
Out[46]: 
(array([ -0.3723+0.j,  5.3723+0.j]), 
 array([[0.8246+0.j    ,  0.0000+0.416j     ], 
        [-0.0000+0.5658j,  0.9094+0.j    ]])) 

但是，如果输入ndarray是实数，则计算出的特征值也将是实数； 因此，在计算时，我们应注意舍入错误，如以下示例所示：

In [47]: z = np.array([[1 + 1e-10, -1e-10],[1e-10, 1 - 1e-10]]) 
In [48]: np.linalg.eig(z) 
Out[48]: 
(array([ 1.,  1.]), array([[0.70710678,  0.707106], 
        [0.70710678,  0.70710757]])) 

ndarrayz是实型（numpy.float64），因此在计算特征值时会自动四舍五入。 从理论上讲，特征值应为1 ± 1e-10，但从第一个np.linalg.eig()可以看出特征值都向上舍入为1。

svd可以认为是特征值的扩展。 我们可以使用numpy.linalg.svd()分解M x N数组，所以让我们从一个简单的例子开始：

In [51]: np.set_printoptions(precision = 4) 
In [52]: A = np.array([3,1,4,1,5,9,2,6,5]).reshape(3,3) 
In [53]: u, sigma, vh = np.linalg.svd(A) 
In [54]: u 
Out[54]: 
array([[-0.3246,  0.799 ,  0.5062], 
       [-0.7531,  0.1055, -0.6494], 
       [-0.5723, -0.592 ,  0.5675]]) 
In [55]: vh 
Out[55]: 
array([[-0.2114, -0.5539, -0.8053], 
       [ 0.4633, -0.7822,  0.4164], 
       [ 0.8606,  0.2851, -0.422 ]]) 
In [56]: sigma 
Out[56]: array([ 13.5824,   2.8455,   2.3287]) 

在此示例中，numpy.linalg.svd()返回了三个元组数组，我们将其解压缩为三个变量：u，sigma和vh，其中u代表A的左奇异向量（AA-1的特征向量），vh是A的右奇异向量（A-1A的特征向量的逆矩阵），sigma是A的非零奇异值（AA-1和A-1A的特征值）。 在该示例中，存在三个特征值，它们按顺序返回。 您可能会对结果感到怀疑，所以让我们做一些数学运算来验证它：

In [57]: diag_sigma = np.diag(sigma) 
In [58]: diag_sigma 
Out[58]: 
array([[ 13.5824,   0\.    ,   0\.    ], 
       [  0\.    ,   2.8455,   0\.    ], 
       [  0\.    ,   0\.    ,   2.3287]]) 
In [59]: Av = u.dot(diag_sigma).dot(vh) 
In [60]: Av 
Out[60]: 
array([[ 3.,  1.,  4.], 
       [ 1.,  5.,  9.], 
       [ 2.,  6.,  5.]]) 
In [61]: np.allclose(A, Av) 
Out[61]: True 

输入数组A可以转换为svd中的U ∑ V*，其中∑是奇异向量的值。 但是，从 NumPy 返回的sigma是具有非零值的数组，我们需要将其设为向量，因此在此示例中，形状为(3, 3)。 我们首先使用numpy.diag()制作sigma对角矩阵diag_sigma。 然后我们只需在u，diag_sigma和vh之间执行矩阵乘法，以检查计算结果（Av）是否与原始输入A相同，这意味着我们验证了 svd 结果。

QR 分解（有时称为极坐标分解）可用于任何M x N数组，并将其分解为正交矩阵（Q）和上三角矩阵（R）。 让我们尝试使用它来解决先前的Ax = b问题：

In [62]: b = np.array([1,2,3]).reshape(3,1) 
In [63]: q, r = np.linalg.qr(A) 
In [64]: x = np.dot(np.linalg.inv(r), np.dot(q.T, b)) 
In [65]: x 
Out[65]: 
array([[ 0.2667], 
       [ 0.4667], 
       [-0.0667]]) 

我们使用numpy.linalg.qr()分解A以获得q和r。 因此现在将原始方程式转换为(q * r) x = b。 我们可以使用r和q和b的逆矩阵乘法（点积）获得x。 由于q是一个单位矩阵，因此我们使用了转置而不是逆。 如您所见，结果x与我们使用矩阵和numpy.linalg.solve()时的结果相同; 这是解决线性问题的另一种方法。

注意

通常，三角矩阵逆的计算效率更高，因为您可以创建一个大型数据集并比较不同解决方案之间的性能。

多项式运算

NumPy 还提供了使用多项式的方法，并包括一个名为numpy.polynomial的包，用于创建，操作和拟合多项式。 常见的应用是内插和外推。 在本节中，我们的重点仍然是将ndarray与 NumPy 函数一起使用，而不是使用多项式实例。 （不用担心，我们仍将向您展示多项式类的用法。）

正如我们在矩阵类部分所述，将ndarray与 NumPy 函数结合使用是首选，因为ndarray可以在任何函数中接受，而矩阵和多项式对象则需要转换，尤其是在与其他程序通信时。 它们都提供了方便的属性，但是在大多数情况下，ndarray足够好。

在本节中，我们将介绍如何基于一组根来计算系数，以及如何求解多项式方程，最后我们将求值积分和导数。 让我们从计算多项式的系数开始：

In [66]: root = np.array([1,2,3,4]) 
In [67]: np.poly(root) 
Out[67]: array([  1, -10,  35, -50,  24]) 

numpy.poly()返回一阶多项式系数数组，其根是从最高到最低指数的给定数组root。 在此示例中，我们采用根数组[1,2,3,4]并返回多项式，等效于x^4 - 10x^3 + 35x^2 - 50x + 24。

我们需要注意的一件事是输入根数组应该是一维或正方形二维数组，否则会触发ValueError。 当然，我们也可以执行相反的操作：使用numpy.roots()根据系数计算根：

In [68]: np.roots([1,-10,35,-50,24]) 
Out[68]: array([ 4.,  3.,  2.,  1.]) 

现在，假设我们有方程式y = x^4 - 10x^3 + 35x^2 - 50x + 24，当x = 5时我们想知道y的值。 我们可以使用numpy.polyval()来计算：

In [69]: np.polyval([1,-10,35,-50,24], 5) 
Out[69]: 24 

numpy.polyval()具有两个输入参数，第一个是多项式的系数数组，第二个是用于求值给定多项式的特定点值。 我们也可以输入x的序列，结果将返回ndarray，其值对应于给定的x序列。

接下来我们将讨论积分和导数。 我们将继续以x^4 - 10x^3 + 35x^2 - 50x + 24的示例为例：

In [70]: coef = np.array([1,-10,35,-50,24]) 
In [71]: integral = np.polyint(coef) 
In [72]: integral 
Out[72]: array([  0.2 ,  -2.5 ,  11.6667, -25\.  ,  24\.  ,  0\.  ]) 
In [73]: np.polyder(integral) == coef 
Out[73]: array([ True,  True,  True,  True,  True], dtype=bool) 
In [74]: np.polyder(coef, 5) 
Out[74]: array([], dtype=int32) 

在此示例中，我们对整数演算使用numpy.polyint()，结果等于：

默认积分常数为 0，但我们可以使用输入参数k进行指定。 您可以自己做一些练习，以获得不同的k的积分。

让我们回到前面的示例-完成积分后，我们立即使用numpy.polyder()执行了微积分，并将导数与原始coef数组进行了比较。 我们得到了五个True布尔数组，它们验证了两个数组是否相同。

我们还可以在numpy.polyder()中指定区分的顺序（默认为 1）。 如我们所料，当我们计算四阶多项式的五阶导数时，它将返回一个空数组。

现在，我们将使用多项式类的实例重复这些示例，以查看用法的差异。 使用numpy.polynomial类的第一步是初始化多项式实例。 开始吧：

In [75]: from numpy.polynomial import polynomial 
In [76]: p = polynomial.Polynomial(coef) 
In [77]: p 
Out[77]: Polynomial([  1., -10.,  35., -50.,  24.], [-1,  1], [-1,  1]) 

请注意，在返回的p类型旁边是Polynomial类的实例，并且返回了三部分。 第一部分是多项式的系数。 第二个是domain，它表示多项式中的输入值间隔（默认为[-1, 1]）。 第三个是window，它根据多项式将域间隔映射到相应的间隔（默认也是[-1, 1]）：

In [78]: p.coef 
Out[78]: array([  1., -10.,  35., -50.,  24.]) 
In [79]: p.roots() 
Out[79]: array([ 0.25  ,  0.3333,  0.5   ,  1\.    ]) 

使用Polynomial实例，我们可以简单地调用coef属性以显示系数的ndarray。 roots()方法将显示根。 接下来，我们将求值特定值的多项式5：

In [80]: polynomial.polyval(p, 5) 
Out[80]: Polynomial([ 5.], [-1.,  1.], [-1.,  1.]) 

集成和派生还使用Polynomial类中的内置函数roots()完成，但函数名称更改为integ()和derive()：

In [81]: p.integ() 
Out[81]: Polynomial([  0\.    ,   1\.    ,  -5\.    ,  11.6667, -12.5   ,   4.8   ], [-1.,  1.], [-1.,  1.]) 
In [82]: p.integ().deriv() == p 
Out[82]: True 

多项式包还提供了特殊的多项式，例如 Chebyshev，Legendre 和 Hermite。 有关这些内容的更多详细信息，请参考这里。

总之，在大多数情况下，ndarray和 NumPy 函数可以解决与多项式有关的问题。 它们也是一种更可取的方式，因为程序中类型之间的转换较少，这意味着较少的潜在问题。 但是，在处理特殊多项式时，我们仍然需要多项式包。 我们几乎完成了数学部分。 在下一节中，我们将讨论线性代数的应用。

应用 -- 回归和曲线拟合

由于我们在讨论线性代数的应用，因此我们的经验来自实际案例。 让我们从线性回归开始。 因此，假设我们对一个人的年龄与其睡眠质量之间的关系感到好奇。 我们将使用 2012 年大不列颠睡眠调查在线提供的数据。

有 20,814 人参加了调查，年龄范围从 20 岁以下到 60 岁以上，他们通过 4 到 6 分对他们的睡眠质量进行了评估。

在这种情况下，我们将只使用 100 作为总人口，并模拟年龄和睡眠得分，其分布与调查结果相同。 我们想知道他们的年龄在增长，睡眠质量（分数）增加还是减少？ 如您所知，这是一个隐藏的线性回归实践。 一旦我们绘制了年龄和睡眠分数的回归线，通过观察该线的斜率，就可以得出答案。

但是在讨论应该使用哪个 NumPy 函数以及如何使用它之前，让我们首先创建数据集。 根据调查，我们知道 20 岁以下的参与者占 7%，21 岁至 30 岁的参与者占 24%，31 岁至 40 岁的参与者占 21%，60 岁以上的参与者占 21% 。因此，我们首先创建一组列表，代表每个年龄组的人数，并使用numpy.random.randint()模拟我们 100 个人口中的实际年龄，以查看年龄变量。 现在我们知道了每个年龄段的睡眠分数分布，我们称其为scores：这是[5.5, 5.7, 5.4, 4.9, 4.6, 4.4]的列表，[5.5, 5.7, 5.4, 4.9, 4.6, 4.4]是根据年龄段从最小到最大的顺序排列的。 在这里，我们还使用np.random.rand()函数以及均值（来自分数列表）和标准方差（均设置为0.01）来模拟分数分布（当然，如果您有一个好的数据集，则可以使用） ，最好只使用上一章介绍的numpy.genfromtxt()函数）：

In [83]: groups = [7, 24, 21, 19, 17, 12] 

In [84]: age = np.concatenate([np.random.randint((ind + 1)*10, (ind + 2)*10, group) for ind, group in enumerate(groups)]) 

In [85]: age 
Out[85]: 
array( 
[11, 15, 12, 17, 17, 18, 12, 26, 29, 24, 28, 25, 27, 25, 26, 24, 23,  27, 26, 24, 27, 20, 28, 20, 22, 21, 23, 25, 27, 24, 25, 35, 39, 33, 35, 30, 32, 32, 36, 38, 31, 35, 38, 31, 37, 36, 39, 30, 36, 33, 36, 37, 45, 41, 44, 48, 45, 40, 44, 42, 47, 46, 47, 42, 42, 42, 44, 40, 40, 47, 47, 57, 56, 53, 53, 57, 54, 55, 53, 52, 54, 57, 53, 58, 58, 54, 57, 55, 64, 67, 60, 63, 68, 65, 66, 63, 67, 64, 68, 66] 
) 
In [86]: scores = [5.5, 5.7, 5.4, 4.9, 4.6, 4.4] 
In [87]: sim_scores = np.concatenate([.01 * np.random.rand(group) + scores[ind] for ind, group in enumerate(groups)] ) 

In [88]: sim_scores 
Out[88]: 
array([ 
5.5089,  5.5015,  5.5024,  5.5   ,  5.5033,  5.5019,  5.5012, 
5.7068,  5.703 ,  5.702 ,  5.7002,  5.7084,  5.7004,  5.7036, 
5.7055,  5.7024,  5.7099,  5.7009,  5.7013,  5.7093,  5.7076, 
5.7029,  5.702 ,  5.7067,  5.7007,  5.7004,  5.7   ,  5.7017, 
5.702 ,  5.7031,  5.7087,  5.4079,  5.4082,  5.4083,  5.4025, 
5.4008,  5.4069,  5.402 ,  5.4071,  5.4059,  5.4037,  5.4004, 
5.4024,  5.4058,  5.403 ,  5.4041,  5.4075,  5.4062,  5.4014, 
5.4089,  5.4003,  5.4058,  4.909 ,  4.9062,  4.9097,  4.9014, 
4.9097,  4.9023,  4.9   ,  4.9002,  4.903 ,  4.9062,  4.9026, 
4.9094,  4.9099,  4.9071,  4.9058,  4.9067,  4.9005,  4.9016, 
4.9093,  4.6041,  4.6031,  4.6016,  4.6021,  4.6079,  4.6046, 
4.6055,  4.609 ,  4.6052,  4.6005,  4.6017,  4.6091,  4.6073, 
4.6029,  4.6012,  4.6062,  4.6098,  4.4014,  4.4043,  4.4013, 
4.4091,  4.4087,  4.4087,  4.4027,  4.4017,  4.4067,  4.4003, 
4.4021,  4.4061]) 

现在我们有了年龄和睡眠得分，每个变量都有 100 次事件。 接下来，我们将计算回归线：y = mx + c，其中y代表sleeping_score，而x代表age。 回归线的 NumPy 函数为numpy.linalg.lstsq()，它将系数矩阵和因变量值作为输入。 因此，我们要做的第一件事是将变量年龄打包到一个系数矩阵中，我们称之为AGE：

In [87]: AGE = np.vstack([age, np.ones(len(age))]).T 
In [88]: m, c = np.linalg.lstsq(AGE, sim_scores)[0] 
In [89]: m 
Out[90]: -0.029435313781 
In [91]: c 
Out[92]: 6.30307651938 

现在我们有斜率m和常数c。 我们的回归线是y = -0.0294x + 6.3031，这表明，随着年龄的增长，人的睡眠分数/质量会略有下降，如下图所示：

您可能认为回归线方程看起来很熟悉。 还记得我们在矩阵部分求解的第一个线性方程吗？ 是的，您也可以使用numpy.linalg.lstsq()来求解Ax = b方程，实际上这将是本章的第四个解决方案。 自己尝试； 用法与您使用numpy.linalg.solve()时非常相似。

但是，并非每个问题都能简单地通过绘制回归线来回答，例如按年的房价。 它显然不是线性关系，可能是平方或三次关系。 那么我们如何解决这个问题呢？ 让我们使用房价指数（国家统计局中的统计数据），然后选择 2004 年至 2013 年。我们将平均房价（英镑）调整为通胀因素； 我们想知道明年的平均价格。

在寻求解决方案之前，让我们首先分析问题。 问题的背后是多项式曲线拟合问题； 我们想找到最适合我们的问题的多项式，但是我们应该为它选择哪个 NumPy 函数？ 但是在此之前，让我们创建两个变量：每年的价格price和房屋的年份year：

In [93]: year = np.arange(1,11) 
In [94]: price = np.array([129000, 133000, 138000, 144000, 142000, 141000, 150000, 135000, 134000, 139000]). 
In [95]: year 
Out[95]: array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10]) 

现在我们有了年份和价格数据，我们假设它们之间的关系是平方的。 我们的目标是找到多项式：y = ax^2 + bx + c表示关系（一种典型的最小二乘法）。y代表price，x代表year。 这里我们将使用numpy.polyfit()帮助我们找到该多项式的系数：

In [97]: a, b, c = np.polyfit(year, price, 2) 
In [98]: a 
Out[98]: -549.242424242 
In [99]: b 
Out[99]: 6641.66666667 
In [100]: c 
Out[100]: 123116.666667 
In [101]: a*11**2 + b*11 + c 
Out[101]: 129716.66666666642 

我们从numpy.polyfit()获得了多项式的所有系数，它采用三个输入参数：第一个代表自变量：year; 第二个是因变量：price; 最后一个是多项式的阶数，在这种情况下为 2。现在我们只需要使用year = 11（从 2004 年起 11 年），就可以计算出估算价格。 您可以在下图中看到结果：

NumPy 可以从线性代数中获得许多应用，例如插值和外推，但是我们不能在本章中全部介绍它们。 我们希望本章为您使用 NumPy 解决线性或多项式问题提供一个良好的开端。

总结

在本章中，我们介绍了线性代数的矩阵类和多项式类。 我们研究了两个类提供的高级功能，还看到了ndarray在进行基本转置时的性能优势。 我们还介绍了numpy.linalg类，它提供了许多函数来处理ndarray的线性或多项式计算。

在本章中，我们做了很多数学运算，但同时也发现了如何使用 NumPy 帮助我们回答一些现实问题。

在下一章中，我们将了解傅立叶变换及其在 NumPy 中的应用。

六、NumPy 中的傅立叶分析

除其他事项外，傅立叶分析通常用于数字信号处理。 这要归功于它在将输入信号（时域）分离为以离散频率（频域）起作用的分量方面如此强大。 开发了另一种快速算法来计算离散傅里叶变换（DFT），这就是众所周知的快速傅里叶变换（FFT），它为分析及其应用提供了更多可能性。 NumPy 针对数字计算，也支持 FFT。 让我们尝试使用 NumPy 在应用上进行一些傅立叶分析！ 注意，本章假定不熟悉信号处理或傅立叶方法。

本章将涉及的主题是：

• 傅立叶分析的基础
• 一维和二维傅立叶变换
• 频谱密度估计
• 时频分析

开始之前

众所周知，傅里叶分析将函数表示为周期分量的总和（正弦和余弦函数的组合），并且这些分量能够恢复原始函数。 它在数字信号处理（例如滤波，插值等）中具有广泛的应用，因此我们在不提供 NumPy 的任何应用细节的情况下就不想谈论 NumPy 中的傅立叶分析。 为此，我们需要一个可视化的模块。

Matplotlib 是我们将在本章中使用的可视化模块。 请从官方网站下载并安装它。 或者，如果您正在使用 Anaconda 之类的 Scientific Python 发行版，则应该已经包括了 matplotlib。

我们将编写一个名为show()的简单显示函数，以帮助我们了解本章中的练习示例。 函数输出如下图所示：

上图显示原始函数（信号），下图显示傅立叶变换。 请在 IPython 命令提示符下键入以下代码，或将其保存到.py文件并将其加载到提示符：

##### The Plotting Functions ####import matplotlib.pyplot as plt 
import numpy as np 
def show(ori_func, ft, sampling_period = 5): 
    n = len(ori_func) 
    interval = sampling_period / n 
    plt.subplot(2, 1, 1) 
    plt.plot(np.arange(0, sampling_period, interval), ori_func, 'black') 
    plt.xlabel('Time'), plt.ylabel('Amplitude') 
    plt.subplot(2,1,2) 
    frequency = np.arange(n / 2) / (n * interval) 
    nfft = abs(ft[range(int(n / 2))] / n ) 
    plt.plot(frequency, nfft, 'red') 
    plt.xlabel('Freq (Hz)'), plt.ylabel('Amp. Spectrum') 
    plt.show() 

这是一个名为show()的显示函数，它具有两个输入参数：第一个是原始信号函数（ori_func），第二个是其傅里叶变换（ft）。 此方法将使用matplotlib.pyplot模块创建两个折线图：顶部带有黑线的原始信号，其中 x 轴表示时间间隔（在我们所有的示例中，我们设置了默认值，信号采样周期为 5 秒）， y 轴代表信号的幅度。 图表的下部是带有红线的傅里叶变换，其中 x 轴表示频率， y 轴代表振幅频谱。

在下一节中，我们将简单地介绍不同类型的信号波，并使用numpy.fft模块计算傅立叶变换。 然后我们调用show()函数以提供它们之间的视觉比较。

信号处理

在本节中，我们将使用 NumPy 函数来模拟多个信号函数并将其转换为傅立叶变换。 我们将重点介绍numpy.fft及其相关函数。 我们希望在本节之后，您将对在 NumPy 中使用傅立叶变换有所了解。 理论部分将在下一部分中介绍。

我们要使用的第一个示例是心跳信号，它是一系列正弦波。 频率为每分钟 60 次（1Hz），我们的采样周期为 5 秒长，采样间隔为 0.005 秒。 首先创建正弦波：

In [1]: time = np.arange(0, 5, .005) 
In [2]: x = np.sin(2 * np.pi * 1 * time) 
In [3]: y = np.fft.fft(x) 
In [4]: show(x, y) 

在此示例中，我们首先创建了采样时间间隔并将其保存到名为time的ndarray中。 然后，我们将time数组乘以2π并将其频率设为 1Hz 传递给numpy.sin()方法，以创建正弦波（x）。 然后将傅立叶变换应用于x并将其保存到y。 最后，我们使用预定义的方法show()与正弦波及其归一化的 Fourier 变换进行视觉比较，如下图所示：

上方的绿线代表心跳波； 因为我们使用 1Hz 持续 5 秒钟，所以获得了 5 个连续的正弦波。 这里要注意的一件事是，我们的采样间隔为 0.005 秒，这意味着我们使用 200 点（1 / 0.005）来模拟一个波形，因此它看起来相对平滑。 如果增加采样间隔（减少每个波的点数），我们将获得更强烈的正弦波。 图表的下部是基于频率（所谓的频谱）的标准化傅里叶变换的绝对值。 我们可以看到在 1Hz 处有一个高点，与我们的原始波频率相匹配。

接下来，我们将尝试计算多频正弦波并对其傅里叶变换进行计算。 在此之后，我们可能对傅立叶变换有了更清晰的了解。 以下代码显示了如何执行此操作：

In [8]: x2 = np.sin(2 * np.pi * 20 * time) 
In [9]: x3 = np.sin(2 * np.pi * 60 * time) 
In [10]: x += x2 + x3 
In [11]: y = np.fft.fft(x) 
In [12]: show(x, y) 

首先，我们再创建两个频率不同的正弦波，频率为 20Hz 的x2和 60Hz 的x3，然后将它们添加到原始的 1Hz 正弦波x中。 然后我们将修改后的x传递给傅立叶变换，并使用预定义的show()绘制图形。 您可以在下图中看到它：

在上方的绿色折线图中，我们可以看到正弦波组合了不同的频率，但实际上很难区分它们。 但是，在计算傅立叶变换并将其转换到频域之后，我们可以在下部的红色折线图中清楚地看到已识别出三个高点，分别是 1Hz，20Hz 和 60Hz。 这与我们原始的正弦波频率匹配。

从这两个示例中，您必须能够对傅立叶变换有所了解。 接下来，我们将演示另外三种信号处理：一种用于正方形信号，一种用于脉冲，另一种用于随机信号。

首先，我们使用numpy.zeros()以相同的时间间隔（time）创建方波信号。 我们希望方波频率为 10Hz，幅度为 1，因此我们将每 20 个时间间隔（200/10）设置为值 1，来模拟波浪并将其传递给傅立叶变换，如下面的代码块所示：

In [13]: x = np.zeros(len(time)) 
In [14]: x[::20] = 1 
In [15]: y = np.fft.fft(x) 
In [16]: show(x, y) 

此代码生成以下图形：

从上方的绿色折线图中，我们可以在 5 秒（10Hz）中看到 50 个连续的方波，但是当我们计算其傅立叶变换时，我们在频谱中获得了几个红色的高点，而不是 10Hz 时的一个红色高点。 您可能想知道方波是否也是周期函数，但是傅立叶变换为什么与正弦波有很大不同？ 请记住，傅立叶变换将时域转换为频域，但是在引擎盖下，有一系列正弦和余弦函数可以分解原始函数。 我们仍然可以看到红色的高点是规则间隔的，间隔为 10Hz。

接下来，我们将生成一个没有任何频率的单脉冲信号，并且将要计算其傅立叶变换。 将其与以前的方波进行比较，您可能会对傅里叶变换有更好的了解：

In [17]: x = np.zeros(len(time)) 
In [18]: x[380:400] = np.arange(0, 1, .05) 
In [19]: x[400:420] = np.arange(1, 0, -.05) 
In [20]: y = np.fft.fft(x) 
In [21]: show(x, y) 

首先我们创建一个与时间变量大小相同的全零ndarray，然后生成一个单脉冲信号，该信号在 2 秒时发生（x数组的第 400 个元素 ）。 我们在 2 秒左右的时间内占用了 40 个元素来模拟脉冲：20 个元素从 0 增加到 1，另一半从 1 减少到 0。我们将一个脉冲信号传递给傅里叶变换，并使用show()进行视觉比较。

下图中上部的绿色折线图是我们模拟的一脉冲信号，下部的红色折线图是其傅里叶变换。 我们可以看到下面图表中的最高点出现在等于 0 的频率上，这很有意义，因为我们的模拟信号没有任何频率。 但是，在零频率之后，我们仍然可以看到来自变换过程的两个不同频率的高点。

本节的最后一个示例是随机信号处理。 与前面的示例一样，我们还将 5 秒作为 100 个随机信号的总采样周期，该随机信号没有任何固定的频率与之关联。 然后，我们将随机信号传递给傅立叶变换，以获得其频域。 代码块如下：

In [22]: x = np.random.random(100) 
In [23]: y = np.fft.fft(x) 
In [24]: show(x, y) 

以下是由代码生成的图形：

看完这些示例之后，我们知道如何在 NumPy（简称为numpy.fft.fft()）中使用傅立叶变换-并且对傅立叶变换的外观有了一些了解。 在下一节中，我们将重点介绍理论部分。

傅立叶分析

定义 DFT 的方法很多。 但是，在 NumPy 实现中，DFT 定义为以下公式：

A[k]代表离散傅里叶变换，a[m]代表原始函数。 a[m] -> A[k]的转换是从配置空间到频率空间的转换。 让我们手动计算此方程，以更好地了解转换过程。 我们将使用具有 500 个值的随机信号：

In [25]: x = np.random.random(500) 
In [26]: n = len(x) 
In [27]: m = np.arange(n) 
In [28]: k = m.reshape((n, 1)) 
In [29]: M = np.exp(-2j * np.pi * k * m / n) 
In [30]: y = np.dot(M, x) 

在此代码块中，x是我们的模拟随机信号，其中包含 500 个值，并且在公式中对应于a[m]。 根据x的大小，我们计算出：

然后将其保存到M。 最后一步是使用M和x之间的矩阵乘法来生成 DFT 并将其保存到y中。

让我们通过将其与内置numpy.fft进行比较来验证结果：

In [31]: np.allclose(y, np.fft.fft(x)) 
Out[31]: True 

如我们所料，手动计算的 DFT 与 NumPy 内置模块相同。 当然，numpy.fft就像 NumPy 中的任何其他内置模块一样-已经过优化，并且已应用 FFT 算法。 让我们比较一下我们的手动 DFT 和numpy.fft的性能：

In [32]: %timeit np.dot(np.exp(-2j * np.pi * np.arange(n).reshape((n, 1)) * np.arange(n) / n), x) 
10 loops, best of 3: 18.5 ms per loop 
In [33]: %timeit np.fft.fft(x) 
100000 loops, best of 3: 10.9 µs per loop 

首先，我们将此方程式实现代码放在一行上，以测量执行时间。 我们可以看到它们之间的巨大性能差异。 在引擎盖下，NumPy 使用FFTPACK库执行傅立叶变换，该傅立叶变换在性能和准确性上都是非常稳定的库。

提示

如果您仍然觉得FFTPACK不够快，通常有一个FFTW库的性能要比FFTPACK好，但是从FFTPACK到FFTW的加速将不会那么快。

接下来，我们将计算逆 DFT。 iDFT 将频率序列映射回原始时间序列，该序列由以下公式定义：

我们可以看到反方程与 DFT 方程的不同之处在于指数参数的符号和通过1 / n进行归一化。 让我们再次进行手动计算。 由于指数参数的符号更改，我们可以重复使用先前代码中的m，k和n变量，只需重新计算M：

In [34]: M2 = np.exp(2j * np.pi * k * m / n) 
In [35]: x2 = np.dot(y, M2) / n 

再次，让我们用原始随机信号x验证计算出的逆 DFT 结果x2。 两者ndarray应该相同：

In [36]: np.allclose(x, x2) 
Out[36]: True 

当然，numpy.fft模块还支持逆 DFT 调用numpy.fft.ifft()以执行计算，如以下示例所示：

In [37]: np.allclose(x, np.fft.ifft(y)) 
Out[37]: True 

您可能会注意到，在前面的示例中，我们始终使用一维数组作为输入信号。 这是否意味着numpy.fft仅处理一维数据？ 当然不是; numpy.fft也可以处理二维或多维数据。 在开始这一部分之前，我们想谈谈返回的 FFT 数组的顺序和numpy.fft中的shift方法。

让我们创建一个包含 10 个随机整数的简单信号数组，并计算其傅里叶变换：

In [38]: a = np.random.randint(10, size = 10) 
In [39]: a 
Out[39]: array([7, 4, 9, 9, 6, 9, 2, 6, 8, 3]) 
In [40]: a.mean() 
Out[40]: 6.2999999999999998 
In [41]: A = np.fft.fft(a) 
In [42]: A 
Out[42]: 
array([ 63.00000000 +0.00000000e+00j,   
        -2.19098301 -6.74315233e+00j, 
        -5.25328890 +4.02874005e+00j, 
        -3.30901699 -2.40414157e+00j, 
        13.75328890 -1.38757276e-01j,    
      1.00000000 -2.44249065e-15j, 
        13.75328890 +1.38757276e-01j, 
     -3.30901699 +2.40414157e+00j, 
        -5.25328890 -4.02874005e+00j, 
     -2.19098301 +6.74315233e+00j]) 
In [43]: A[0] / 10 
Out[43]: (6.2999999999999998+0j) 
In [44]: A[int(10 / 2)] 
Out[44]: (1-2.4424906541753444e-15j) 

在此示例中，a是我们的原始随机信号，A是a的傅立叶变换。 当我们调用numpy.fft.fft(a)时，结果ndarray遵循“标准”顺序，其中第一个值A[0]包含零频率项（信号的均值）。 进行归一化处理（将其除以原始信号数组的长度A[0] / 10)时，我们得到的值与计算信号数组的平均值（a.mean()）时的值相同。

然后A[1:n/2]包含正频率项，A[n/2 + 1: n]包含负频率项。 在我们的示例中，当输入为偶数时，A[n/2]代表正数和负数。 如果要将零频率分量移到频谱中心，可以使用numpy.fft.fftshift() 例程。 请参见以下示例：

In [45]: np.fft.fftshift(A) 
Out[45]: 
array([  1.00000000 -2.44249065e-15j,   
     13.75328890 +1.38757276e-01j, 
        -3.30901699 +2.40414157e+00j, 
        -5.25328890 -4.02874005e+00j, 
        -2.19098301 +6.74315233e+00j, 
        63.00000000 +0.00000000e+00j, 
        -2.19098301 -6.74315233e+00j, 
     -5.25328890 +4.02874005e+00j, 
        -3.30901699 -2.40414157e+00j,   
     13.75328890 -1.38757276e-01j]) 

从此示例中，您可以看到numpy.fft.fftshift()交换了数组中的半空间，因此零频率分量移到了中间。 numpy.fft.ifftshift是反函数，将顺序移回“标准”。

现在，我们要谈谈多维 DFT。 让我们从二维开始。 您可能会看到以下等式与一维 DFT 非常相似，而第二维以明显的方式扩展。 多维 DFT 的思想是相同的，较高维中的逆函数也是如此。 您也可以尝试修改先前的代码，以将一维 DFT 计算为二维或多维 DFT，以更好地理解过程。 但是，现在我们仅要演示如何将numpy.fft用于二维和多维傅里叶变换：

In [46]: x = np.random.random(24) 
In [47]: x.shape = 2,12 
In [48]: y2 = np.fft.fft2(x) 
In [49]: x.shape = 1,2,12 
In [50]: y3 = np.fft.fftn(x, axes = (1, 2)) 

从这些示例中，您可以看到我们将numpy.fft.fft2()用于二维傅立叶变换，将numpy.fft.fftn()称为多维。 axes参数是可选的； 它指示要计算 FFT 的轴。 对于二维，如果未指定轴，则使用最后两个轴。 对于多维，模块使用所有轴。 在前面的示例中，我们仅应用了最后两个轴，因此傅立叶变换结果将与二维轴相同。 让我们来看看：

In [51]: np.allclose(y2, y3) 
Out[51]: True 

傅立叶变换应用

在前面的部分中，您学习了如何将numpy.fft用于一个一维和多维ndarray，并在幕后了解了实现细节。 现在是一些应用的时候了。 在本节中，我们将使用傅立叶变换进行一些图像处理。 我们将分析频谱，然后对图像进行插值以将其放大到两倍大小。 首先，让我们从 Packt Publishing 网站博客文章中下载练习图像。 将映像另存为本地目录scientist.png。

这是 RGB 图像，这意味着，当我们将其转换为ndarray时，它将是三维的。 为了简化练习，我们使用matplotlib中的图像模块读取图像并将其转换为灰度：

In [52]: from matplotlib import image 
In [53]: img = image.imread('./scientist.png') 
In [54]: gray_img = np.dot(img[:,:,:3], [.21, .72, .07]) 
In [55]: gray_img.shape 
Out[55]: (317L, 661L) 
In [56]: plt.imshow(gray_img, cmap = plt.get_cmap('gray')) 
Out[56]: <matplotlib.image.AxesImage at 0xa6165c0> 
In [57]: plt.show() 

您将得到以下图像作为结果：

预处理部分完成。 我们将图像读取到三维ndarray（img）中，并应用[亮度]公式使用0.21R + 0.72G + 0.07B将 RGB 图像转换为灰度图像 。 我们使用matplotlib中的pyplot模块显示灰度图像。 这里我们在图中未应用任何轴标签，但是从轴比例可以看到ndarray gray_img代表317 x 661像素的图像。

接下来，我们将进行傅立叶变换并显示频谱：

In [58]: fft = np.fft.fft2(gray_img) 
In [59]: amp_spectrum = np.abs(fft) 
In [60]: plt.imshow(np.log(amp_spectrum)) 
Out[60]: <matplotlib.image.AxesImage at 0xcdeff60> 
In [61]: plt.show() 

此代码将给出以下输出：

首先，我们对gray_img使用二维傅立叶变换，并使用对数刻度色图绘制幅度谱。 我们可以看到，由于零频率分量，拐角有所不同。 请记住，当我们使用numpy.fft.fft2()时，该顺序遵循标准的顺序，并且我们希望将零频分量置于中心。 因此，让我们使用shift例程：

In [62]: fft_shift = np.fft.fftshift(fft) 
In [63]: plt.imshow(np.log(np.abs(fft_shift))) 
Out[63]: <matplotlib.image.AxesImage at 0xd201dd8> 
In [64]: plt.show() 

这会将图像更改为：

现在我们可以看到零频率分量位于中心。 让我们转到本练习的最后一步：对图像进行插值以扩大尺寸。 我们在这里使用的技术非常简单。 我们将零频率插值到fft_shift数组中，并使它变成两倍大小。 然后我们将fft_shift逆转为标准阶数，并进行另一次逆转换回到原始域：

In [65]: m, n = fft_shift.shape 
In [66]: b = np.zeros((int(m / 2), n)) 
In [67]: c = np.zeros((2 * m - 1, int(n / 2))) 
In [68]: fft_shift = np.concatenate((b, fft_shift, b), axis = 0) 
In [69]: fft_shift = np.concatenate((c, fft_shift, c), axis = 1) 
In [70]: ifft = np.fft.ifft2(np.fft.ifftshift(fft_shift)) 
In [71]: ifft.shape 
Out[71]: (633L, 1321L) 
In [72]: ifft = np.real(ifft) 
In [73]: plt.imshow(ifft, cmap = plt.get_cmap('gray')) 
Out[73]: <matplotlib.image.AxesImage at 0xf9a0f98> 
In [74]: plt.show() 

在上一个代码块中，我们首先获取了fft_shift数组的形状（大小与gray_img相同）。 然后我们创建了两个零ndarrays并将它们沿四个方向填充到fft_shift数组中以将其放大。 因此，当我们将修改后的fft_shift数组逆回到标准阶数时，零频率将完美地位于中间。 当我们进行逆变换时，您可以看到形状已经加倍。 为了让pyplot模块绘制新数组，我们需要将数组转换为实数。 绘制新数组后，我们可以看到轴刻度是其大小的两倍。 而且，我们几乎不会丢失任何细节或图像模糊。 已使用傅立叶变换对图像进行插值。

总结

在本章中，我们介绍了一维和多维傅立叶变换的用法以及它们在信号处理中的应用方式。 现在，您了解了 NumPy 中离散傅立叶变换的实现，并且我们在手动实现的脚本与 NumPy 内置模块之间进行了性能比较。

我们还完成了图像插值的实际应用，并且由于了解matplotlib包的一些基础知识而获得了加号。

在下一章中，我们将看到如何使用numpy.distutils()子模块分发代码。

七、构建和分发 NumPy 代码

在现实世界中，您将编写一个应用，以将其分发到世界上或在其他各种计算机上重用。 为此，您希望应用以标准方式打包，以便社区中的每个人都能理解和遵循。 正如您现在已经注意到的那样，Python 用户主要使用名为pip的包管理器来自动安装其他程序员创建的模块。 Python 具有一个称为 PyPI（Python 包索引）的打包平台，该平台是 50,000 多个 Python 包的官方中央存储库。 一旦在 PyPi（又名 Cheese Shop）中注册了包，世界各地的其他用户都可以在使用pip等包管理系统对其进行配置后进行安装。 Python 随附了许多解决方案，可帮助您构建代码以准备分发给 Cheese Shop ，并且在本章中，我们将重点介绍两个此类工具，setuptools 和Distutils 除了这两个工具之外，我们还将研究 NumPy 提供的称为numpy.distutils的特定模块。 该模块使程序员更容易构建和分发特定于 NumPy 的代码。 该模块还提供了其他函数，例如用于编译 Fortran 代码，调用f2py等方法。 在本章中，我们将通过以下步骤来学习包装工作流程：

• 我们将建立一个小的但可行的设置
• 我们将说明将 NumPy 模块集成到您的设置中的步骤
• 我们将说明如何在互联网上注册和分发您的应用

Distutils 和 Setuptools 简介

在开始之前，首先让我们了解这些工具是什么以及为什么我们偏爱另一个工具。 Distutils是 Python 默认提供的框架，setuptools建立在标准Distutils的基础上，以提供增强的功能和特性。 在现实世界中，您将永远不会使用Distutils。 您可能想单独使用Distutils的唯一情况是setuptools不可用。 （良好的设置脚本应在继续之前检查setuptools的可用性。）在大多数情况下，用户最好安装setuptools，因为当今大多数包都是基于它们构建的。 在接下来的章节中，我们将使用setuptools来构建 Cython 代码; 因此，出于我们的目的，我们现在将安装setuptools并从现在开始广泛使用它。

接下来，让我们从安装必需的工具开始，以构建我们的第一个虚拟（但有效）安装程序。 安装程序正常运行后，我们将在 Pandas 脚本模块的真实脚本中深入介绍 NumPy 的更多功能。 我们将研究脚本中进行的检查，以使其更强大，以及在发生故障时如何提供更多信息。

准备工具

要在您的系统上安装setuptools，您需要先从这里下载系统中的ez_setup.py，然后从命令提示符处执行以下操作 ：

    $ python ez_setup.py

要测试setuptools的安装，请打开 Python shell 并键入以下内容：

> import setuptools 

如果前面的导入没有给出任何错误，则说明我们已成功安装setuptools。

建立第一个有效的发行版

我们前面提到的所有工具（setuptools，Distutils和numpy.distutils）都围绕setup函数。 为了了解大多数包装要求，我们将研究一个简单的setup函数，然后研究一个成熟的安装程序。 要创建基本的安装程序，我们需要使用有关包的元数据调用setup函数。 让我们叫第一个包py_hello，它只有一个函数greeter，并且在调用时只打印一条消息。 可从 Bitbucket 存储库下载该包。该项目的目录结构如下：

py_hello 
├── README 
├── MANIFEST.in 
├── setup.py 
├── bin 
│   └── greeter.bat 
└── greeter 
    ├── __init__.py 
    ├── greeter.py 

让我们在这里看一些标准文件：

• README - 此文件用于存储有关您的项目的信息。 该文件不是系统所需的文件，如果没有该文件，您仍将获得安装程序的构建，但是将其保留在此处是一种很好的做法。
• MANIFEST.in - 这是一个文本文件，Distutils使用该文本文件来收集项目中的所有文件。 这非常重要，只有此处列出的文件才会进入最终安装程序tar存档。 除了指定最终安装程序中应包含的文件之外，manifest还可以用于从项目目录中排除某些文件。 manifest文件是必需的； 如果不存在，则在使用setup.py时会出现错误。 如果您具有svn设置，则可以使用sdist命令通过解析.svn文件并构建manifest.in文件来自动包含文件。
• __init__.py - 该文件对于 Python 将该目录识别为模块很重要。 创建后可以将其留空。

要为此安装程序创建安装程序，我们在根目录中有setup.py，它使用setuptools中的setup函数：

from setuptools import setup 
import os 
description = open(os.path.join(os.path.dirname(__file__), 'README'), 'r').read() 
setup( 
    name = "py_hello", 
    packages = ["greeter"], 
    scripts = ["bin/greeter.bat"], 
    include_package_data = True, 
    package_data = { 
        "py_hello":[] 
        }, 
    version = "0.1.0", 
    description = "Simple Application", 
    author = "packt", 
    author_email = "packt@packt.com", 
    url = "https://bitbucket.org/tdatta/book/py_hello", 
    download_url = "https://bitbucket.org/tdatta/book/py_hello/zipball/master", 
    keywords = ["tanmay", "example_seutp", "packt"  "app"], 
    install_requires=[ 
        "setup >= 0.1"], 
    license='LICENSE', 
    classifiers = [ 
        "Programming Language :: Python", 
        "Development Status :: release 0.1", 
        "Intended Audience ::  new users", 
        "License :: Public", 
        "Operating System :: POSIX :: Linux", 
        "Topic :: Demo", 
        ], 
    long_description = description 
    ) 

以下是安装程序中使用的选项：

• name - 这是安装 TAR 归档文件的名称。
• packages - 这是一个列出要包含的包的列表。
• scripts - 这是要安装到/usr/bin.等标准位置的脚本的列表。在此特定情况下，仅存在一个 echo 脚本。 这样做的目的是向读者展示如何将包附带脚本。
• package_data - 这是字典，具有与文件列表相关联的键（包）。
• version - 这是您的项目的版本。 这将附加到安装程序名称的末尾。
• long_description - 在 PyPI 网站上显示时，它将转换为 HTML。 它应该包含有关您的项目打算提供的信息。 您可以直接在脚本中编写它； 但是，最佳实践是维护README文件并从此处读取说明。
• install_required - 这是用于添加安装依赖项的列表。 您将添加代码中使用的第三方模块的名称和版本。 请注意遵循约定以在此处指定版本。
• classifiers - 当您在 PyPI 网站上上传包时，将选中此选项。 您应该从以下网站提供的选项中进行选择。

现在，使用build选项运行setup.py应该不会给您任何错误，并生成带有.egg-info后缀的文件夹。 此时，您可以使用sdist选项运行setup.py，并创建一个可以与世界共享的包。

您应该看到最终消息为创建 tar 归档文件，如下所示：

要测试该包，可以按照以下步骤将其安装在本地计算机上：

python setup.py install 

并按以下所示检查它：

这时，在cmd/bash提示符下写greeter，您将看到一条消息does nothing。 该回显消息来自greeter.bat，我们将其放置在安装文件的scripts键中。

下一部分可以添加到此框架setup.py以包括NumPy特定的功能。

添加 NumPy 和非 Python 源代码

接下来，我们将研究一些特定于 NumPy 的代码，并了解如何提高设置的错误处理能力； 通常，我们将探索一些良好的编程习惯。 我们还将展示如何将非 Python 源（c,fortran或f2py）添加到安装程序中。 以下分析显示了完整代码的一部分，您可以在随附的代码文件中或在这个页面中找到这些代码：

if sys.version_info[0] < 3: 
   import __builtin__ as builtins 
else: 
   import builtins 

..... 
..... 
..... 
For full sample look for setup.py file with the accompanying CD
..... 
..... 

##define a function to import numpy if available and return true else false 

def is_numpy_installed(): 
   try: 
        import numpy 
   except ImportError: 
        return False 
   return True  
## next we will import setuptools feature here  
## We need to do this here because setuptools will "Monkey patch" the setup function 
## 
SETUPTOOLS_COMMANDS = set([ 
    'develop', 'release', 'bdist_egg','bdist_rpm', 
    'bdist_wininst', 'install_egg_info', 'build_sphinx', 
    'easy_install', 'upload', 'bdist_wheel', 
    '--single-version-externally-managed', 
]) 

if SETUPTOOLS_COMMANDS.intersection(sys.argv): 
    import setuptools 
    extra_setuptools_args = dict( 
        zip_safe-False # Custom clean command to remove build artifacts 

## The main function where we link everything 

def setup_package(): 

   # check NumPy and raise excpetions 

    if is_numpy_installed() is False: 
        raise ImportError("Numerical Python (NumPy) is not installed. The package requires it be installed. Installation instruction available at the NumPy website") 

    from numpy.distutils.core import setup, Extension 

    # add extension from Fortran 

    ext1 = Extension(name = "firstExt", 
                     sources = ['firstExt.f']) 
    ext2  = Extension(name = "convolutedExt", 
                      sources = ['convolutedExt.pyf, stc2.f'], 
                      include_dir = ['paths to include'], 
                      extra_objects = "staticlib.a") 

    metadata = dict(name = "yourPackage", 
                    description="short desc", 
                    license = "licence info here", 
                    ext_modules = [ext1, ext2] 
                    .. 
                    # metadata as we set previously 
                    .. 
                    **extra_setuptools_args 
    ) 

    setup(**metadata) 

if __name__ == "__main__": 
    setup_package() 

上面的脚本从完整的工作设置中删除，着重于几乎所有设置脚本中都可以找到的某些方面。 这些任务确保您已完成足够的错误处理，并且脚本在不解释/提示下一步操作的情况下不会失败：

1. 检查是否已安装 NumPy。 此处用来确保已安装 NumPy 的模式是一种标准模式，您可以将其用于计划使用的所有模块，并且是安装程序所必需的。 为了执行此任务，我们首先构建一个函数is_numpy_installed尝试导入numpy并返回一个布尔值。 您可能会为安装文件可能使用的所有外部包创建类似的函数。 高级用户可以使用 Python 装饰器来以更优雅的方式进行处理。 如果此函数返回错误值，则安装程序应输出警告/信息，以防没有此包无法完成安装。
2. 将Extensions添加到设置文件中。
3. Extension类是使我们能够向安装程序中添加非 Python 代码的对象。 sources参数可能包含 Fortran 源文件列表。 但是，列表源可能最多包含一个f2py签名文件，然后扩展模块的名称必须与签名文件中使用的<module>匹配。 f2py签名文件必须恰好包含一个 Python 模块块，否则安装程序将无法构建。 您可以决定不在sources参数中添加签名文件。 在这种情况下，f2py将扫描 Fortran 源文件以获取常规签名，以构造 Fortran 代码的包装器。 可以使用Extension类参数f2py_options来指定f2py进程的其他选项。 这些选项不在本书的讨论范围内，大多数读者不会使用它们。 有关更多详细信息，用户可以参考api文档中的numpy.distutils扩展类。

可以按以下方式测试安装文件：

 $ python <setup.py file> build_src build_ext --help 

这里的build_src参数用于构造 Fortran 包装器扩展模块。 这里假定用户在其计算机上安装了 C/C++ 和 Fortran 编译器。

测试您的包

非常重要的一点是，所构建的包可以在用户的​​计算机上正常运行/安装。 因此，您应该花时间测试包。 测试安装背后的总体思路是创建一个 VirtualEnv 并尝试安装该包或完全使用另一个系统。 在此阶段遇到的任何错误都应删除，并且作者应尝试确保更容易遵循这些异常。 异常也应尝试提供解决方案。 此阶段的常见错误是：

• 关于预装模块和库的假设。
• 开发人员可能会忘记在安装文件中包含依赖项。 如果使用新的 VirtualEnv 来测试安装程序，则会捕获此错误。
• 权限和提升权限的要求。
• 某些用户可能对计算机具有只读访问权限。 这很容易被忽略，因为大多数开发人员在自己的机器上都没有这种情况。 如果包的提供者遵循正确的方法来选择要写入的目录，则应该不会出现此问题。 通常，通过使用没有管理员访问权限的用户测试脚本来检查这种情况是一种很好的做法。

分发您的应用

完成模块/应用的所有开发并准备好完整的正常工作的应用和设置文件后，下一个任务就是与世界分享您的辛勤工作，使他人受益。 使用 PyPI 将其发布到全世界的步骤非常简单。 作为包作者，您需要做的第一件事就是注册自己。 您可以直接从命令行执行以下操作：

 **$ python setup.py register
    running register
    running egg_info
    ....
    ....
    We need to know who you are, so please choose either:
     1\. use your existing login,
     2\. register as a new user,
     3\. have the server generate a new password for you (and email it to you), or
     4\. quit
    Your selection [default 1]:** 

提示

如果setup.py中缺少任何文件的正确元数据信息，则此过程将失败。 确保首先工作setup.py。

最后，您可以通过执行以下操作在 PyPI 上上传您的发行版：

 $ python setup.py sdist upload 

希望，如果您正确键入了所有内容，您的应用将被打包并在 PyPI 上供世界使用。

总结

在本章中，我们介绍了用于打包和分发应用的工具。 我们首先看了一个更简单的setup.py文件。 您研究了setup函数的属性以及这些参数如何链接到最终安装程序。 接下来，我们添加了与 NumPy 相关的代码，并添加了一些异常处理代码。 最后，我们构建了安装程序并学习了如何在 Cheese Shop（PyPI 网站）上上传它。 在下一章中，您将研究通过将 Python 代码的一部分转换为 Cython 来进一步加速 Python 代码的方法。

八、使用 Cython 加速 NumPy

Python 与 NumPy 库相结合为用户提供了编写高度复杂的函数和分析的工具。 随着代码的大小和复杂性的增长，代码库中的低效率问题开始蔓延。一旦项目进入完成阶段，开发人员就应开始关注代码的性能并分析瓶颈。 Python 提供了许多工具和库来创建优化且性能更快的代码。

在本章中，我们将研究一种名为 Cython 的工具。 Cython 是 Python 和“Cython”语言的静态编译器，在从事科学库/数值计算的开发人员中特别流行。 许多用 Python 编写的著名分析库都大量使用 Cython（Pandas，SciPy，scikit-learn 等）。

Cython 编程语言是 Python 的超集，用户仍然喜欢 Python 所提供的所有功能和更高层次的结构。 在本章中，我们将研究 Cython 起作用的许多原因，并且您将学习如何将 Python 代码转换为 Cython。 但是，本章不是 Cython 的完整指南。

在本章中，我们将介绍以下主题：

• 在我们的计算机上安装 Cython
• 将少量 Python 代码重写为 Cython 版本并进行分析
• 学习在 Cython 中使用 NumPy

优化代码的第一步

每个开发人员在优化其代码时应注意的问题如下：

• 您的代码执行多少个函数调用？
• 有多余的调用吗？
• 该代码使用了多少内存？
• 是否存在内存泄漏？
• 瓶颈在哪里？

前四个问题主要由分析器工具回答。 建议您至少学习一种分析工具。 分析工具将不在本章中介绍。 在大多数情况下，建议先尝试优化函数调用和内存使用，然后再使用低级方法，例如 Cython 或汇编语言（使用 C 衍生语言）。

一旦确定了瓶颈，并且解决了算法和逻辑的所有问题，Python 开发人员便可以进入 Cython 的世界，以提高应用的速度。

设置 Cython

Cython 是一个将类型定义的 Python 代码转换为 C 代码的编译器，该代码仍在 Python 环境中运行。 最终输出是本机代码，其运行速度比 Python 生成的字节码快得多。 在大量使用循环的代码中，Python 代码加速的幅度更加明显。 为了编译 C 代码，首要条件是在计算机上安装 C/C++ 编译器，例如gcc（Linux）或mingw（Windows）。

第二步是安装 Cython。 Cython 与其他带有 Python 模块的库一样，可以使用任何首选的方法（PIP，EasyInstall 等）进行安装。 完成这两个步骤后，您可以通过尝试从 Shell 调用 Cython 来测试设置。 如果收到错误消息，则说明您错过了第二步，需要重新安装 Cython 或从 Cython 官方网站下载 TAR 归档文件，然后从这次下载的root文件夹中运行以下命令：

 python setup.py install 

正确完成所有操作后，您可以继续使用 Cython 编写第一个程序。

Cython 中的 Helloworld

Cython 程序看起来与 Python 程序非常相似，但大多带有附加的类型信息。 让我们看一个简单的程序，该程序计算给定n的第n个斐波那契数：

defcompute_fibonacchi(n): 
    """ 
    Computes fibonacchi sequence 

    """ 

    a = 1 
    b = 1 
    intermediate = 0  
    for x in xrange(n): 
intermediate = a 
        a = a + b 
        b = intermediate 
    return a 

让我们研究一下该程序，以了解在调用带有某些数字输出的函数时幕后的情况。 假设compute_fibonacchi(3)。

众所周知，Python 是一种解释性和动态语言，这意味着您无需在使用变量之前声明变量。 这意味着在函数调用开始时，Python 解释器无法确定n将保留的值的类型。 当您使用某个整数值调用函数时，Python 会通过名为装箱和拆箱的过程自动为您进行类型推断。

在 Python 中，一切都是对象。 因此，当您输入1或hello时，Python 解释器将在内部将其转换为对象。 在许多在线材料中，此过程也称为拳击。 该过程可以可视化为：

那么当您将函数应用于对象时会发生什么呢？ Python 解释器必须做一些额外的工作来推断类型并应用函数。 在一般意义上，下图说明了add函数在 Python 中的应用。 Python 是一种解释型语言，它在优化函数调用方面做得并不出色，但是可以使用 C 或 Cython 很好地优化它们：

这种装箱和拆箱不是免费的，需要花费宝贵的计算时间。 当这样的操作被循环执行多次时，效果变得更加显着。

在n = 20上运行时，以下程序在 IPython 笔记本上每个循环大约需要 1.8 微秒：

现在让我们将该程序重写为 Cython：

defcompute_fibonacchi_cython(int n): 
cdefint a, b, intermediate, x 
    a, b= 1, 1 
    intermediate, x  = 0, 0 
    for x in xrange(n): 
        intermediate = a 
        a = a+b 
        b = intermediate 
    return a 

该程序每个循环花费64.5纳秒：

提示

尽管在此示例代码中提高速度非常重要，但这不是您将遇到的实际代码，因此您应始终记住首先在代码上运行分析器并确定需要优化的部分。 同样，在使用 Cython 时，开发人员应考虑在使用静态类型和灵活性之间进行权衡。 使用类型会降低灵活性，有时甚至会降低可读性。

通过删除xrange并改用for循环，可以进一步改进此代码。 当您对模块的所有组件/功能都满意并且没有错误后，用户可以将这些函数/过程存储在扩展名为.pyx的文件中。 这是 Cython 使用的扩展名。 将此代码与您的应用集成的下一步是在安装文件中添加信息。

在这里，出于说明目的，我们将代码存储在名为fib.pyx的文件中，并创建了一个构建该模块的安装文件：

from distutils.core import setup, Extension 
from Cython.Build import cythonize 
from Cython.Distutils import build_ext 

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

在这里，请注意扩展名first 的名称与模块的名称完全匹配。 如果您无法保持相同的名称，则将收到一个神秘的错误：

多线程代码

您的应用可能会使用多线程代码。 由于全局解释器锁（GIL），Python 不适合多线程代码。 好消息是，在 Cython 中，您可以显式解锁 GIL，并使您的代码真正成为多线程。 只需在您的代码中放置一个with nogil:语句即可。 您以后可以在代码中使用with gil获取 GIL：

with nogil: 
<The code block here> 
function_name(args) with gil:  
<function body> 

NumPy 和 Cython

Cython 具有内置支持，可提供对 NumPy 数组的更快访问。 这些功能使 Cython 成为优化 NumPy 代码的理想人选。 在本节中，我们将研究用于计算欧式期权价格的代码，欧式期权是一种使用蒙特卡洛技术的金融工具。 不期望有金融知识； 但是，我们假设您对蒙特卡洛模拟有基本的了解：

defprice_european(strike = 100, S0 = 100, time = 1.0,  
rate = 0.5, mu = 0.2, steps = 50,  
N = 10000, option = "call"): 

dt = time / steps 
rand = np.random.standard_normal((steps + 1, N)) 
S = np.zeros((steps+1, N)); 
S[0] = S0 

for t in range(1,steps+1): 
S[t] = S[t-1] * np.exp((rate-0.5 * mu ** 2) * dt 
+ mu * np.sqrt(dt) * rand[t]) 
price_call = (np.exp(-rate * time) 
* np.sum(np.maximum(S[-1] - strike, 0))/N) 
price_put = (np.exp(-rate * time) 
* np.sum(np.maximum(strike - S[-1], 0))/N) 

returnprice_call if option.upper() == "CALL" else price_put 

以下是前面示例的 Cython 化代码：

import numpy as np 
def price_european_cython(double strike = 100,doubleS0 = 100, 
                          double time = 1.0, double rate = 0.5, 
                          double mu = 0.2, int steps = 50,  
                   long N = 10000, char* option = "call"): 
   cdef double dt = time / steps 
   cdefnp.ndarray rand = np.random.standard_normal((steps + 1, N)) 
   cdefnp.ndarray S = np.zeros([steps+1, N], dtype=np.float) 
       #cdefnp.ndarrayprice_call = np.zeroes([steps+1,N],     dtype=np.float) 
       S[0] = S0 

   for t in xrange(1,steps+1): 
           S[t] = S[t-1] * np.exp((rate-0.5 * mu ** 2) * dt 
                               + mu * np.sqrt(dt) * rand[t]) 
           price_call = (np.exp(-rate * time) 
                  * np.sum(np.maximum(S[-1] - strike, 0))/N) 
           price_put = (np.exp(-rate * time) 
                  * np.sum(np.maximum(strike - S[-1], 0))/N) 

   return price_call if option.upper() == "CALL" else price_put 

与此相关的安装文件如下所示：

from distutils.core import setup, Extension 
from Cython.Build import cythonize 
from Cython.Distutils import build_ext 
import numpy.distutils.misc_util 

include_dirs = numpy.distutils.misc_util.get_numpy_include_dirs() 

setup( 

name="numpy_first", 
    version="0.1", 
ext_modules=[Extension('dynamic_BS_MC', 
                           ['dynamic_BS_MC.pyx'], 
include_dirs = include_dirs)], 
cmdclass={'build_ext': build_ext} 
) 

虽然通过 Cython 代码获得的加速效果非常好，并且您可能会倾向于在 Cython 中编写大多数代码，但建议仅将性能至关重要的部分转换为 Cython。 NumPy 在优化对数组的访问和执行更快的计算方面做得非常出色。 该代码可以视为描述该代码的理想候选者。 前面的代码有很多“松散的结果”，可以当作练习来解决 Python 中的性能问题，并在采用 Cython 方式之前先最佳地使用 NumPy。 由于 Python 的动态特性，盲目地对 NumPy 代码进行 Cython 化的速度提升可能不如具有真正问题的最优编写代码那样快。

最后，我们介绍在 Cython 中开发模块时应遵循的以下内容：

1. 用纯 Python 编写代码并进行测试。
2. 运行分析器并确定要关注的关键区域。
3. 创建一个新模块以保存 Cython 代码（<module_name>.pyx）。
4. 将这些区域中的所有变量和循环索引转换为它们的 C 对应物。
5. 使用以前的测试设置进行测试。
6. 将扩展添加到安装文件中。

总结

在本章中，我们了解了如何将 Python 代码隐蔽到 Cython 中。 我们还研究了一些涉及 NumPy 数组的示例 Python 代码。 我们简要介绍了 Python 语言中装箱和拆箱的概念以及它们如何影响代码性能。 我们还说明了如何显式解锁臭名昭著的 GIL。 为了进一步深入研究 Cython，我们建议《Cython 编程学习手册》，Philip Herron，Packt Publishing。 在下一章中，您将了解 NumPy C API 以及如何使用它。

九、NumPy C-API 简介

NumPy 是一个通用库，旨在满足科学应用开发人员的大多数需求。 但是，随着应用的代码库和覆盖范围的增加，计算也随之增加，有时用户需要更具体的操作和优化的代码段。 我们已经展示了 NumPy 和 Python 如何具有诸如 F2PY 和 Cython 之类的工具来满足这些需求。 这些工具可能是将函数重写为本地编译代码以提高速度的绝佳选择。 但是在某些情况下（利用 C 库，例如 NAG 编写一些分析），您可能想做一些更根本的事情，例如为您自己的库专门创建新的数据结构。 这将要求您有权访问 Python 解释器中的低级控件。 在本章中，我们将研究如何使用 Python 及其扩展名 NumPy C-API 提供的 C-API 进行此操作。 C-API 本身是一个非常广泛的主题，可能需要一本书才能完全涵盖它。 在这里，我们将提供简短的介绍和示例，以帮助您开始使用 NumPy C-API。

本章将涉及的主题是：

• Python C-API 和 NumPy C-API
• 扩展模块的基本结构
• 一些特定于 NumPy 的 C-API 函数的简介
• 使用 C-API 创建函数
• 创建一个可调用的模块
• 通过 Python 解释器和其他模块使用模块

Python 和 NumPy C-API

我们使用的 Python 实现是 Python 解释器的基于 C 的实现。 NumPy 专用于此基于 C 的 Python 实现。 Python 的此实现带有 C-API，它是解释器的基础，并向其用户提供低级控制。 NumPy 通过提供丰富的 C-API 进一步增强了这一功能。

用 C/C++ 编写函数可以为开发人员提供灵活性，以利用这些语言提供的一些高级库。 但是，就必须在解析输入周围编写太多样板代码以构造返回值而言，代价显而易见。 此外，开发人员在引用/解引用对象时必须格外小心，因为这最终可能会导致讨厌的错误和内存泄漏。 随着 C-API 的不断发展，还存在代码未来兼容性的问题。 因此，如果开发人员想要迁移到更高版本的 Python，则他们可能需要为这些基于 C-API 的扩展进行大量维护工作。 由于这些困难，大多数开发人员选择尝试其他优化技术。 （例如 Cython 或 F2PY），然后再探索这条路径。 但是，在某些情况下，您可能想重用 C/C++ 中的其他现有库，这可能适合您的特定目的。 在这些情况下，最好为现有函数编写包装并公开 Python 项目。

接下来，我们将看一些示例代码，并在本章继续介绍时解释其关键函数和宏。 此处提供的代码与 Python 2.X 版本兼容，可能不适用于 Python 3.X。 但是，转换过程应该相似。

提示

开发人员可以尝试使用称为 cpychecker 的工具来检查模块中的引用计数时的常见错误。 请访问这里了解更多详细信息。

扩展模块的基本结构

用 C 编写的扩展模块将包含以下部分：

• 标头段，其中包含所有外部库和Python.h
• 初始化段，您可以在其中定义模块名称和 C 模块中的函数
• 方法结构数组，用于定义模块中的所有函数
• 一个实现部分，您在其中定义要公开的所有函数

标头段

标题片段是非常标准的，就像普通的 C 模块一样。 我们需要包括Python.h头文件，以使我们的 C 代码可以访问 C-API 的内部。 该文件位于<path_to_python>/include中。 我们将在示例代码中使用数组对象，因此我们也包含了numpy/arrayobject.h头文件。 我们不需要在此处指定头文件的完整路径，因为路径解析是在setup.py中处理的，我们将在后面进行介绍：

/* 
Header Segment 
*/ 

##include <Python.h> 
##include <math.h> 
##include <numpy/arrayobject.h> 
Initialization Segment 

初始化段

初始化段从以下内容开始：

1. 调用PyMODINIT_FUNC宏。 此宏在 Python 标头中定义，并且在开始定义模块之前总是会被调用。
2. 下一行定义了初始化函数，并在加载该函数时由 Python 解释器调用。 函数名称必须为init<module_name>格式，C 代码将要公开的模块和函数的名称。

该函数的主体包含对Py_InitModule3的调用，该调用定义模块的名称和模块中的函数。 该函数的一般结构如下：

(void)Py_InitModule3(name_of_module, method_array, Docstring) 

对import_array()的最终调用是特定于 NumPy 的函数，如果您的函数正在使用 Numpy 数组对象，则需要此函数。 这样可以确保加载 C-API，以便如果您的 C++ 代码使用 C-API，则 API 表可用。 未能调用此函数和使用其他 NumPy API 函数将很可能导致分段错误错误。 建议您阅读 NumPy 文档中的import_array()和import_ufunc()：

/* 
Initialization module 
*/ 

PyMODINIT_FUNC 
initnumpy_api_demo(void) 
{ 
(void)Py_InitModule3("numpy_api_demo", Api_methods, 
         "A demo to show Python and Numpy C-API"); 
import_array(); 
} 

方法结构数组

在此部分中，您将定义模块将要公开给 Python 的方法数组。 我们在这里定义了两个函数以求其平方。 一种方法将普通的 Python double值作为输入，第二种方法对 Numpy 数组进行操作。 PyMethodDef结构可以在 C 中定义如下：

Struct PyMethodDef { 
char *method_name; 
PyCFunction method_function; 
int method_flags; 
char *method_docstring; 
}; 

这是此结构的成员的描述：

• method_name：函数的名称在此处。 这将是函数向 Python 解释器公开的名称。
• method_function：此变量保存在 Python 解释器中调用method_name时实际调用的 C 函数的名称。
• method_flags：这告诉解释器我们的函数正在使用三个签名中的哪个。 该标志的值通常为METH_VARARGS。 如果要允许关键字参数进入函数，可以将该标志与METH_KEYWORDS组合。 它也可以具有METH_NOARGS的值，这表明您不想接受任何参数。
• method_docstring：这是函数的文档字符串。

该结构需要以一个由NULL和 0 组成的标记终止，如以下示例所示：

/* 
Method array structure definition 
*/ 
static PyMethodDefApi_methods[] = 
{ 
{"py_square_func", square_func, METH_VARARGS, "evaluate the squares"}, 
{"np_square", square_nparray_func, METH_VARARGS,  "evaluates the square in numpy array"}, 
{NULL, NULL, 0, NULL} 
}; 

实现部分

实现部分是最直接的部分。 这就是方法的 C 定义所要去的地方。 在这里，我们将研究两个函数来平方它们的输入值。 这些函数的复杂度保持在较低水平，以便您专注于方法的结构。

使用 Python C-API 创建数组平方函数

Python 函数将对自身的引用作为第一个参数，然后是赋予该函数的真实参数。 PyArg_ParseTuple函数用于将 Python 函数中的值解析为 C 函数中的局部变量。 在此函数中，我们将值强制转换为双精度，因此我们将d用作第二个参数。 您可以在这个页面上查看此函数接受的字符串的完整列表。

使用Py_Buildvalue返回计算的最终结果，它使用类似类型的格式字符串从您的答案中创建 Python 值。 我们在这里使用f表示浮点数，以证明对double和float的处理方式类似：

/* 
Implementation of the actual C funtions 
*/ 

static PyObject* square_func(PyObject* self, PyObject* args) 
{ 
double value; 
double answer; 

/*  parse the input, from python float to c double */ 
if (!PyArg_ParseTuple(args, "d", &value)) 
return NULL; 
/* if the above function returns -1, an appropriate Python exception will 
* have been set, and the function simply returns NULL 
*/ 

answer = value*value; 

return Py_BuildValue("f", answer); 
} 

使用 NumPy C-API 创建数组平方函数

在本节中，我们将创建一个函数以对 NumPy 数组的所有值求平方。 这里的目的是演示如何在 C 语言中获取 NumPy 数组，然后对其进行迭代。 在现实世界中，可以使用映射或通过向量化平方函数以更简单的方式完成此操作。 我们正在使用与O!格式字符串相同的PyArg_ParseTuple函数。 该格式字符串具有(object) [typeobject, PyObject *]签名，并以 Python 类型对象作为第一个参数。 用户应阅读官方 API 文档，以查看允许使用其他格式的字符串以及哪种字符串适合他们的需求：

注意

如果传递的值的类型不同，则引发TypeError。

以下代码段说明了如何使用PyArg_ParseTuple解析参数。

// Implementation of square of numpy array 

static PyObject* square_nparray_func(PyObject* self, PyObject* args) 
{ 

// variable declarations 
PyArrayObject *in_array; 
PyObject      *out_array; 
NpyIter *in_iter; 
NpyIter *out_iter; 
NpyIter_IterNextFunc *in_iternext; 
NpyIter_IterNextFunc *out_iternext; 

// Parse the argument tuple by specifying type "object" and putting the reference in in_array 
if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &in_array)) 
return NULL; 
...... 
...... 

下一步是创建一个数组以存储其输出值和迭代器，以便在 Numpy 数组上进行迭代。 请注意，创建对象时，每个步骤都有一个{handle failure}代码。 这是为了确保如果发生任何错误，我们可以通过调试来确定错误代码的位置：

//Construct the output from the new constructed input array 
out_array = PyArray_NewLikeArray(in_array, NPY_ANYORDER, NULL, 0); 
// Test it and if the input is nothing then just return nothing. 
{handle failure} 

//  Create the iterators 
in_iter = NpyIter_New(in_array, NPY_ITER_READONLY, NPY_KEEPORDER, 
NPY_NO_CASTING, NULL); 

// {handle failure} 

out_iter = NpyIter_New((PyArrayObject *)out_array, NPY_ITER_READWRITE, 
NPY_KEEPORDER, NPY_NO_CASTING, NULL); 
{handle failure} 

in_iternext = NpyIter_GetIterNext(in_iter, NULL); 
out_iternext = NpyIter_GetIterNext(out_iter, NULL); 
{handle failure} 

double ** in_dataptr = (double **) NpyIter_GetDataPtrArray(in_iter); 
double ** out_dataptr = (double **) NpyIter_GetDataPtrArray(out_iter); 

A simple handle failure module is like 
// {Start handling failure} 
if (in_iter == NULL) 
// remove the ref and return null 
Py_XDECREF(out_array); 
return NULL; 
// {End handling failure} 

看了前面的样板代码之后，我们终于来到了发生所有实际动作的部分。 那些熟悉 C++ 的人会发现迭代方法与向量迭代相似。 我们之前定义的in_iternext函数在这里派上用场，用于迭代 Numpy 数组。 在while循环之后，我们确保在两个迭代器上都调用了NpyIter_Deallocate，在输出数组上调用了Py_INCREF; 未能调用这些函数是导致内存泄漏的最常见错误类型。 内存泄漏问题通常非常微妙，通常在具有长时间运行的代码（例如服务或守护程序）时才会出现。 要抓住这些问题，不幸的是，没有比使用调试器更深入的方法容易的方法了。 有时，只需要编写几个printf语句即可输出总内存使用情况：

/*  iterate over the arrays */ 
do { 
    out_dataptr =pow(**in_dataptr,2); 
} while(in_iternext(in_iter) && out_iternext(out_iter)); 

/*  clean up and return the result */ 
NpyIter_Deallocate(in_iter); 
NpyIter_Deallocate(out_iter); 
Py_INCREF(out_array); 
return out_array; 

构建和安装扩展模块

成功编写函数后，下一步是构建模块并在我们的 Python 模块中使用它。 setup.py文件看起来像以下代码片段：

from distutils.core import setup, Extension 
import numpy 
## define the extension module 
demo_module = Extension('numpy_api_demo', sources=['numpy_api.c'], 
include_dirs=[numpy.get_include()]) 

## run the setup 
setup(ext_modules=[demo_module]) 

由于我们使用特定于 NumPy 的标头，因此我们需要在include_dirs变量中具有numpy.get_include函数。 要运行此安装文件，我们将使用一个熟悉的命令：

python setup.py build_ext -inplace 

前面的命令将在目录中创建一个numpy_api_demo.pyd文件，供我们在 Python 解释器中使用。

为了测试我们的模块，我们将打开一个 Python 解释器测试，并尝试像我们对用 Python 编写的模块所做的一样，从该模块调用这些函数：

>>>import numpy_api_demo as npd 
>>> import numpy as np 
>>>npd.py_square_func(4) 
>>> 16.0 
>>> x = np.arange(0,10,1) 
>>> y = npd.np_square(x) 

总结

在本章中，我们向您介绍了另一种使用 Python 和 NumPy 提供的 C-API 优化或集成 C/C++ 代码的方法。 我们解释了该代码的基本结构以及其他示例代码，开发人员必须编写这些代码才能创建扩展模块。 之后，我们创建了两个函数，这些函数计算出一个数字的平方，并将该平方函数从math.h库映射到一个 Numpy 数组。 这里的目的是使您熟悉如何利用 C/C++ 编写的数字库，以最少的代码重写来创建自己的模块。 编写 C 代码的范围比这里描述的要广泛得多。 但是，我们希望本章使您有信心在需要时利用 C-API。