Python-数学应用-三-

Python 数学应用（三）

原文：zh.annas-archive.org/md5/123a7612a4e578f6816d36f968cfec22

译者：飞龙

协议：CC BY-NC-SA 4.0

第九章：几何问题

本章描述了关于二维几何的几个问题的解决方案。几何是数学的一个分支，涉及点、线和其他图形（形状）的特征，这些图形之间的相互作用以及这些图形的变换。在本章中，我们将重点关注二维图形的特征以及这些对象之间的相互作用。

在 Python 中处理几何对象时，我们必须克服几个问题。最大的障碍是表示问题。大多数几何对象占据二维平面中的一个区域，因此不可能存储区域内的每个点。相反，我们必须找到一种更紧凑的方式来表示可以存储为相对较少的点的区域。例如，我们可以存储沿对象边界的一些点，从而可以重建边界和对象本身。此外，我们将几何问题重新表述为可以使用代表性数据回答的问题。

第二个最大的问题是将纯几何问题转化为可以使用软件理解和解决的形式。这可能相对简单-例如，找到两条直线相交的点是解决矩阵方程的问题-或者可能非常复杂，这取决于所提出的问题类型。解决这些问题的常见技术是使用更简单的对象表示所讨论的图形，并使用每个简单对象解决（希望）更容易的问题。然后，这应该给我们一个关于原始问题的解决方案的想法。

我们将首先向您展示如何可视化二维形状，然后学习如何确定一个点是否包含在另一个图形中。然后，我们将继续查看边缘检测、三角剖分和寻找凸包。最后，我们将通过构造贝塞尔曲线来结束本章。

本章包括以下教程：

• 可视化二维几何形状

• 寻找内部点

• 在图像中查找边缘

• 对平面图形进行三角剖分

• 计算凸包

• 构造贝塞尔曲线

让我们开始吧！

技术要求

对于本章，我们将像往常一样需要numpy包和matplotlib包。我们还需要 Shapely 包和scikit-image包，可以使用您喜欢的软件包管理器（如pip）进行安装：

          python3.8 -m pip install numpy matplotlib shapely scikit-image

该章节的代码可以在 GitHub 存储库的Chapter 08文件夹中找到：github.com/PacktPublishing/Applying-Math-with-Python/tree/master/Chapter%2008。

查看以下视频以查看代码实际操作：bit.ly/3hpeKEF。

可视化二维几何形状

本章的重点是二维几何，因此我们的第一个任务是学习如何可视化二维几何图形。这里提到的一些技术和工具可能适用于三维几何图形，但通常需要更专门的软件包和工具。

几何图形，至少在本书的上下文中，是指边界是一组线和曲线的任何点、线、曲线或封闭区域（包括边界）的集合。简单的例子包括点和线（显然）、矩形、多边形和圆。

在本教程中，我们将学习如何使用 Matplotlib 可视化几何图形。

准备工作

对于本教程，我们需要导入 NumPy 包作为np，导入 Matplotlib pyplot模块作为plt。我们还需要从 Matplotlib patches模块导入Circle类和从 Matplotlib collections模块导入PatchCollection类。可以使用以下命令完成这些操作：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from matplotlib.collections import PatchCollection

我们还需要本章的代码库中的swisscheese-grid-10411.csv数据文件。

如何做...

以下步骤向您展示了如何可视化一个二维几何图形：

1. 首先，我们从本书的代码库中的swisscheese-grid-10411.csv文件加载数据：

data = np.loadtxt("swisscheese-grid-10411.csv")

1. 我们创建一个代表绘图区域的新补丁对象。这将是一个圆（圆盘），其中心在原点，半径为1。我们创建一个新的轴集，并将这个补丁添加到其中：

fig, ax = plt.subplots()
outer = Circle((0.0, 0.0), 1.0, zorder=0, fc="k")
ax.add_patch(outer)

1. 接下来，我们从步骤 1中加载的数据创建一个PatchCollection对象，其中包含了一些其他圆的中心和半径。然后我们将这个PatchCollection添加到步骤 2中创建的轴上：

col = PatchCollection(
    (Circle((x, y), r) for x, y, r in data),
    facecolor="white", zorder=1, linewidth=0.2, 
    ls="-", ec="k"
)
ax.add_collection(col)

1. 最后，我们设置*x-和y-*轴的范围，以便整个图像都能显示出来，然后关闭轴线：

ax.set_xlim((-1.1, 1.1))
ax.set_ylim((-1.1, 1.1))
ax.set_axis_off()

结果图像是一个瑞士奶酪，如下所示：

图 8.1：瑞士奶酪的绘图

它是如何工作的...

这个食谱的关键是Circle和PatchCollection对象，它们代表了 Matplotlib Axes上的绘图区域的区域。在这种情况下，我们创建了一个大的圆形补丁，它位于原点，半径为1，具有黑色的面颜色，并使用zorder=0将其放在其他补丁的后面。这个补丁被添加到Axes对象中使用add_patch方法。

下一步是创建一个对象，它将呈现从 CSV 文件中加载的数据所代表的圆。这些数据包括中心（x，y）和半径r的值，用于表示个别圆的中心和半径（总共 10,411 个）。PatchCollection对象将一系列补丁组合成一个单一的对象，可以添加到Axes对象中。在这里，我们为我们的数据中的每一行添加了一个Circle，然后使用add_collection方法将其添加到Axes对象中。请注意，我们已经将面颜色应用到整个集合，而不是每个单独的Circle成员。我们将面颜色设置为白色（使用facecolor="w"参数），边缘颜色设置为黑色（使用ec="k"），边缘线宽设置为 0.2（使用linewidth=0.2），边缘样式设置为连续线。所有这些放在一起，就得到了我们的图像。

我们在这里创建的图像被称为“瑞士奶酪”。这些最初是由爱丽丝·罗斯在 1938 年在有理逼近理论中使用的；随后它们被重新发现，并且类似的构造自那时以来已经被多次使用。我们使用这个例子是因为它由一个大的个体部分和一个大量的较小的个体部分组成。罗斯的瑞士奶酪是平面上具有正面积但没有拓扑内部的一个集合的例子。（这样一个集合甚至能存在是相当惊人的！）更重要的是，有一些连续函数在这个瑞士奶酪上是不能被有理函数逼近的。这个特性使得类似的构造在均匀 代数理论中非常有用。

Circle类是更一般的Patch类的子类。还有许多其他Patch类，代表不同的平面图形，比如Polygon和PathPatch，它们代表了由路径（曲线或曲线集合）所界定的区域。这些可以用来生成可以在 Matplotlib 图中呈现的复杂补丁。集合可以用来同时应用设置到多个补丁对象上，这在本例中特别有用，因为你有大量的对象都将以相同的样式呈现。

还有更多...

Matplotlib 中有许多不同的补丁类型。在本教程中，我们使用了Circle补丁类，它表示坐标轴上的圆形区域。还有Polygon补丁类，它表示多边形（规则或其他）。还有PatchPath对象，它们是由不一定由直线段组成的曲线包围的区域。这类似于许多矢量图形软件包中可以构造阴影区域的方式。

除了 Matplotlib 中的单个补丁类型外，还有许多集合类型，它们将许多补丁聚集在一起，以便作为单个对象使用。在本教程中，我们使用了PatchCollection类来收集大量的Circle补丁。还有更多专门的补丁集合，可以用来自动生成这些内部补丁，而不是我们自己生成它们。

另请参阅

在数学中，可以在以下传记文章中找到关于瑞士奶酪的更详细的历史：Daepp,U., Gauthier, P., Gorkin, P. and Schmieder, G., 2005. Alice in Switzerland: The life and mathematics of Alice Roth. The Mathematical Intelligencer, 27(1), pp.41-54。

查找内部点

在编程环境中处理二维图形的一个问题是，您不可能存储所有位于图形内的点。相反，我们通常存储表示图形的远少于的点。在大多数情况下，这将是一些点（通过线连接）来描述图形的边界。这在内存方面是有效的，并且可以使用 MatplotlibPatches轻松在屏幕上可视化它们。但是，这种方法使确定点或其他图形是否位于给定图形内变得更加困难。这是许多几何问题中的一个关键问题。

在本教程中，我们将学习如何表示几何图形并确定点是否位于图形内。

做好准备

对于本教程，我们需要将matplotlib包（整体）导入为mpl，将pyplot模块导入为plt：

import matplotlib as mpl
import matplotlib.pyplot as plt

我们还需要从 Shapely 包的geometry模块中导入Point和Polygon对象。Shapely 包包含许多用于表示、操作和分析二维几何图形的例程和对象：

from shapely.geometry import Polygon, Point

如何操作...

以下步骤向您展示如何创建多边形的 Shapely 表示，然后测试点是否位于此多边形内：

1. 创建一个样本多边形进行测试：

polygon = Polygon(
    [(0, 2), (-1, 1), (-0.5, -1), (0.5, -1), (1, 1)],
)

1. 接下来，我们在新图上绘制多边形。首先，我们需要将多边形转换为可以添加到图中的 MatplotlibPolygon补丁：

fig, ax = plt.subplots()
poly_patch = mpl.patches.Polygon(polygon.exterior, ec="k", 
   lw="1", alpha=0.5)
ax.add_patch(poly_patch)
ax.set(xlim=(-1.05, 1.05), ylim=(-1.05, 2.05))
ax.set_axis_off()

1. 现在，我们需要创建两个测试点，其中一个将位于多边形内，另一个将位于多边形外：

p1 = Point(0.0, 0.0)
p2 = Point(-1.0, -0.75)

1. 我们在多边形上方绘制并注释这两个点，以显示它们的位置：

ax.plot(0.0, 0.0, "k*")
ax.annotate("p1", (0.0, 0.0), (0.05, 0.0))
ax.plot(-0.8, -0.75, "k*")
ax.annotate("p2", (-0.8, -0.75), (-0.8 + 0.05, -0.75))

1. 最后，我们使用contains方法测试每个点是否位于多边形内，然后将结果打印到终端：

print("p1 inside polygon?", polygon.contains(p1))
print("p2 inside polygon?", polygon.contains(p2))

结果显示，第一个点p1包含在多边形内，而第二个点p2不包含。这也可以在下图中看到，清楚地显示了一个点包含在阴影多边形内，而另一个点不包含：

图 8.2：多边形区域内外的点

工作原理...

ShapelyPolygon类是多边形的表示，它将其顶点存储为点。外边界包围的区域-存储顶点之间的五条直线-对我们来说是明显的，并且很容易被眼睛识别，但是“在”边界内的概念很难以一种计算机容易理解的方式定义。甚至很难给出关于“在”给定曲线内的含义的正式数学定义。

确定点是否位于简单闭合曲线内有两种主要方法 - 即从同一位置开始并结束且不包含任何自交点的曲线。第一种方法使用数学概念称为绕数，它计算曲线“绕”点的次数，以及射线交叉计数方法，其中我们计算从点到无穷远处的点的射线穿过曲线的次数。幸运的是，我们不需要自己计算这些数字，因为我们可以使用 Shapely 包中的工具来为我们执行这些计算。这就是多边形的contains方法所做的。（在底层，Shapely 使用 GEOS 库执行此计算。）

Shapely 的Polygon类可用于计算与这些平面图形相关的许多数量，包括周长和面积。contains方法用于确定点或一组点是否位于对象表示的多边形内（该类有关多边形的表示存在一些限制）。实际上，您可以使用相同的方法来确定一个多边形是否包含在另一个多边形内，因为正如我们在这个示例中看到的，多边形由一组简单的点表示。

在图像中找到边缘

在图像中找到边缘是减少包含大量噪音和干扰的复杂图像为包含最突出轮廓的非常简单图像的好方法。这可以作为我们分析过程的第一步，例如在图像分类中，或者作为将线轮廓导入计算机图形软件包的过程。

在这个示例中，我们将学习如何使用scikit-image包和 Canny 算法来找到复杂图像中的边缘。

准备工作

对于这个示例，我们需要将 Matplotlib 的pyplot模块导入为plt，从skimage.io模块导入imread例程，以及从skimage.feature模块导入canny例程：

import matplotlib.pyplot as plt
from skimage.io import imread
from skimage.feature import canny

如何做到...

按照以下步骤学习如何使用scikit-image包在图像中找到边缘：

1. 从源文件加载图像数据。这可以在本章的 GitHub 存储库中找到。关键是，我们传入as_gray=True以以灰度加载图像：

image = imread("mandelbrot.png", as_gray=True)

以下是原始图像，供参考。集合本身由白色区域显示，如您所见，边界由较暗的阴影表示，非常复杂：

图 8.3：使用 Python 生成的 Mandelbrot 集合的绘图

1. 接下来，我们使用canny例程，需要从scikit-image包的features模块导入。对于这个图像，sigma值设置为 0.5：

edges = canny(image, sigma=0.5)

1. 最后，我们将edges图像添加到一个新的图中，使用灰度（反转）色图：

fig, ax = plt.subplots()
ax.imshow(edges, cmap="gray_r")
ax.set_axis_off()

已检测到的边缘可以在以下图像中看到。边缘查找算法已经识别出 Mandelbrot 集合边界的大部分可见细节，尽管并不完美（毕竟这只是一个估计）：

图 8.4：使用 scikit-image 包的 Canny 边缘检测算法找到的 Mandelbrot 集合的边缘

它是如何工作的...

scikit-image包提供了各种用于操作和分析从图像中导出的数据的实用程序和类型。正如其名称所示，canny例程使用 Canny 边缘检测算法来找到图像中的边缘。该算法使用图像中的强度梯度来检测边缘，其中梯度较大。它还执行一些过滤以减少它找到的边缘中的噪音。

我们提供的sigma关键字值是应用于图像的高斯平滑的标准偏差，用于计算边缘检测。这有助于我们去除图像中的一些噪音。我们设置的值（0.5）小于默认值（1），但在这种情况下可以提供更好的分辨率。较大的值会遮盖 Mandelbrot 集边界的一些细节。

对平面图形进行三角剖分

正如我们在第三章中看到的，微积分和微分方程，我们经常需要将连续区域分解为更小、更简单的区域。在之前的示例中，我们将实数区间缩小为一系列长度较小的小区间。这个过程通常称为离散化。在本章中，我们正在处理二维图形，因此我们需要这个过程的二维版本。为此，我们将一个二维图形（在这个示例中是一个多边形）分解为一系列更小和更简单的多边形。所有多边形中最简单的是三角形，因此这是二维离散化的一个很好的起点。找到一组"铺砌"几何图形的三角形的过程称为三角剖分。

在这个示例中，我们将学习如何使用 Shapely 包对多边形（带有孔）进行三角剖分。

准备工作

在这个示例中，我们需要将 NumPy 包导入为np，将 Matplotlib 包导入为mpl，并将pyplot模块导入为plt。

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

我们还需要从 Shapely 包中获取以下项目：

from shapely.geometry import Polygon
from shapely.ops import triangulate

如何做...

以下步骤向您展示了如何使用 Shapely 包对带有孔的多边形进行三角剖分：

1. 首先，我们需要创建一个代表我们希望进行三角剖分的图形的Polygon对象：

polygon = Polygon(
    [(2.0, 1.0), (2.0, 1.5), (-4.0, 1.5), (-4.0, 0.5), 
       (-3.0, -1.5), (0.0, -1.5), (1.0, -2.0), (1.0, -0.5), 
         (0.0, -1.0), (-0.5, -1.0), (-0.5, 1.0)],

    holes=[np.array([[-1.5, -0.5], [-1.5, 0.5], [-2.5, 0.5], 
       [-2.5, -0.5]])]
)

1. 现在，我们应该绘制图形，以便了解我们将在其中工作的区域：

fig, ax = plt.subplots()
plt_poly = mpl.patches.Polygon(polygon.exterior, 
   ec="k", lw="1", alpha=0.5, zorder=0)
ax.add_patch(plt_poly)
plt_hole = mpl.patches.Polygon(polygon.interiors[0], 
   ec="k", fc="w")
ax.add_patch(plt_hole)
ax.set(xlim=(-4.05, 2.05), ylim=(-2.05, 1.55))
ax.set_axis_off()

可以在下图中看到这个多边形。正如我们所看到的，这个图形中有一个"孔"，必须仔细考虑：

图 8.5：带有孔的示例多边形

1. 我们使用triangulate例程生成多边形的三角剖分。这个三角剖分包括外部边缘，这是我们在这个示例中不想要的：

triangles = triangulate(polygon)

1. 为了去除位于原始多边形外部的三角形，我们需要使用内置的filter例程，以及contains方法（在本章前面已经看到）：

filtered = filter(lambda p: polygon.contains(p), triangles) 

1. 将三角形绘制在原始多边形上，我们需要将 Shapely 三角形转换为 Matplotlib Patch对象，然后将其存储在PatchCollection中：

patches = map(lambda p: mpl.patches.Polygon(p.exterior), filtered)
col = mpl.collections.PatchCollection(patches, fc="none", ec="k")

1. 最后，我们将三角形补丁的集合添加到之前创建的图形中：

ax.add_collection(col)

在原始多边形上绘制的三角剖分可以在下图中看到。在这里，我们可以看到每个顶点都连接到另外两个顶点，形成了覆盖整个原始多边形的三角形系统：

图 8.6：带有孔的示例多边形的三角剖分

它是如何工作的...

triangulate例程使用一种称为Delaunay 三角剖分的技术将一组点连接到一组三角形中。在这种情况下，这组点是多边形的顶点。Delaunay 方法以这样一种方式找到这些三角形，即没有任何点包含在任何三角形的外接圆内。这是该方法的技术条件，但这意味着三角形被有效地选择，因为它避免了非常长、细的三角形。得到的三角剖分利用了原始多边形中存在的边缘，并连接了一些外部边缘。

为了去除原多边形外的三角形，我们使用内置的filter例程，它通过移除标准函数失败的项目来创建一个新的可迭代对象。这与 Shapely Polygon对象上的contains方法一起使用，以确定每个三角形是否位于原始图形内。正如我们之前提到的，我们需要将这些 Shapely 项目转换为 Matplotlib 补丁，然后才能将它们添加到图中。

还有更多...

三角剖分通常用于将复杂的几何图形简化为一组三角形，这些三角形对于某种计算任务来说要简单得多。然而，它们也有其他用途。三角剖分的一个特别有趣的应用是解决“艺术画廊问题”。这个问题涉及找到必要的“守卫”艺术画廊的最大数量。三角剖分是 Fisk 对艺术画廊定理的简单证明的重要部分，这个定理最初是由 Chvátal 证明的。

假设这个食谱中的多边形是一个艺术画廊的平面图，并且一些守卫需要放置在顶点上。一点工作就会表明，你需要在多边形的顶点处放置三个守卫，整个博物馆才能被覆盖。在下面的图像中，我们绘制了一个可能的布局：

图 8.7：在顶点上放置守卫的艺术画廊问题的一个可能解决方案。

点由点表示，并且它们相应的视野范围被阴影表示。

每个顶点都放置了一个守卫，并且他们的视野范围由相应的阴影区域表示。在这里，你可以看到整个多边形至少被一种颜色覆盖。艺术画廊问题的解决方案——实际上是原问题的一个变体——告诉我们，最多需要四名守卫。

另请参阅

关于艺术画廊问题的更多信息可以在 O'Rourke 的经典著作中找到：ORourke, J. (1987). Art gallery theorems and algorithms. New York: Oxford University Press.

计算凸包

如果图形内的每一对点都可以使用一条直线连接，并且这条直线也包含在图形内，那么几何图形被称为凸。凸体的简单例子包括点、直线、正方形、圆（圆盘）、正多边形等。图 8.5 中显示的几何图形不是凸的，因为孔的对面的点不能通过保持在图形内的直线连接起来。

从某种角度来看，凸图形是简单的，这意味着它们在各种应用中都很有用。一个特别的问题涉及找到包含一组点的最小凸集。这个最小凸集被称为这组点的凸包。

在这个食谱中，我们将学习如何使用 Shapely 包找到一组点的凸包。

准备工作

对于这个食谱，我们需要导入 NumPy 包作为np，导入 Matplotlib 包作为mpl，并导入plt模块：

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

我们还需要从 NumPy 导入默认的随机数生成器。我们可以这样导入：

from numpy.random import default_rng
rng = default_rng(12345)

最后，我们需要从 Shapely 导入MultiPoint类：

from shapely.geometry import MultiPoint

操作方法...

按照以下步骤找到一组随机生成点的凸包：

1. 首先，我们生成一个二维数组的随机数：

raw_points = rng.uniform(-1.0, 1.0, size=(50, 2))

1. 接下来，我们创建一个新图形，并在这个图形上绘制这些原始样本点：

fig, ax = plt.subplots()
ax.plot(raw_points[:, 0], raw_points[:, 1], "k.")
ax.set_axis_off()

这些随机生成的点可以在下图中看到。这些点大致分布在一个正方形区域内：

图 8.8：平面上的一组点

1. 接下来，我们构建一个MultiPoint对象，收集所有这些点并将它们放入一个单一对象中：

points = MultiPoint(raw_points)

1. 现在，我们使用convex_hull属性获取这个MultiPoint对象的凸包：

convex_hull = points.convex_hull

1. 然后，我们创建一个 MatplotlibPolygon补丁，可以在我们的图中绘制，以显示找到的凸包的结果：

patch = mpl.patches.Polygon(convex_hull.exterior, alpha=0.5,
   ec="k", lw=1.2)

1. 最后，我们将Polygon补丁添加到图中，以显示凸包：

ax.add_patch(patch)

随机生成的点的凸包可以在下图中看到：

图 8.9：平面上一组点的凸包

工作原理...

Shapely 包是围绕 GEOS 库的 Python 包装器，用于几何分析。Shapely 几何对象的convex_hull属性调用 GEOS 库中的凸包计算例程，从而产生一个新的 Shapely 对象。从这个教程中，我们可以看到一组点的凸包是一个多边形，其顶点是离“中心”最远的点。

构造贝塞尔曲线

贝塞尔曲线，或B 样条，是一族曲线，在矢量图形中非常有用-例如，它们通常用于高质量的字体包中。这是因为它们由少量点定义，然后可以用来廉价地计算沿曲线的大量点。这允许根据用户的需求来缩放细节。

在本教程中，我们将学习如何创建一个表示贝塞尔曲线的简单类，并计算沿其路径的若干点。

准备工作

在本教程中，我们将使用导入为np的 NumPy 包，导入为plt的 Matplotlib pyplot模块，以及 Python 标准库math模块中导入为binom的comb例程：

from math import comb as binom
import matplotlib.pyplot as plt
import numpy as np

如何做...

按照以下步骤定义一个表示贝塞尔曲线的类，该类可用于计算沿曲线的点：

1. 第一步是设置基本类。我们需要为实例属性提供控制点（节点）和一些相关的数字：

class Bezier:
    def __init__(self, *points):
        self.points = points
        self.nodes = n = len(points) - 1
        self.degree = l = points[0].size

1. 仍然在__init__方法中，我们生成贝塞尔曲线的系数，并将它们存储在实例属性的列表中：

        self.coeffs = [binom(n, i)*p.reshape((l, 1)) for i, 
           p in enumerate(points)]

1. 接下来，我们定义一个__call__方法，使类可调用。我们将实例中的节点数加载到本地变量中，以便清晰明了：

    def __call__(self, t):
        n = self.nodes

1. 接下来，我们重新整理输入数组，使其包含单行：

        t = t.reshape((1, t.size))

1. 现在，我们使用实例的coeffs属性中的每个系数生成值数组的列表：

        vals = [c @ (t**i)*(1-t)**(n-i) for i, 
           c in enumerate(self.coeffs)]

1. 最后，我们对步骤 5中构造的所有数组进行求和，并返回结果数组：

       return np.sum(vals, axis=0)

1. 现在，我们将通过一个示例来测试我们的类。我们将为此示例定义四个控制点：

p1 = np.array([0.0, 0.0])
p2 = np.array([0.0, 1.0])
p3 = np.array([1.0, 1.0])
p4 = np.array([1.0, 3.0])

1. 接下来，我们为绘图设置一个新的图形，并用虚线连接线绘制控制点：

fig, ax = plt.subplots()
ax.plot([0.0, 0.0, 1.0, 1.0], [0.0, 1.0, 1.0, 3.0], "*--k")
ax.set(xlabel="x", ylabel="y", title="Bezier curve with 
    4 nodes, degree 3")

1. 然后，我们使用步骤 7中定义的四个点创建我们的Bezier类的新实例：

b_curve = Bezier(p1, p2, p3, p4)

1. 现在，我们可以使用linspace创建 0 到 1 之间等间距点的数组，并计算沿着贝塞尔曲线的点：

t = np.linspace(0, 1)
v = b_curve(t)

1. 最后，我们在之前绘制的控制点上绘制这条曲线：

ax.plot(v[0,:], v[1, :])

我们绘制的贝塞尔曲线可以在下图中看到。正如你所看到的，曲线从第一个点（0, 0）开始，结束于最终点（1, 3）：

图 8.10：使用四个节点构造的三次贝塞尔曲线

工作原理...

贝塞尔曲线由一系列控制点描述，我们以递归方式构造曲线。一个点的贝塞尔曲线是一个保持在该点的常数曲线。具有两个控制点的贝塞尔曲线是这两个点之间的线段：

当我们添加第三个控制点时，我们取对应点之间的线段，这些点是由一个较少点构成的贝塞尔曲线的曲线。这意味着我们使用以下公式构造具有三个控制点的贝塞尔曲线：

这种构造可以在下图中看到：

图 8.11：使用递归定义构造二次贝塞尔曲线。黑色虚线显示了两条线性贝塞尔曲线。

这种构造方式继续定义了任意数量控制点上的贝塞尔曲线。幸运的是，在实践中我们不需要使用这种递归定义，因为我们可以将公式展开成曲线的单一公式，即以下公式：

在这里，p[i]元素是控制点，t是一个参数，而

是二项式系数。请记住，t参数是生成曲线点的变化量。我们可以分离前述求和中涉及t的项和不涉及t的项。这定义了我们在步骤 2中定义的系数，每个系数由以下代码片段给出：

binom(n, i)*p.reshape((l, 1))

我们在这一步中对每个点p进行了 reshape，以确保它被排列为列向量。这意味着每个系数都是一个列向量（作为 NumPy 数组），由二项式系数缩放的控制点组成。

现在，我们需要指定如何在不同的t值上评估贝塞尔曲线。这就是我们利用 NumPy 包中的高性能数组操作的地方。在形成系数时，我们将控制点 reshape 为列向量。在步骤 4中，我们将输入t值 reshape 为行向量。这意味着我们可以使用矩阵乘法运算符将每个系数乘以相应的（标量）值，具体取决于输入的t。这就是步骤 5中列表推导式中发生的情况。在下一行中，我们将l×1数组乘以1×N数组，得到一个l×N数组：

c @ (t**i)*(1-t)**(n-i)

我们为每个系数都得到一个这样的数组。然后，我们可以使用np.sum例程来对这些l×N数组中的每一个进行求和，以得到贝塞尔曲线上的值。在本示例中，输出数组的顶行包含曲线的x值，底行包含曲线的y值。在指定axis=0关键字参数时，我们必须小心确保sum例程对我们创建的列表进行求和，而不是对该列表包含的数组进行求和。

我们定义的类是使用贝塞尔曲线的控制点进行初始化的，然后用于生成系数。曲线值的实际计算是使用 NumPy 完成的，因此这种实现应该具有相对良好的性能。一旦创建了这个类的特定实例，它的功能就非常像一个函数，正如你所期望的那样。但是，这里没有进行类型检查，所以我们只能用 NumPy 数组作为参数来调用这个“函数”。

还有更多...

贝塞尔曲线是使用迭代构造定义的，其中具有n个点的曲线是使用连接由第一个和最后一个n-1点定义的曲线来定义的。使用这种构造跟踪每个控制点的系数将很快导致我们用来定义前述曲线的方程。这种构造还导致贝塞尔曲线的有趣和有用的几何特性。

正如我们在这个配方的介绍中提到的，贝塞尔曲线出现在许多涉及矢量图形的应用程序中，比如字体。它们也出现在许多常见的矢量图形软件包中。在这些软件包中，通常会看到二次贝塞尔曲线，它们由三个点的集合定义。然而，你也可以通过提供两个端点以及这些点上的梯度线来定义一个二次贝塞尔曲线。这在图形软件包中更常见。生成的贝塞尔曲线将沿着梯度线离开每个端点，并在这些点之间平滑地连接曲线。

我们在这里构建的实现对于小型应用程序来说性能相对较好，但对于涉及在大量t值上渲染具有大量控制点的曲线的应用程序来说是不够的。对于这一点，最好使用一个用编译语言编写的低级软件包。例如，bezier Python 软件包使用编译的 Fortran 后端进行计算，并提供比我们在这里定义的类更丰富的接口。

当然，贝塞尔曲线可以自然地扩展到更高的维度。结果是一个贝塞尔曲面，使它们成为非常有用的通用工具，用于高质量、可伸缩的图形。

进一步阅读

• 计算几何中一些常见算法的描述可以在以下书籍中找到：Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., 2007. Numerical recipes: the art of scientific computing*. 3rd ed. Cambridge: Cambridge University Press*。

• 有关计算几何中一些问题和技术的更详细描述，请查阅以下书籍：O'Rourke, J., 1994. Computational geometry in C*. Cambridge: Cambridge University Press*。

第十章：寻找最优解

在本章中，我们将讨论寻找给定情况下最佳结果的各种方法。这被称为优化，通常涉及最小化或最大化目标函数。目标函数是一个接受多个参数作为参数并返回代表给定参数选择的成本或回报的单个标量值的函数。关于最小化和最大化函数的问题实际上是相互等价的，因此我们只会在本章讨论最小化目标函数。最小化函数f(x)等同于最大化函数*-f*(x)。在我们讨论第一个配方时将提供更多细节。

我们可以利用的算法来最小化给定函数取决于函数的性质。例如，包含一个或多个变量的简单线性函数与具有许多变量的非线性函数相比，可用的算法不同。线性函数的最小化属于线性规划范畴，这是一个发展完善的理论。对于非线性函数，我们通常利用函数的梯度（导数）来寻找最小点。我们将讨论几种不同类型函数的最小化方法。

寻找单变量函数的极小值和极大值特别简单，如果函数的导数已知，可以轻松完成。如果不知道导数，则适用于适当配方的方法。最小化非线性函数配方中的注释提供了一些额外细节。

我们还将提供一个非常简短的介绍博弈论。广义上讲，这是一个围绕决策制定的理论，并在经济学等学科中具有广泛的影响。特别是，我们将讨论如何在 Python 中将简单的双人游戏表示为对象，计算与某些选择相关的回报，并计算这些游戏的纳什均衡。

我们将首先看如何最小化包含一个或多个变量的线性和非线性函数。然后，我们将继续研究梯度下降方法和使用最小二乘法进行曲线拟合。最后，我们将通过分析双人游戏和纳什均衡来结束本章。

在本章中，我们将涵盖以下配方：

• 最小化简单线性函数

• 最小化非线性函数

• 使用梯度下降法进行优化

• 使用最小二乘法拟合数据的曲线

• 分析简单的双人游戏

• 计算纳什均衡

让我们开始吧！

技术要求

在本章中，我们将像往常一样需要 NumPy 包、SciPy 包和 Matplotlib 包。我们还将需要 Nashpy 包用于最后两个配方。这些包可以使用您喜欢的包管理器（如pip）进行安装：

          python3.8 -m pip install numpy scipy matplotlib nashpy

本章的代码可以在 GitHub 存储库的Chapter 09文件夹中找到，网址为github.com/PacktPublishing/Applying-Math-with-Python/tree/master/Chapter%2009。

查看以下视频以查看代码的实际操作：bit.ly/2BjzwGo。

最小化简单线性函数

在优化中我们面临的最基本问题是找到函数取得最小值的参数。通常，这个问题受到参数可能值的一些限制的约束，这增加了问题的复杂性。显然，如果我们要最小化的函数也很复杂，那么这个问题的复杂性会进一步增加。因此，我们必须首先考虑线性函数，它们的形式如下：

为了解决这些问题，我们需要将约束转化为计算机可用的形式。在这种情况下，我们通常将它们转化为线性代数问题（矩阵和向量）。一旦完成了这一步，我们就可以使用 NumPy 和 SciPy 中的线性代数包中的工具来找到我们所寻求的参数。幸运的是，由于这类问题经常发生，SciPy 有处理这种转化和随后求解的例程。

在这个配方中，我们将使用 SciPy optimize模块的例程来解决以下受限线性最小化问题：

这将受到以下条件的约束：

准备工作

对于这个配方，我们需要在别名np下导入 NumPy 包，以plt的名称导入 Matplotlib pyplot模块，以及 SciPy optimize模块。我们还需要从mpl_toolkits.mplot3d导入Axes3D类，以使 3D 绘图可用：

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

如何做...

按照以下步骤使用 SciPy 解决受限线性最小化问题：

1. 将系统设置为 SciPy 可以识别的形式：

A = np.array([
    [2, 1],   # 2*x0 + x1 <= 6
    [-1, -1]  # -x0 - x1 <= -4
])
b = np.array([6, -4])
x0_bounds = (-3, 14) # -3 <= x0 <= 14
x1_bounds = (2, 12)  # 2 <= x1 <= 12
c = np.array([1, 5])

1. 接下来，我们需要定义一个评估线性函数在x值处的例程，这是一个向量（NumPy 数组）：

def func(x):
    return np.tensordot(c, x, axes=1)

1. 然后，我们创建一个新的图并添加一组3d轴，我们可以在上面绘制函数：

fig = plt.figure()
ax = fig.add_subplot(projection="3d")
ax.set(xlabel="x0", ylabel="x1", zlabel="func")
ax.set_title("Values in Feasible region")

1. 接下来，我们创建一个覆盖问题区域的值网格，并在该区域上绘制函数的值：

X0 = np.linspace(*x0_bounds)
X1 = np.linspace(*x1_bounds)
x0, x1 = np.meshgrid(X0, X1)
z = func([x0, x1])
ax.plot_surface(x0, x1, z, alpha=0.3)

1. 现在，我们在函数值平面上绘制与临界线2*x0 + x1 == 6对应的线，并在我们的图上绘制落在范围内的值：

Y = (b[0] - A[0, 0]*X0) / A[0, 1]
I = np.logical_and(Y >= x1_bounds[0], Y <= x1_bounds[1])
ax.plot(X0[I], Y[I], func([X0[I], Y[I]]), "r", lw=1.5)

1. 我们重复这个绘图步骤，对第二条临界线x0 + x1 == -4：

Y = (b[1] - A[1, 0]*X0) / A[1, 1]
I = np.logical_and(Y >= x1_bounds[0], Y <= x1_bounds[1])
ax.plot(X0[I], Y[I], func([X0[I], Y[I]]), "r", lw=1.5)

1. 接下来，我们着色位于两条临界线之间的区域，这对应于最小化问题的可行区域：

B = np.tensordot(A, np.array([x0, x1]), axes=1)
II = np.logical_and(B[0, ...] <= b[0], B[1, ...] <= b[1]) 
ax.plot_trisurf(x0[II], x1[II], z[II], color="b", alpha=0.5)

函数值在可行区域上的图可以在以下图片中看到：

图 9.1：突出显示了可行区域的线性函数值

正如我们所看到的，位于这个着色区域内的最小值发生在两条临界线的交点处。

1. 接下来，我们使用linprog来解决带有我们在步骤 1中创建的边界的受限最小化问题。我们在终端中打印出结果对象：

res = optimize.linprog(c, A_ub=A, b_ub=b, bounds=
    (x0_bounds, x1_bounds))
print(res)

1. 最后，我们在可行区域上绘制最小函数值：

ax.plot([res.x[0]], [res.x[1]], [res.fun], "k*")

更新后的图可以在以下图片中看到：

图 9.2：在可行区域上绘制的最小值

在这里，我们可以看到linprog例程确实发现了最小值在两条临界线的交点处。

工作原理...

受限线性最小化问题在经济情况中很常见，您尝试在保持参数的其他方面的同时最小化成本。实际上，优化理论中的许多术语都反映了这一事实。解决这些问题的一个非常简单的算法称为单纯形法，它使用一系列数组操作来找到最小解。从几何上讲，这些操作代表着改变单纯形的不同顶点（我们在这里不定义），正是这一点赋予了算法其名称。

在我们继续之前，我们将简要概述单纯形法用于解决受限线性优化问题的过程。我们所面临的问题不是一个矩阵方程问题，而是一个矩阵不等式问题。我们可以通过引入松弛变量来解决这个问题，将不等式转化为等式。例如，通过引入松弛变量s[1]，可以将第一个约束不等式重写如下：

只要*s[1]不是负数，这就满足了所需的不等式。第二个约束不等式是大于或等于类型的不等式，我们必须首先将其更改为小于或等于类型。我们通过将所有项乘以-1 来实现这一点。这给了我们在配方中定义的矩阵A的第二行。引入第二个松弛变量s[2]*后，我们得到了第二个方程：

从中，我们可以构建一个矩阵，其列包含两个参数变量*x[1]和x[2]以及两个松弛变量s[1]和s[2]的系数。该矩阵的行代表两个边界方程和目标函数。现在可以使用对该矩阵进行初等行操作来解决这个方程组，以获得最小化目标函数的x[1]和x[2]*的值。由于解决矩阵方程很容易且快速，这意味着我们可以快速高效地最小化线性函数。

幸运的是，我们不需要记住如何将不等式系统化简为线性方程组，因为诸如linprog之类的例程会为我们完成这一工作。我们只需将边界不等式提供为矩阵和向量对，包括每个系数，以及定义目标函数的单独向量。linprog例程负责制定和解决最小化问题。

实际上，单纯形法并不是linprog例程用于最小化函数的算法。相反，linprog使用内点算法，这更有效率。（可以通过提供method关键字参数并使用适当的方法名称将方法设置为simplex或revised-simplex。在打印的输出结果中，我们可以看到只用了五次迭代就达到了解决方案。）该例程返回的结果对象包含最小值发生的参数值存储在x属性中，该最小值存储在fun属性中，以及有关解决过程的各种其他信息。如果方法失败，那么status属性将包含一个数值代码，描述方法失败的原因。

在本文的步骤 2中，我们创建了一个代表此问题的目标函数的函数。该函数以单个数组作为输入，该数组包含应在其上评估函数的参数空间值。在这里，我们使用了 NumPy 的tensordot例程（带有axes=1）来评估系数向量c与每个输入x的点积。在这里我们必须非常小心，因为我们传递给函数的值在后续步骤中将是一个 2×50×50 的数组。普通的矩阵乘法（np.dot）在这种情况下不会给出我们所期望的 50×50 数组输出。

在步骤 5和6中，我们计算了临界线上的点，这些点具有以下方程：

然后我们计算了相应的z值，以便绘制在由目标函数定义的平面上的线。我们还需要“修剪”这些值，以便只包括在问题中指定范围内的值。

还有更多...

本文介绍了受约束的最小化问题以及如何使用 SciPy 解决它。然而，相同的方法也可以用于解决受约束的最大化问题。这是因为最大化和最小化在某种意义上是对偶的，即最大化函数f(x)等同于最小化函数*-f*(x)，然后取其负值。事实上，我们在本文中使用了这一事实，将第二个约束不等式从≥改为≤。

在这个示例中，我们解决了一个只有两个参数变量的问题，但是相同的方法将适用于涉及两个以上这样的变量的问题（除了绘图步骤）。我们只需要向每个数组添加更多的行和列，以考虑这增加的变量数量 - 这包括提供给例程的边界元组。在处理非常大量的变量时，例程也可以与稀疏矩阵一起使用，以获得额外的效率。

linprog例程的名称来自线性规划，用于描述这种类型的问题 - 找到满足一些矩阵不等式的x的值，受其他条件的限制。由于与矩阵理论和线性代数的理论有非常紧密的联系，因此在线性规划问题中有许多非常快速和高效的技术，这些技术在非线性情境中是不可用的。

最小化非线性函数

在上一个示例中，我们看到了如何最小化一个非常简单的线性函数。不幸的是，大多数函数都不是线性的，通常也没有我们希望的良好性质。对于这些非线性函数，我们不能使用为线性问题开发的快速算法，因此我们需要设计可以在这些更一般情况下使用的新方法。我们将使用的算法称为 Nelder-Mead 算法，这是一种健壮且通用的方法，用于找到函数的最小值，并且不依赖于函数的梯度。

在这个示例中，我们将学习如何使用 Nelder-Mead 单纯形法来最小化包含两个变量的非线性函数。

准备工作

在这个示例中，我们将使用导入为np的 NumPy 包，导入为plt的 Matplotlib pyplot模块，从mpl_toolkits.mplot3d导入的Axes3D类以启用 3D 绘图，以及 SciPy optimize模块：

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

如何做...

以下步骤向您展示如何使用 Nelder-Mead 单纯形法找到一般非线性目标函数的最小值：

1. 定义我们将最小化的目标函数：

def func(x):
    return ((x[0] - 0.5)**2 + (x[1] + 0.5)**2)*
       np.cos(0.5*x[0]*x[1])

1. 接下来，创建一个值网格，我们可以在上面绘制我们的目标函数：

x_r = np.linspace(-1, 1)
y_r = np.linspace(-2, 2)
x, y = np.meshgrid(x_r, y_r)

1. 现在，我们在这个点网格上评估函数：

z = func([x, y])

1. 接下来，我们创建一个带有3d轴对象的新图，并设置轴标签和标题：

fig = plt.figure(tight_layout=True)
ax = fig.add_subplot(projection="3d")
ax.tick_params(axis="both", which="major", labelsize=9)
ax.set(xlabel="x", ylabel="y", zlabel="z")
ax.set_title("Objective function")

1. 现在，我们可以在我们刚刚创建的轴上将目标函数绘制为表面：

ax.plot_surface(x, y, z, alpha=0.7)

1. 我们选择一个初始点，我们的最小化例程将从该点开始迭代，并在表面上绘制这个点：

x0 = np.array([-0.5, 1.0])
ax.plot([x0[0]], [x0[1]], func(x0), "r*")

可以在以下图像中看到目标函数表面的绘图，以及初始点。在这里，我们可以看到最小值似乎出现在 x 轴上约 0.5，y 轴上约-0.5 的位置：

图 9.3：具有起始值的非线性目标函数

1. 现在，我们使用optimize包中的minimize例程来找到最小值，并打印它产生的result对象：

result = optimize.minimize(func, x0, tol=1e-6, method=
    "Nelder-Mead")
print(result)

1. 最后，我们在目标函数表面上绘制minimize例程找到的最小值：

ax.plot([result.x[0]], [result.x[1]], [result.fun], "r*")

包括minimize例程找到的最小点的目标函数的更新绘图可以在以下图像中看到：

图 9.4：具有起始点和最小点的目标函数

它是如何工作的...

Nelder-Mead 单纯形法-不要与线性优化问题的单纯形法混淆-是一种简单的算法，用于找到非线性函数的最小值，即使目标函数没有已知的导数也可以工作。（这不适用于此示例中的函数；使用基于梯度的方法的唯一收益是收敛速度。）该方法通过比较单纯形的顶点处的目标函数值来工作，在二维空间中是一个三角形。具有最大函数值的顶点通过相反的边被“反射”，并执行适当的扩展或收缩，实际上将单纯形“下坡”移动。

SciPyoptimize模块中的minimize例程是许多非线性函数最小化算法的入口点。在这个示例中，我们使用了 Nelder-Mead 单纯形算法，但还有许多其他可用的算法。其中许多算法需要对函数的梯度有所了解，该梯度可能会被算法自动计算。可以通过向method关键字参数提供适当的名称来使用该算法。

“最小化”例程返回的result对象包含有关求解器找到的解决方案的大量信息-或者如果发生错误，则未找到-。特别是，计算出的最小值发生的期望参数存储在结果的x属性中，而函数值存储在fun属性中。

“最小化”例程需要函数和x0的起始值。在这个示例中，我们还提供了一个容差值，最小值应该使用tol关键字参数计算。更改此值将修改计算解的准确度。

还有更多...

Nelder-Mead 算法是“无梯度”最小化算法的一个例子，因为它不需要任何关于目标函数的梯度（导数）的信息。有几种这样的算法，通常涉及在指定点评估目标函数，然后使用这些信息朝向最小值移动。一般来说，无梯度方法的收敛速度比梯度下降模型慢。但是，它们可以用于几乎任何目标函数，即使很难精确计算梯度或通过近似手段计算。

通常，优化单变量函数比多维情况更容易，并且在 SciPyoptimize库中有其专用函数。minimize_scalar例程对单变量函数执行最小化，并且在这种情况下应该使用而不是minimize。

在优化中使用梯度下降方法

在上一个示例中，我们使用 Nelder-Mead 单纯形算法最小化包含两个变量的非线性函数。这是一种相当健壮的方法，即使对目标函数了解甚少也可以工作。然而，在许多情况下，我们对目标函数了解更多，这一事实使我们能够设计更快和更有效的最小化函数的算法。我们可以通过利用函数的梯度等属性来做到这一点。

多于一个变量的函数的梯度描述了函数在各个分量方向上的变化率。这是函数对每个变量的偏导数的向量。从这个梯度向量中，我们可以推断出函数在哪个方向上增长最快，反之亦然，从任意给定位置开始，函数在哪个方向上下降最快。这为梯度下降方法最小化函数提供了基础。算法非常简单：给定一个起始位置x，我们计算在这个x处的梯度以及梯度最快下降的相应方向，然后沿着那个方向迈出一小步。经过几次迭代，这将从起始位置移动到函数的最小值。

在这个示例中，我们将学习如何实现基于最陡下降算法的算法，以在有界区域内最小化目标函数。

准备工作

对于这个示例，我们需要导入 NumPy 包作为np，导入 Matplotlib 的pyplot模块作为plt，并从mpl_toolkits.mplot3d导入Axes3D对象：

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

如何做...

在接下来的步骤中，我们将实现一个简单的梯度下降方法，以最小化具有已知梯度函数的目标函数（实际上我们将使用一个生成器函数，以便我们可以看到方法的工作方式）：

1. 我们将首先定义一个descend例程，它将执行我们的算法。函数声明如下：

def descend(func, x0, grad, bounds, tol=1e-8, max_iter=100):

1. 接下来，我们需要实现这个例程。我们首先定义将在方法运行时保存迭代值的变量：

xn = x0
xnm1 = np.inf
grad_xn = grad(x0)

1. 然后，我们开始我们的循环，这将运行迭代。我们立即检查是否在继续之前正在取得有意义的进展：

for i in range(max_iter):
    if np.linalg.norm(xn - xnm1) < tol:
        break

1. 方向是梯度向量的负。我们计算一次并将其存储在direction变量中：

direction = -grad_xn

1. 现在，我们更新先前和当前的值，分别为xnm1和xn，准备进行下一次迭代。这结束了descend例程的代码：

xnm1 = xn
xn = xn + 0.2*direction

1. 现在，我们可以计算当前值的梯度并产生所有适当的值：

grad_xn = grad(xn)
yield i, xn, func(xn), grad_xn

这结束了descend例程的定义。

1. 我们现在可以定义一个要最小化的样本目标函数：

def func(x):
    return ((x[0] - 0.5)**2 + (x[1] + 0.5)**2)*np.cos(0.5*x[0]*x[1])

1. 接下来，我们创建一个网格，我们将评估并绘制目标函数：

x_r = np.linspace(-1, 1)
y_r = np.linspace(-2, 2)
x, y = np.meshgrid(x_r, y_r)

1. 一旦网格创建完成，我们可以评估我们的函数并将结果存储在z变量中：

z = func([x, y])

1. 接下来，我们创建目标函数的三维表面图：

surf_fig = plt.figure(tight_layout=True)
surf_ax = surf_fig.add_subplot(projection="3d")
surf_ax.tick_params(axis="both", which="major", labelsize=9)
surf_ax.set(xlabel="x", ylabel="y", zlabel="z")
surf_ax.set_title("Objective function")
surf_ax.plot_surface(x, y, z, alpha=0.7)

1. 在我们开始最小化过程之前，我们需要定义一个初始点x0。我们在前一步中创建的目标函数图上绘制这一点：

x0 = np.array([-0.8, 1.3])
surf_ax.plot([x0[0]], [x0[1]], func(x0), "r*")

目标函数的表面图，以及初始值，可以在以下图像中看到：

图 9.5：具有初始位置的目标函数表面

1. 我们的descend例程需要一个评估目标函数梯度的函数，因此我们将定义一个：

def grad(x):
    c1 = x[0]**2 - x[0] + x[1]**2 + x[1] + 0.5
    cos_t = np.cos(0.5*x[0]*x[1])
    sin_t = np.sin(0.5*x[0]*x[1])
    return np.array([
        (2*x[0]-1)*cos_t - 0.5*x[1]*c1*sin_t,
        (2*x[1]+1)*cos_t - 0.5*x[0]*c1*sin_t
    ])

1. 我们将在轮廓图上绘制迭代，因此我们设置如下：

cont_fig, cont_ax = plt.subplots()
cont_ax.set(xlabel="x", ylabel="y")
cont_ax.set_title("Contour plot with iterates")
cont_ax.contour(x, y, z, levels=30)

1. 现在，我们创建一个变量，将x和y方向的边界作为元组的元组保存。这些是步骤 10中linspace调用中的相同边界：

bounds = ((-1, 1), (-2, 2))

1. 我们现在可以使用for循环驱动descend生成器来产生每个迭代，并将步骤添加到轮廓图中：

xnm1 = x0
for i, xn, fxn, grad_xn in descend(func, x0, grad, bounds):
    cont_ax.plot([xnm1[0], xn[0]], [xnm1[1], xn[1]], "k*--")
    xnm1, grad_xnm1 = xn, grad_xn

1. 循环完成后，我们将最终值打印到终端：

print(f"iterations={i}")
print(f"min val at {xn}")
print(f"min func value = {fxn}")

前面打印语句的输出如下：

iterations=37
min val at [ 0.49999999 -0.49999999]
min func value = 2.1287163880894953e-16

在这里，我们可以看到我们的例程使用了 37 次迭代来找到大约在(0.5，-0.5)处的最小值，这是正确的。

可以在以下图像中看到带有其迭代的轮廓图：

图 9.6：梯度下降迭代到最小值的目标函数轮廓图

在这里，我们可以看到每次迭代的方向 - 由虚线表示 - 是目标函数下降最快的方向。最终迭代位于目标函数“碗”的中心，这是最小值出现的地方。

它是如何工作的...

这个配方的核心是descend例程。在这个例程中定义的过程是梯度下降法的一个非常简单的实现。在给定点计算梯度由grad参数处理，然后用direction = -grad推断迭代的旅行方向。我们将这个方向乘以一个固定的比例因子（有时称为学习率）0.2 的值来获得缩放步骤，然后通过添加0.2*direction到当前位置来进行这一步。

这个配方的解决方案需要 37 次迭代才能收敛，这比最小化非线性函数配方中的 Nelder-Mead 单纯形算法需要 58 次迭代要好一点。（这并不是一个完美的比较，因为我们改变了这个配方的起始位置。）这种性能在很大程度上取决于我们选择的步长。在这种情况下，我们将最大步长固定为方向向量的 0.2 倍。这使得算法简单，但并不特别高效。

在这个配方中，我们选择将算法实现为生成器函数，以便我们可以看到每一步的输出，并在迭代中绘制在等高线图上。在实践中，我们可能不想这样做，而是在迭代完成后返回计算出的最小值。为此，我们可以简单地删除yield语句，并在函数的最后，即主函数的缩进位置（即不在循环内），用return xn替换它。如果你想防止不收敛，你可以使用for循环的else特性来捕捉循环完成的情况，因为它已经到达了迭代器的末尾而没有触发break关键字。这个else块可以引发一个异常，以表明算法未能稳定到一个解决方案。在这个配方中，我们用来结束迭代的条件并不能保证该方法已经达到了最小值，但这通常是这种情况。

还有更多...

在实践中，你通常不会为自己实现梯度下降算法，而是使用来自库（如 SciPy optimize模块）的通用例程。我们可以使用与前一个配方中使用的相同的minimize例程来执行多种不同算法的最小化，包括几种梯度下降算法。这些实现可能比这样的自定义实现具有更高的性能和更强大的鲁棒性。

我们在这个配方中使用的梯度下降方法是一个非常天真的实现，可以通过允许例程在每一步选择步长来大大改进。（允许选择自己步长的方法有时被称为自适应方法。）这种改进的困难部分是选择在这个方向上采取的步长大小。为此，我们需要考虑单变量函数，它由以下方程给出：

在这里，x*[n]表示当前点，d[n]表示当前方向，t是一个参数。为了简单起见，我们可以使用 SciPy optimize模块中的minimize_scalar最小化例程来处理标量值函数。不幸的是，这并不像简单地传入这个辅助函数并找到最小值那样简单。我们必须限定t*的可能值，以便计算出的最小化点x[n]+ td[n]位于我们感兴趣的区域内。

要理解我们如何限制t的值，我们必须首先从几何上看这个构造。我们引入的辅助函数在给定方向上沿着一条直线评估目标函数。我们可以将其想象为通过当前x*[n]点在d[n]方向上穿过的表面的单个横截面。算法的下一步是找到最小化沿着这条直线的目标函数值的步长t*，这是一个标量函数，要最小化它要容易得多。然后边界应该是t值的范围，在这个范围内，这条直线位于由x和y边界值定义的矩形内。我们确定这条直线穿过这些x和y边界线的四个值，其中两个将是负值，另外两个将是正值。（这是因为当前点必须位于矩形内。）我们取两个正值中的最小值和两个负值中的最大值，并将这些边界传递给标量最小化例程。这是通过以下代码实现的：

alphas = np.array([
        (bounds[0][0] - xn[0]) / direction[0], # x lower
        (bounds[1][0] - xn[1]) / direction[1], # y lower
        (bounds[0][1] - xn[0]) / direction[0], # x upper
        (bounds[1][1] - xn[1]) / direction[1] # y upper
])

alpha_max = alphas[alphas >= 0].min()
alpha_min = alphas[alphas < 0].max()
result = minimize_scalar(lambda t: func(xn + t*direction), 
        method="bounded", bounds=(alpha_min, alpha_max))
amount = result.x

一旦选择了步长，唯一剩下的步骤就是更新当前的xn值，如下所示：

xn = xn + amount * direction

使用这种自适应步长会增加例程的复杂性，但性能大大提高。使用这个修改后的例程，该方法在只进行了三次迭代后就收敛了，远少于本食谱中使用的朴素代码（37 次迭代）或上一个食谱中使用的 Nelder-Mead 单纯形算法（58 次迭代）。通过在梯度函数中提供更多信息，减少了迭代次数，这正是我们预期的。

我们创建了一个函数，返回给定点处函数的梯度。我们在开始之前手动计算了这个梯度，这并不总是容易或甚至可能的。相反，更常见的是用数值计算的梯度替换这里使用的“解析”梯度，这个数值梯度是使用有限差分或类似算法估计的。这对性能和准确性有影响，就像所有近似一样，但鉴于梯度下降方法提供的收敛速度的提高，这些问题通常是次要的。

梯度下降类型的算法在机器学习应用中特别受欢迎。大多数流行的 Python 机器学习库，包括 PyTorch、TensorFlow 和 Theano，都提供了用于自动计算数据数组梯度的实用工具。这使得梯度下降方法可以在后台使用以提高性能。

梯度下降方法的一种流行变体是随机梯度下降，其中梯度是通过随机抽样来估计的，而不是使用整个数据集。这可以显著减少方法的计算负担，但会以更慢的收敛速度为代价，特别是对于高维问题，比如在机器学习应用中常见的问题。随机梯度下降方法通常与反向传播结合，构成了机器学习应用中训练人工神经网络的基础。

基本随机梯度下降算法有几种扩展。例如，动量算法将先前的增量纳入到下一个增量的计算中。另一个例子是自适应梯度算法，它纳入每个参数的学习率，以改善涉及大量稀疏参数的问题的收敛速度。

使用最小二乘法拟合数据曲线

最小二乘法是一种强大的技术，用于从相对较小的潜在函数族中找到最能描述特定数据集的函数。这种技术在统计学中特别常见。例如，最小二乘法用于线性回归问题-在这里，潜在函数族是所有线性函数的集合。通常，我们尝试拟合的这个函数族具有相对较少的参数，可以调整以解决问题。

最小二乘法的思想相对简单。对于每个数据点，我们计算残差的平方-点的值与给定函数的期望值之间的差异-并尝试使这些残差的平方和尽可能小（因此最小二乘法）。

在这个配方中，我们将学习如何使用最小二乘法来拟合一组样本数据的曲线。

做好准备

对于这个配方，我们将需要导入 NumPy 包，通常作为np，并导入 Matplotlib pyplot模块作为plt：

import numpy as np
import matplotlib.pyplot as plt

我们还需要从 NumPy random模块中导入默认随机数生成器的实例，如下所示：

from numy.random import default_rng
rng = default_rng(12345)

最后，我们需要从 SciPy optimize模块中的curve_fit例程：

from scipy.optimize import curve_fit

如何做...

以下步骤向您展示如何使用curve_fit例程来拟合一组数据的曲线：

1. 第一步是创建样本数据：

SIZE = 100
x_data = rng.uniform(-3.0, 3.0, size=SIZE)
noise = rng.normal(0.0, 0.8, size=SIZE)
y_data = 2.0*x_data**2 - 4*x_data + noise

1. 接下来，我们生成数据的散点图，以查看是否可以识别数据中的潜在趋势：

fig, ax = plt.subplots()
ax.scatter(x_data, y_data)
ax.set(xlabel="x", ylabel="y", title="Scatter plot of sample data")

我们生成的散点图可以在下面的图像中看到。在这里，我们可以看到数据显然不遇线性趋势（直线）。由于我们知道趋势是多项式，我们的下一个猜测将是二次趋势。这就是我们在这里使用的：

图 9.7：样本数据的散点图。我们可以看到数据不遵循线性趋势

1. 接下来，我们创建一个代表我们希望拟合的模型的函数：

def func(x, a, b, c):
    return a*x**2 + b*x + c

1. 现在，我们可以使用curve_fit例程将模型函数拟合到样本数据中：

coeffs, _ = curve_fit(func, x_data, y_data)
print(coeffs)
# [ 1.99611157 -3.97522213 0.04546998]

1. 最后，我们在散点图上绘制最佳拟合曲线，以评估拟合曲线描述数据的效果如何：

x = np.linspace(-3.0, 3.0, SIZE)
y = func(x, coeffs[0], coeffs[1], coeffs[2])
ax.plot(x, y, "k--")

更新后的散点图可以在下面的图像中看到：

图 9.8：散点图，使用最小二乘法找到的最佳拟合曲线

在这里，我们可以看到我们找到的曲线相当合理地拟合了数据。

它是如何工作的...

curve_fit例程执行最小二乘拟合，将模型的曲线拟合到样本数据中。在实践中，这相当于最小化以下目标函数：

这里，配对（x[i]，y[i]）是样本数据中的点。在这种情况下，我们正在优化一个三维参数空间，每个参数都有一个维度。该例程返回估计的系数-参数空间中的点，其中目标函数被最小化-和包含拟合的协方差矩阵的估计的第二个变量。在这个配方中，我们忽略了这一点。

从curve_fit例程返回的估计协方差矩阵可以用于为估计的参数提供置信区间。这是通过取对角线元素的平方根除以样本大小（在这个配方中为 100）来完成的。这给出了估计的标准误差，当乘以与置信度对应的适当值时，给出了置信区间的大小。（我们在第六章中讨论了置信区间，使用数据和统计。）

您可能已经注意到，curve_fit例程估计的参数接近，但并非完全等于我们用于定义步骤 1中的样本数据的参数。这些参数并非完全相等的原因是我们向数据中添加了正态分布的噪声。在这个示例中，我们知道数据的基本结构是二次的 - 也就是二次多项式 - 而不是其他更神秘的函数。实际上，我们不太可能知道数据的基本结构，这就是我们向样本添加噪声的原因。

还有更多...

SciPy optimize模块中还有另一个用于执行最小二乘拟合的例程，名为least_squares。这个例程的签名略微不太直观，但确实返回了一个包含有关优化过程的更多信息的结果对象。然而，这个例程的设置方式可能更类似于我们在*它是如何工作的...*部分中构造基础数学问题的方式。要使用这个例程，我们定义目标函数如下：

def func(params, x, y):
    return y - (params[0]*x**2 + params[1]*x + params[2])

我们将这个函数与参数空间中的起始估计x0一起传递，例如(1, 0, 0)。目标函数的额外参数func可以使用args关键字参数传递；例如，我们可以使用args=(x_data, y_data)。这些参数被传递到目标函数的x和y参数中。总之，我们可以使用以下调用least_squares来估计参数：

results = least_squares(func, [1, 0, 0], args=(x_data, y_data))

从least_squares例程返回的results对象实际上与本章中描述的其他优化例程返回的对象相同。它包含诸如使用的迭代次数、过程是否成功、详细的错误消息、参数值以及最小值处的目标函数值等细节。

分析简单的两人游戏

博弈论是数学的一个分支，涉及决策和策略分析。它在经济学、生物学和行为科学中有应用。许多看似复杂的情况可以简化为一个相对简单的数学游戏，可以以系统的方式进行分析，找到“最优”解决方案。

博弈论中的一个经典问题是囚徒困境，其原始形式如下：两个同谋被抓住，必须决定是保持沉默还是对对方作证。如果两人都保持沉默，他们都要服刑 1 年；如果一个作证而另一个不作证，作证者将获释，而另一个将服刑 3 年；如果两人都相互作证，他们都将服刑 2 年。每个同谋应该怎么做？事实证明，在对对方有任何合理的不信任的情况下，每个同谋可以做出的最佳选择是作证。采用这种策略，他们将不会服刑或最多服刑 2 年。

由于这本书是关于 Python 的，我们将使用这个经典问题的变体来说明这个问题的普遍性。考虑以下问题：你和你的同事必须为客户编写一些代码。你认为你可以用 Python 更快地编写代码，但你的同事认为他们可以用 C 更快地编写代码。问题是，你应该选择哪种语言来进行项目？

你认为你可以用 Python 代码比 C 代码快 4 倍，所以你用速度 1 写 C，用速度 4 写 Python。你的同事说他们可以比 Python 稍微更快地写 C，所以他们用速度 3 写 C，用速度 2 写 Python。如果你们两个都同意一种语言，那么你们按照你预测的速度编写代码，但如果你们意见不一致，那么更快的程序员的生产力会降低 1。我们可以总结如下：

同事/你	C	Python
C	3 / 1	3 / 2
Python	2 / 1	2 / 4

在这个配方中，我们将学习如何在 Python 中构建一个对象来表示这个简单的双人游戏，然后对这个游戏的结果进行一些基本分析。

准备工作

对于这个配方，我们需要导入 NumPy 包为np，导入 Nashpy 包为nash：

import numpy as np
import nashpy as nash

如何做...

以下步骤向您展示了如何使用 Nashpy 创建和执行一些简单的双人游戏分析：

1. 首先，我们需要创建矩阵，用于保存每个玩家（在这个例子中是您和您的同事）的支付信息。

you = np.array([[1, 3], [1, 4]])
colleague = np.array([[3, 2], [2, 2]])

1. 接下来，我们创建一个Game对象，它保存了由这些支付矩阵表示的游戏：

dilemma = nash.Game(you, colleague)

1. 我们使用索引表示法计算给定选择的效用：

print(dilemma[[1, 0], [1, 0]])  # [1 3]
print(dilemma[[1, 0], [0, 1]])  # [3 2]
print(dilemma[[0, 1], [1, 0]])  # [1 2]
print(dilemma[[0, 1], [0, 1]])  # [4 2]

1. 我们还可以根据做出特定选择的概率计算预期效用：

print(dilemma[[0.1, 0.9], [0.5, 0.5]]) # [2.45 2.05]

它是如何工作的...

在这个配方中，我们构建了一个代表非常简单的双人战略游戏的 Python 对象。 这里的想法是有两个“玩家”需要做决定，每个玩家的选择组合都会给出一个特定的支付值。 我们的目标是找到每个玩家可以做出的最佳选择。 假设玩家同时进行一次移动，即没有人知道对方的选择。 每个玩家都有一个确定他们所做选择的策略。

在步骤 1中，我们创建两个矩阵 - 每个玩家一个 - 分配给每个选择组合的支付值。 这两个矩阵由 Nashpy 的Game类包装，它提供了一个方便且直观（从博弈论的角度来看）的接口来处理游戏。 通过使用索引表示法传递选择，我们可以快速计算给定选择组合的效用。

我们还可以根据一种策略提供预期效用的计算，其中选择是根据某种概率分布随机选择的。 语法与先前描述的确定性情况相同，只是我们为每个选择提供了一个概率向量。 我们根据您选择 Python 的概率计算预期效用为 90％，而您的同事选择 Python 的概率为 50％。 预期速度分别为 2.45 和 2.05。

还有更多...

在 Python 中，还有一种计算博弈论的替代方法。 Gambit 项目是一组用于博弈论计算的工具，具有 Python 接口（www.gambit-project.org/）。 这是一个成熟的项目，建立在 C 库周围，并提供比 Nashpy 更高的性能。

计算纳什均衡

纳什均衡是一个双人战略游戏 - 类似于我们在分析简单的双人游戏配方中看到的游戏 - 代表每个玩家都看到“最佳可能”结果的“稳态”。 但是，这并不意味着与纳什均衡相关联的结果是最好的。 纳什均衡比这更微妙。 纳什均衡的非正式定义如下：在其中没有个别玩家可以改善他们的结果的行动配置，假设所有其他玩家都遵守该配置。

我们将探讨纳什均衡的概念，使用经典的猜拳游戏。 规则如下。 每个玩家可以选择以下选项之一：石头，纸或剪刀。 石头打败剪刀，但输给纸； 纸打败石头，但输给剪刀； 剪刀打败纸，但输给石头。 如果两名玩家做出相同选择的游戏是平局。 在数值上，我们用+1 表示赢，-1 表示输，0 表示平局。 从中，我们可以构建一个双人游戏，并计算该游戏的纳什均衡。

在这个配方中，我们将计算猜拳游戏的纳什均衡。

准备工作

对于这个配方，我们需要导入 NumPy 包为np，导入 Nashpy 包为nash：

import numpy as np
import nashpy as nash

如何做...

以下步骤向您展示了如何计算简单的两人游戏的纳什均衡：

1. 首先，我们需要为每个玩家创建一个收益矩阵。我们将从第一个玩家开始：

rps_p1 = np.array([
    [ 0, -1,  1],  # rock payoff
    [ 1,  0, -1],   # paper payoff
    [-1,  1,  0]   # scissors payoff
])

1. 第二个玩家的收益矩阵是rps_p1的转置：

rps_p2 = rps_p1.transpose()

1. 接下来，我们创建Game对象来表示游戏：

rock_paper_scissors = nash.Game(rps_p1, rps_p2)

1. 我们使用支持枚举算法计算游戏的纳什均衡：

equilibria = rock_paper_scissors.support_enumeration()

1. 我们遍历均衡，并打印每个玩家的策略：

for p1, p2 in equilibria:
    print("Player 1", p1)
    print("Player 2", p2)

这些打印语句的输出如下：

Player 1 [0.33333333 0.33333333 0.33333333]
Player 2 [0.33333333 0.33333333 0.33333333]

它是如何工作的...

纳什均衡在博弈论中非常重要，因为它们允许我们分析战略游戏的结果，并确定有利的位置。它们最早由约翰·F·纳什在 1950 年描述，并在现代博弈论中发挥了重要作用。一个两人游戏可能有许多纳什均衡，但任何有限的两人游戏必须至少有一个。问题在于找到给定游戏的所有可能的纳什均衡。

在这个示例中，我们使用了支持枚举，它有效地枚举了所有可能的策略，并筛选出了那些是纳什均衡的策略。在这个示例中，支持枚举算法只找到了一个纳什均衡，这是一个混合策略。这意味着唯一没有改进的策略是随机选择其中一个选项，每个选项的概率为 1/3。这对于任何玩过石头剪刀布的人来说都不足为奇，因为无论我们做出什么选择，我们的对手都有 1/3 的机会（随机）选择打败我们选择的动作。同样，我们有 1/3 的机会平局或赢得比赛，因此我们在所有这些可能性中的期望值如下：

在不知道我们的对手会选择哪一个选项的情况下，没有办法改进这个预期结果。

还有更多...

Nashpy 软件包还提供了其他计算纳什均衡的算法。具体来说，vertex_enumeration方法在Game对象上使用顶点枚举算法，而lemke_howson_enumeration方法使用Lemke Howson算法。这些替代算法具有不同的特性，对于某些问题可能更有效。

另请参阅

Nashpy 软件包的文档包含了有关所涉及的算法和博弈论的更详细信息。其中包括了一些关于博弈论的文本参考。这些文档可以在nashpy.readthedocs.io/en/latest/找到。

进一步阅读

通常情况下，Numerical Recipes书是数值算法的良好来源。第十章，其他主题，涉及函数的最大化和最小化：

• Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., 2017. Numerical recipes: the art of scientific computing. 3rd ed. Cambridge: Cambridge University Press。

有关优化的更具体信息可以在以下书籍中找到：

• Boyd, S.P. and Vandenberghe, L., 2018. Convex optimization*. Cambridge: Cambridge University Press*。

• *Griva, I., Nash, S., and Sofer, A., 2009. Linear and nonlinear optimization.*2nd ed. Philadelphia: Society for Industrial and Applied Mathematics。

最后，以下书籍是博弈论的良好入门：

• Osborne, M.J., 2017. An introduction to game theory*. Oxford: Oxford University Press*。

第十一章：其他主题

在本章中，我们将讨论一些在本书前几章中没有涉及的主题。这些主题大多涉及不同的计算方式以及优化代码执行的其他方式。其他主题涉及处理特定类型的数据或文件格式。

在前两个内容中，我们将介绍帮助跟踪计算中的单位和不确定性的软件包。这些对于涉及具有直接物理应用的数据的计算非常重要。在下一个内容中，我们将讨论如何从 NetCDF 文件加载和存储数据。NetCDF 通常用于存储天气和气候数据的文件格式。（NetCDF 代表网络通用数据格式。）在第四个内容中，我们将讨论处理地理数据，例如可能与天气或气候数据相关的数据。之后，我们将讨论如何可以在不必启动交互式会话的情况下从终端运行 Jupyter 笔记本。接下来的两个内容涉及验证数据和处理从 Kafka 服务器流式传输的数据。我们最后两个内容涉及两种不同的方式，即使用诸如 Cython 和 Dask 等工具来加速我们的代码。

在本章中，我们将涵盖以下内容：

• 使用 Pint 跟踪单位

• 在计算中考虑不确定性

• 从 NetCDF 文件加载和存储数据

• 处理地理数据

• 将 Jupyter 笔记本作为脚本执行

• 验证数据

• 处理数据流

• 使用 Cython 加速代码

• 使用 Dask 进行分布式计算

让我们开始吧！

技术要求

由于本章包含的内容的性质，需要许多不同的软件包。我们需要的软件包列表如下：

• Pint

• 不确定性

• NetCDF4

• xarray

• GeoPandas

• Geoplot

• Papermill

• Cerberus

• Faust

• Cython

• Dask

所有这些软件包都可以使用您喜欢的软件包管理器（如pip）进行安装：

          python3.8 -m pip install pint uncertainties netCDF4 xarray geopandas
   geoplot papermill cerberus faust cython

安装 Dask 软件包，我们需要安装与软件包相关的各种额外功能。我们可以在终端中使用以下pip命令来执行此操作：

          python3.8 -m pip install dask[complete]

除了这些 Python 软件包，我们还需要安装一些支持软件。对于处理地理数据的内容，GeoPandas 和 Geoplot 库可能需要单独安装许多低级依赖项。详细说明在 GeoPandas 软件包文档中给出，网址为geopandas.org/install.html。

对于处理数据流的内容，我们需要安装 Kafka 服务器。如何安装和运行 Kafka 服务器的详细说明可以在 Apache Kafka 文档页面上找到，网址为kafka.apache.org/quickstart。

对于Cython 加速代码的内容，我们需要安装 C 编译器。如何获取GNU C 编译器（GCC）的说明在 Cython 文档中给出，网址为cython.readthedocs.io/en/latest/src/quickstart/install.html。

本章的代码可以在 GitHub 存储库的Chapter 10文件夹中找到，网址为github.com/PacktPublishing/Applying-Math-with-Python/tree/master/Chapter%2010。

查看以下视频以查看代码的实际操作：bit.ly/2ZMjQVw。

使用 Pint 跟踪单位

在计算中正确跟踪单位可能非常困难，特别是在可以使用不同单位的地方。例如，很容易忘记在不同单位之间进行转换 – 英尺/英寸转换成米 – 或者公制前缀 – 比如将 1 千米转换成 1,000 米。

在这个内容中，我们将学习如何使用 Pint 软件包来跟踪计算中的测量单位。

准备工作

对于这个示例，我们需要 Pint 包，可以按如下方式导入：

import pint

如何做...

以下步骤向您展示了如何使用 Pint 包在计算中跟踪单位：

1. 首先，我们需要创建一个UnitRegistry对象：

ureg = pint.UnitRegistry(system="mks")

1. 要创建带有单位的数量，我们将数字乘以注册对象的适当属性：

distance = 5280 * ureg.feet

1. 我们可以使用其中一种可用的转换方法更改数量的单位：

print(distance.to("miles"))
print(distance.to_base_units())
print(distance.to_base_units().to_compact())

这些print语句的输出如下：

0.9999999999999999 mile
1609.3439999999998 meter
1.6093439999999999 kilometer

1. 我们包装一个例程，使其期望以秒为参数并输出以米为结果：

@ureg.wraps(ureg.meter, ureg.second)
def calc_depth(dropping_time):
    # s = u*t + 0.5*a*t*t
    # u = 0, a = 9.81
    return 0.5*9.81*dropping_time*dropping_time

1. 现在，当我们使用分钟单位调用calc_depth例程时，它会自动转换为秒进行计算：

depth = calc_depth(0.05 * ureg.minute)
print("Depth", depth)
# Depth 44.144999999999996 meter

它是如何工作的...

Pint 包为数字类型提供了一个包装类，为类型添加了单位元数据。这个包装类型实现了所有标准的算术运算，并在这些计算过程中跟踪单位。例如，当我们将长度单位除以时间单位时，我们将得到速度单位。这意味着您可以使用 Pint 来确保在复杂计算后单位是正确的。

UnitRegistry对象跟踪会话中存在的所有单位，并处理不同单位类型之间的转换等问题。它还维护一个度量参考系统，在这个示例中是标准国际系统，以米、千克和秒作为基本单位，表示为mks。

wrap功能允许我们声明例程的输入和输出单位，这允许 Pint 对输入函数进行自动单位转换-在这个示例中，我们将分钟转换为秒。尝试使用没有关联单位或不兼容单位的数量调用包装函数将引发异常。这允许对参数进行运行时验证，并自动转换为例程的正确单位。

还有更多...

Pint 包带有一个大型的预设测量单位列表，涵盖了大多数全球使用的系统。单位可以在运行时定义或从文件加载。这意味着您可以定义特定于您正在使用的应用程序的自定义单位或单位系统。

单位也可以在不同的上下文中使用，这允许在不同单位类型之间轻松转换，这些单位类型通常是不相关的。这可以在需要在计算的多个点之间流畅地移动单位的情况下节省大量时间。

在计算中考虑不确定性

大多数测量设备并不是 100%准确的，通常只能准确到一定程度，通常在 0 到 10%之间。例如，温度计可能准确到 1%，而一对数字卡尺可能准确到 0.1%。在这两种情况下，报告的数值不太可能是真实值，尽管它会非常接近。跟踪数值的不确定性是困难的，特别是当您有多种不同的不确定性以不同的方式组合在一起时。与其手动跟踪这些，最好使用一个一致的库来为您完成。这就是uncertainties包的作用。

在这个示例中，我们将学习如何量化变量的不确定性，并看到这些不确定性如何通过计算传播。

准备工作

对于这个示例，我们将需要uncertainties包，我们将从中导入ufloat类和umath模块：

from uncertainties import ufloat, umath

如何做...

以下步骤向您展示了如何在计算中对数值的不确定性进行量化：

1. 首先，我们创建一个不确定的浮点值为3.0加减0.4：

seconds = ufloat(3.0, 0.4)
print(seconds)  # 3.0+/-0.4

1. 接下来，我们进行涉及这个不确定值的计算，以获得一个新的不确定值：

depth = 0.5*9.81*seconds*seconds
print(depth)  # 44+/-12

1. 接下来，我们创建一个新的不确定浮点值，并在与之前计算相反的方向上应用umath模块的sqrt例程：

other_depth = ufloat(44, 12)
time = umath.sqrt(2.0*other_depth/9.81)
print("Estimated time", time)
# Estimated time 3.0+/-0.4

它是如何工作的...

ufloat类包装了float对象，并在整个计算过程中跟踪不确定性。该库利用线性误差传播理论，使用非线性函数的导数来估计计算过程中传播的误差。该库还正确处理相关性，因此从自身减去一个值会得到 0，没有误差。

要跟踪标准数学函数中的不确定性，您需要使用umath模块中提供的版本，而不是 Python 标准库或第三方包（如 NumPy）中定义的版本。

还有更多...

uncertainties包支持 NumPy，并且前面示例中提到的 Pint 包可以与不确定性结合使用，以确保正确地将单位和误差边界归因于计算的最终值。例如，我们可以从本示例的步骤 2中计算出计算的单位，如下所示：

import pint
from uncertainties import ufloat
g = 9.81*ureg.meters / ureg.seconds ** 2
seconds = ufloat(3.0, 0.4) * ureg.seconds

depth = 0.5*g*seconds**2
print(depth)

如预期的那样，最后一行的print语句给出了我们预期的44+/-12 米。

从 NetCDF 文件加载和存储数据

许多科学应用程序要求我们以稳健的格式开始大量的多维数据。NetCDF 是天气和气候行业用于开发数据的格式的一个例子。不幸的是，数据的复杂性意味着我们不能简单地使用 Pandas 包的实用程序来加载这些数据进行分析。我们需要netcdf4包来能够读取和导入数据到 Python 中，但我们还需要使用xarray。与 Pandas 库不同，xarray可以处理更高维度的数据，同时仍提供类似于 Pandas 的接口。

在这个示例中，我们将学习如何从 NetCDF 文件中加载数据并存储数据。

准备就绪

对于这个示例，我们需要导入 NumPy 包作为np，Pandas 包作为pd，Matplotlib pyplot模块作为plt，以及从 NumPy 导入默认随机数生成器的实例：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from numpy.random import default_rng
rng = default_rng(12345)

我们还需要导入xarray包并使用别名xr。您还需要安装 Dask 包，如“技术要求”部分所述，以及 NetCDF4 包：

import xarray as xr

我们不需要直接导入这两个包。

操作方法...

按照以下步骤加载和存储样本数据到 NetCDF 文件中：

1. 首先，我们需要创建一些随机数据。这些数据包括一系列日期、位置代码列表和随机生成的数字：

dates = pd.date_range("2020-01-01", periods=365, name="date")
locations = list(range(25))
steps = rng.normal(0, 1, size=(365,25))
accumulated = np.add.accumulate(steps)

1. 接下来，我们创建一个包含数据的 xarray Dataset对象。日期和位置是索引，而steps和accumulated变量是数据：

data_array = xr.Dataset({
    "steps": (("date", "location"), steps),
    "accumulated": (("date", "location"), accumulated)
    },
    {"location": locations, "date": dates}
)
print(data_array)

print语句的输出如下所示：

<xarray.Dataset>
Dimensions: (date: 365, location: 25)
Coordinates:
* location (location) int64 0 1 2 3 4 5 6 7 8 ... 17 18 19 20 21 22 23 24
* date (date) datetime64[ns] 2020-01-01 2020-01-02 ... 2020-12-30
Data variables:
steps (date, location) float64 geoplot.pointplot(cities, ax=ax, fc="r", marker="2")
ax.axis((-180, 180, -90, 90))-1.424 1.264 ... -0.4547 -0.4873
accumulated (date, location) float64 -1.424 1.264 -0.8707 ... 8.935 -3.525

1. 接下来，我们计算每个时间索引处所有位置的平均值：

means = data_array.mean(dim="location")

1. 现在，我们在新的坐标轴上绘制平均累积值：

fig, ax = plt.subplots()
means["accumulated"].to_dataframe().plot(ax=ax)
ax.set(title="Mean accumulated values", xlabel="date", ylabel="value")

生成的绘图如下所示：

图 10.1：随时间累积平均值的绘图

1. 使用to_netcdf方法将此数据集保存到新的 NetCDF 文件中：

data_array.to_netcdf("data.nc")

1. 现在，我们可以使用xarray的load_dataset例程加载新创建的 NetCDF 文件：

new_data = xr.load_dataset("data.nc")
print(new_data)

前面代码的输出如下所示：

<xarray.Dataset>
Dimensions: (date: 365, location: 25)
Coordinates:
  * location (location) int64 0 1 2 3 4 5 6 7 8 ... 17 18 19 20 21 22 23 24
  * date (date) datetime64[ns] 2020-01-01 2020-01-02 ... 2020-12-30
Data variables:
    steps (date, location) float64 -1.424 1.264 ... -0.4547 -0.4873
    accumulated (date, location) float64 -1.424 1.264 -0.8707 ... 8.935 -3.525

工作原理...

xarray包提供了DataArray和DataSet类，它们（粗略地说）是 PandasSeries和DataFrame对象的多维等价物。在本例中，我们使用数据集，因为每个索引（日期和位置的元组）都与两个数据相关联。这两个对象都暴露了与它们的 Pandas 等价物类似的接口。例如，我们可以使用mean方法沿着其中一个轴计算平均值。DataArray和DataSet对象还有一个方便的方法，可以将其转换为 PandasDataFrame，称为to_dataframe。我们在这个示例中使用它将其转换为DataFrame进行绘图，这并不是真正必要的，因为xarray内置了绘图功能。

这个配方的真正重点是to_netcdf方法和load_dataset例程。前者将DataSet存储在 NetCDF 格式文件中。这需要安装 NetCDF4 包，因为它允许我们访问相关的 C 库来解码 NetCDF 格式的文件。load_dataset例程是一个通用的例程，用于从各种文件格式（包括 NetCDF，这同样需要安装 NetCDF4 包）将数据加载到DataSet对象中。

还有更多...

xarray包支持除 NetCDF 之外的许多数据格式，如 OPeNDAP、Pickle、GRIB 和 Pandas 支持的其他格式。

处理地理数据

许多应用涉及处理地理数据。例如，当跟踪全球天气时，我们可能希望在地图上以各种传感器在世界各地的位置测量的温度为例进行绘图。为此，我们可以使用 GeoPandas 包和 Geoplot 包，这两个包都允许我们操纵、分析和可视化地理数据。

在这个配方中，我们将使用 GeoPandas 和 Geoplot 包来加载和可视化一些样本地理数据。

准备工作

对于这个配方，我们需要 GeoPandas 包，Geoplot 包和 Matplotlib 的pyplot包作为plt导入：

import geopandas
import geoplot
import matplotlib.pyplot as plt

如何做...

按照以下步骤，使用样本数据在世界地图上创建首都城市的简单绘图：

1. 首先，我们需要从 GeoPandas 包中加载样本数据，其中包含世界地理信息：

world = geopandas.read_file(
        geopandas.datasets.get_path("naturalearth_lowres")
)

1. 接下来，我们需要加载包含世界各个首都城市名称和位置的数据：

cities = geopandas.read_file(
        geopandas.datasets.get_path("naturalearth_cities")
)

1. 现在，我们可以创建一个新的图形，并使用polyplot例程绘制世界地理的轮廓：

fig, ax = plt.subplots()
geoplot.polyplot(world, ax=ax)

1. 最后，我们使用pointplot例程在世界地图上添加首都城市的位置。我们还设置轴限制，以使整个世界可见：

geoplot.pointplot(cities, ax=ax, fc="r", marker="2")
ax.axis((-180, 180, -90, 90))

结果绘制的世界各国首都城市的位置如下：

图 10.2：世界首都城市在地图上的绘图

工作原理...

GeoPandas 包是 Pandas 的扩展，用于处理地理数据，而 Geoplot 包是 Matplotlib 的扩展，用于绘制地理数据。GeoPandas 包带有一些我们在这个配方中使用的样本数据集。naturalearth_lowres包含描述世界各国边界的几何图形。这些数据不是非常高分辨率，正如其名称所示，这意味着地理特征的一些细节可能在地图上不存在（一些小岛根本没有显示）。naturalearth_cities包含世界各国首都城市的名称和位置。我们使用datasets.get_path例程来检索包数据目录中这些数据集的路径。read_file例程将数据导入 Python 会话。

Geoplot 包提供了一些专门用于绘制地理数据的附加绘图例程。polyplot例程从 GeoPandas DataFrame 绘制多边形数据，该数据可能描述一个国家的地理边界。pointplot例程从 GeoPandas DataFrame 在一组轴上绘制离散点，这种情况下描述了首都城市的位置。

将 Jupyter 笔记本作为脚本执行

Jupyter 笔记本是用于编写科学和数据应用的 Python 代码的流行媒介。 Jupyter 笔记本实际上是一个以JavaScript 对象表示（JSON）格式存储在带有ipynb扩展名的文件中的块序列。每个块可以是多种不同类型之一，例如代码或标记。这些笔记本通常通过解释块并在后台内核中执行代码然后将结果返回给 Web 应用程序的 Web 应用程序访问。如果您在个人 PC 上工作，这很棒，但是如果您想在服务器上远程运行笔记本中包含的代码怎么办？在这种情况下，甚至可能无法访问 Jupyter 笔记本软件提供的 Web 界面。papermill 软件包允许我们从命令行参数化和执行笔记本。

在本教程中，我们将学习如何使用 papermill 从命令行执行 Jupyter 笔记本。

准备工作

对于本教程，我们需要安装 papermill 软件包，并且当前目录中需要有一个示例 Jupyter 笔记本。我们将使用本章的代码存储库中存储的sample.ipynb笔记本文件。

如何做...

按照以下步骤使用 papermill 命令行界面远程执行 Jupyter 笔记本：

1. 首先，我们从本章的代码存储库中打开样本笔记本sample.ipynb。笔记本包含三个代码单元格，其中包含以下代码：

import matplotlib.pyplot as plt
from numpy.random import default_rng
rng = default_rng(12345)

uniform_data = rng.uniform(-5, 5, size=(2, 100))

fig, ax = plt.subplots(tight_layout=True)
ax.scatter(uniform_data[0, :], uniform_data[1, :])
ax.set(title="Scatter plot", xlabel="x", ylabel="y")

1. 接下来，我们在终端中打开包含 Jupyter 笔记本的文件夹并使用以下命令：

          papermill --kernel python3 sample.ipynb output.ipynb

1. 现在，我们打开输出文件output.ipynb，该文件现在应该包含已更新为执行代码结果的笔记本。在最终块中生成的散点图如下所示：

图 10.3：在远程使用 papermill 执行的 Jupyter 笔记本中生成的随机数据的散点图

它是如何工作的...

papermill 软件包提供了一个简单的命令行界面，用于解释和执行 Jupyter 笔记本，然后将结果存储在新的笔记本文件中。在本教程中，我们提供了第一个参数 - 输入笔记本文件 - sample.ipynb和第二个参数 - 输出笔记本文件 - output.ipynb。然后工具执行笔记本中包含的代码并生成输出。笔记本文件格式跟踪上次运行的结果，因此这些结果将添加到输出笔记本并存储在所需的位置。在本教程中，这是一个简单的本地文件，但是 papermill 也可以存储到云位置，例如Amazon Web Services（AWS）S3 存储或 Azure 数据存储。

在步骤 2中，我们在使用 papermill 命令行界面时添加了--kernel python3选项。此选项允许我们指定用于执行 Jupyter 笔记本的内核。如果 papermill 尝试使用与用于编写笔记本的内核不同的内核执行笔记本，则可能需要这样做以防止错误。可以使用以下命令在终端中找到可用内核的列表：

          jupyter kernelspec list

如果在执行笔记本时出现错误，您可以尝试切换到不同的内核。

还有更多...

Papermill 还具有 Python 接口，因此您可以从 Python 应用程序内执行笔记本。这对于构建需要能够在外部硬件上执行长时间计算并且结果需要存储在云中的 Web 应用程序可能很有用。它还具有向笔记本提供参数的能力。为此，我们需要在笔记本中创建一个标有默认值的参数标记的块。然后可以通过命令行界面使用-p标志提供更新的参数，后跟参数的名称和值。

验证数据

数据通常以原始形式呈现，可能包含异常或不正确或格式不正确的数据，这显然会给后续处理和分析带来问题。通常最好在处理管道中构建验证步骤。幸运的是，Cerberus 包为 Python 提供了一个轻量级且易于使用的验证工具。

对于验证，我们必须定义一个模式，这是关于数据应该如何以及应该对数据执行哪些检查的技术描述。例如，我们可以检查类型并设置最大和最小值的边界。Cerberus 验证器还可以在验证步骤中执行类型转换，这使我们可以将直接从 CSV 文件加载的数据插入验证器中。

在这个示例中，我们将学习如何使用 Cerberus 验证从 CSV 文件加载的数据。

准备工作

对于这个示例，我们需要从 Python 标准库中导入csv模块，以及 Cerberus 包：

import csv
import cerberus

我们还需要这一章的代码库中的sample.csv文件。

如何做...

在接下来的步骤中，我们将使用 Cerberus 包从 CSV 中加载的一组数据进行验证：

1. 首先，我们需要构建描述我们期望的数据的模式。为此，我们必须为浮点数定义一个简单的模式：

float_schema = {"type": "float", "coerce": float, "min": -1.0,
   "max": 1.0}

1. 接下来，我们为单个项目构建模式。这些将是我们数据的行：

item_schema = {
    "type": "dict",
    "schema": {
        "id": {"type": "string"},
        "number": {"type": "integer", "coerce": int},
        "lower": float_schema,
        "upper": float_schema,
    }
}

1. 现在，我们可以定义整个文档的模式，其中将包含一系列项目：

schema = {
    "rows": {
        "type": "list",
        "schema": item_schema
    }
}

1. 接下来，我们使用刚刚定义的模式创建一个Validator对象：

validator = cerberus.Validator(schema)

1. 然后，我们使用csv模块中的DictReader加载数据：

with open("sample.csv") as f:
    dr = csv.DictReader(f)
    document = {"rows": list(dr)}

1. 接下来，我们使用Validator上的validate方法来验证文档：

validator.validate(document)

1. 然后，我们从Validator对象中检索验证过程中的错误：

errors = validator.errors["rows"][0]

1. 最后，我们可以打印出任何出现的错误消息：

for row_n, errs in errors.items():
    print(f"row {row_n}: {errs}")

错误消息的输出如下：

row 11: [{'lower': ['min value is -1.0']}]
row 18: [{'number': ['must be of integer type', "field 'number' cannot be coerced: invalid literal for int() with base 10: 'None'"]}]
row 32: [{'upper': ['min value is -1.0']}]
row 63: [{'lower': ['max value is 1.0']}]

它是如何工作的...

我们创建的模式是对我们需要根据数据检查的所有标准的技术描述。这通常被定义为一个字典，其中项目的名称作为键，属性字典作为值，例如字典中的值的类型或值的边界。例如，在步骤 1中，我们为浮点数定义了一个模式，限制了数字的范围，使其在-1 和 1 之间。请注意，我们包括coerce键，该键指定在验证期间应将值转换为的类型。这允许我们传入从 CSV 文档中加载的数据，其中只包含字符串，而不必担心其类型。

Validator对象负责解析文档，以便对其进行验证，并根据模式描述的所有标准检查它们包含的数据。在这个示例中，我们在创建Validator对象时向其提供了模式。但是，我们也可以将模式作为第二个参数传递给validate方法。错误存储在一个嵌套字典中，其结构与文档的结构相似。

处理数据流

一些数据是从各种来源以恒定流的形式接收的。例如，我们可能会遇到多个温度探头通过 Kafka 服务器定期报告值的情况。Kafka 是一个流数据消息代理，根据主题将消息传递给不同的处理代理。

处理流数据是异步 Python 的完美应用。这使我们能够同时处理更大量的数据，这在应用程序中可能非常重要。当然，在异步上下文中我们不能直接对这些数据进行长时间的分析，因为这会干扰事件循环的执行。

使用 Python 的异步编程功能处理 Kafka 流时，我们可以使用 Faust 包。该包允许我们定义异步函数，这些函数将充当处理代理或服务，可以处理或以其他方式与来自 Kafka 服务器的数据流进行交互。

在这个食谱中，我们将学习如何使用 Faust 包来处理来自 Kafka 服务器的数据流。

准备工作

与本书中大多数食谱不同，由于我们将从命令行运行生成的应用程序，因此无法在 Jupyter 笔记本中运行此食谱。

对于这个食谱，我们需要导入 Faust 包：

import faust

我们还需要从 NumPy 包中运行默认随机数生成器的实例：

from numpy.random import default_rng
rng = default_rng(12345)

我们还需要在本地机器上运行 Kafka 服务的实例，以便我们的 Faust 应用程序可以与消息代理进行交互。

一旦您下载了 Kafka 并解压了下载的源代码，就导航到 Kafka 应用程序所在的文件夹。在终端中打开此文件夹。使用以下命令启动 ZooKeeper 服务器（适用于 Linux 或 Mac）：

          bin/zookeeper-server-start.sh config/zookeeper.properties

如果您使用 Windows，改用以下命令：

          bin\windows\zookeeper-server-start.bat config\zookeeper.properties

然后，在一个新的终端中，使用以下命令启动 Kafka 服务器（适用于 Linux 或 Mac）：

          bin/kafka-server-start.sh config/server.properties

如果您使用 Windows，改用以下命令：

          bin\windows\kafka-server-start.bat config\server.properties

在每个终端中，您应该看到一些日志信息，指示服务器正在运行。

操作步骤...

按照以下步骤创建一个 Faust 应用程序，该应用程序将读取（和写入）数据到 Kafka 服务器并进行一些简单的处理：

1. 首先，我们需要创建一个 FaustApp实例，它将充当 Python 和 Kafka 服务器之间的接口：

app = faust.App("sample", broker="kafka://localhost")

1. 接下来，我们将创建一个记录类型，模拟我们从服务器期望的数据：

class Record(faust.Record):
    id_string: str
    value: float

1. 现在，我们将向 FaustApp对象添加一个主题，将值类型设置为我们刚刚定义的Record类：

topic = app.topic("sample-topic", value_type=Record)

1. 现在，我们定义一个代理，这是一个包装在App对象上的agent装饰器的异步函数：

@app.agent(topic)
async def process_record(records):
    async for record in records:
        print(f"Got {record.id_string}: {record.value}")

1. 接下来，我们定义两个源函数，将记录发布到我们设置的样本主题的 Kafka 服务器上。这些是异步函数，包装在timer装饰器中，并设置适当的间隔：

@app.timer(interval=1.0)
async def producer1(app):
    await app.send(
        "sample-topic",
        value=Record(id_string="producer 1", value=
            rng.uniform(0, 2))
    )

@app.timer(interval=2.0)
async def producer2(app):
    await app.send(
        "sample-topic",
        value=Record(id_string="producer 2", value=
            rng.uniform(0, 5))
    )

1. 在文件底部，我们启动应用程序的main函数：

app.main()

1. 现在，在一个新的终端中，我们可以使用以下命令启动应用程序的工作进程（假设我们的应用程序存储在working-with-data-streams.py中）：

          python3.8 working-with-data-streams.py worker

在这个阶段，您应该看到代理生成的一些输出被打印到终端中，如下所示：

[2020-06-21 14:15:27,986] [18762] [WARNING] Got producer 1: 0.4546720449343393 
[2020-06-21 14:15:28,985] [18762] [WARNING] Got producer 2: 1.5837916985487643 
[2020-06-21 14:15:28,989] [18762] [WARNING] Got producer 1: 1.5947309146654682 
[2020-06-21 14:15:29,988] [18762] [WARNING] Got producer 1: 1.3525093415019491

这将是由 Faust 生成的一些应用程序信息的下方。

1. 按下Ctrl + C关闭工作进程，并确保以相同的方式关闭 Kafka 服务器和 Zookeeper 服务器。

工作原理...

这是 Faust 应用程序的一个非常基本的示例。通常，我们不会生成记录并通过 Kafka 服务器发送它们，并在同一个应用程序中处理它们。但是，这对于本演示来说是可以的。在生产环境中，我们可能会连接到远程 Kafka 服务器，该服务器连接到多个来源并同时发布到多个不同的主题。

Faust 应用程序控制 Python 代码与 Kafka 服务器之间的交互。我们使用agent装饰器添加一个函数来处理发布到特定通道的信息。每当新数据被推送到样本主题时，将执行此异步函数。在这个食谱中，我们定义的代理只是将Record对象中包含的信息打印到终端中。

timer装饰器定义了一个服务，定期在指定的间隔执行某些操作。在我们的情况下，计时器通过App对象向 Kafka 服务器发送消息。然后将这些消息推送给代理进行处理。

Faust 命令行界面用于启动运行应用程序的工作进程。这些工作进程实际上是在 Kafka 服务器上或本地进程中对事件做出反应的处理者，例如本示例中定义的定时器服务。较大的应用程序可能会使用多个工作进程来处理大量数据。

此外

Faust 文档提供了有关 Faust 功能的更多详细信息，以及 Faust 的各种替代方案：faust.readthedocs.io/en/latest/。

有关 Kafka 的更多信息可以在 Apache Kafka 网站上找到：kafka.apache.org/。

使用 Cython 加速代码

Python 经常因为速度慢而受到批评——这是一个无休止的争论。使用具有 Python 接口的高性能编译库（例如科学 Python 堆栈）可以解决许多这些批评，从而大大提高性能。然而，在某些情况下，很难避免 Python 不是编译语言的事实。在这些（相当罕见的）情况下，改善性能的一种方法是编写 C 扩展（甚至完全重写代码为 C）以加速关键部分。这肯定会使代码运行更快，但可能会使维护软件包变得更加困难。相反，我们可以使用 Cython，这是 Python 语言的扩展，可以转换为 C 并编译以获得更好的性能改进。

例如，我们可以考虑一些用于生成 Mandelbrot 集图像的代码。为了比较，我们假设纯 Python 代码——我们假设这是我们的起点——如下所示：

# mandelbrot/python_mandel.py

import numpy as np

def in_mandel(cx, cy, max_iter):
    x = cx
    y = cy
    for i in range(max_iter):
        x2 = x**2
        y2 = y**2
        if (x2 + y2) >= 4:
            return i
        y = 2.0*x*y + cy
        x = x2 - y2 + cx
    return max_iter

def compute_mandel(N_x, N_y, N_iter):
    xlim_l = -2.5
    xlim_u = 0.5
    ylim_l = -1.2
    ylim_u = 1.2
    x_vals = np.linspace(xlim_l, xlim_u, N_x, dtype=np.float64)
    y_vals = np.linspace(ylim_l, ylim_u, N_y, dtype=np.float64)

    height = np.empty((N_x, N_y), dtype=np.int64)
    for i in range(N_x):
        for j in range(N_y):
            height[i, j] = in_mandel(x_vals[i], y_vals[j], N_iter)
    return height

纯 Python 中这段代码相对较慢的原因是相当明显的：嵌套循环。为了演示目的，让我们假设我们无法使用 NumPy 对这段代码进行矢量化。一些初步测试显示，使用这些函数生成 Mandelbrot 集的 320×240 点和 255 步大约需要 6.3 秒。您的时间可能会有所不同，这取决于您的系统。

在这个示例中，我们将使用 Cython 大大提高前面代码的性能，以生成 Mandelbrot 集图像。

准备工作

对于这个示例，我们需要安装 NumPy 包和 Cython 包。您还需要在系统上安装 GCC 等 C 编译器。例如，在 Windows 上，您可以通过安装 MinGW 来获取 GCC 的版本。

操作步骤

按照以下步骤使用 Cython 大大提高生成 Mandelbrot 集图像的代码性能：

1. 在mandelbrot文件夹中创建一个名为cython_mandel.pyx的新文件。在这个文件中，我们将添加一些简单的导入和类型定义：

# mandelbrot/cython_mandel.pyx

import numpy as np
cimport numpy as np
cimport cython
ctypedef Py_ssize_t Int
ctypedef np.float64_t Double

1. 接下来，我们使用 Cython 语法定义in_mandel例程的新版本。我们在这个例程的前几行添加了一些声明：

cdef int in_mandel(Double cx, Double cy, int max_iter):
    cdef Double x = cx
    cdef Double y = cy
    cdef Double x2, y2
    cdef Int i

1. 函数的其余部分与 Python 版本的函数相同：

    for i in range(max_iter):
        x2 = x**2
        y2 = y**2
        if (x2 + y2) >= 4:
            return i
        y = 2.0*x*y + cy
        x = x2 - y2 + cx
    return max_iter

1. 接下来，我们定义compute_mandel函数的新版本。我们向这个函数添加了 Cython 包的两个装饰器：

@cython.boundscheck(False)
@cython.wraparound(False)
def compute_mandel(int N_x, int N_y, int N_iter):

1. 然后，我们像在原始例程中一样定义常量：

    cdef double xlim_l = -2.5
    cdef double xlim_u = 0.5
    cdef double ylim_l = -1.2
    cdef double ylim_u = 1.2

1. 我们使用 NumPy 包中的linspace和empty例程的方式与 Python 版本完全相同。这里唯一的添加是我们声明了i和j变量，它们是Int类型的：

    cdef np.ndarray x_vals = np.linspace(xlim_l, xlim_u, 
        N_x, dtype=np.float64)
    cdef np.ndarray y_vals = np.linspace(ylim_l, ylim_u, 
        N_y, dtype=np.float64)
    cdef np.ndarray height = np.empty((N_x, N_y), dtype=np.int64)
    cdef Int i, j

1. 定义的其余部分与 Python 版本完全相同：

    for i in range(N_x):
        for j in range(N_y):
            height[i, j] = in_mandel(x_vals[i], y_vals[j], N_iter)
    return height

1. 接下来，在mandelbrot文件夹中创建一个名为setup.py的新文件，并将以下导入添加到此文件的顶部：

# mandelbrot/setup.py

import numpy as np
from setuptools import setup, Extension
from Cython.Build import cythonize

1. 之后，我们使用指向原始python_mandel.py文件的源定义一个扩展模块。将此模块的名称设置为hybrid_mandel：

hybrid = Extension(
    "hybrid_mandel",
    sources=["python_mandel.py"],
    include_dirs=[np.get_include()],
    define_macros=[("NPY_NO_DEPRECATED_API", 
       "NPY_1_7_API_VERSION")]
)

1. 现在，我们定义第二个扩展模块，将源设置为刚刚创建的cython_mandel.pyx文件：

cython = Extension(
    "cython_mandel",
    sources=["cython_mandel.pyx"],
    include_dirs=[np.get_include()],
    define_macros=[("NPY_NO_DEPRECATED_API", 
       "NPY_1_7_API_VERSION")]
)

1. 接下来，将这两个扩展模块添加到列表中，并调用setup例程来注册这些模块：

extensions = [hybrid, cython]
setup(
    ext_modules = cythonize(extensions, compiler_directives=
       {"language_level": "3"}),
)

1. 在mandelbrot文件夹中创建一个名为__init__.py的新空文件，以便将其转换为可以在 Python 中导入的包。

2. 在mandelbrot文件夹中打开终端，并使用以下命令构建 Cython 扩展模块：

          python3.8 setup.py build_ext --inplace

1. 现在，开始一个名为run.py的新文件，并添加以下导入语句：

# run.py

from time import time
from functools import wraps
import matplotlib.pyplot as plt

1. 从我们定义的每个模块中导入各种compute_mandel例程：原始的python_mandel；Cython 化的 Python 代码hybrid_mandel；以及编译的纯 Cython 代码cython_mandel：

from mandelbrot.python_mandel import compute_mandel 
    as compute_mandel_py
from mandelbrot.hybrid_mandel import compute_mandel 
    as compute_mandel_hy
from mandelbrot.cython_mandel import compute_mandel
    as compute_mandel_cy

1. 定义一个简单的计时器装饰器，我们将用它来测试例程的性能：

def timer(func, name):
    @wraps(func)
    def wrapper(*args, **kwargs):
        t_start = time()
        val = func(*args, **kwargs)
        t_end = time()
        print(f"Time taken for {name}: {t_end - t_start}")
        return val
    return wrapper

1. 将timer装饰器应用于每个导入的例程，并定义一些用于测试的常量：

mandel_py = timer(compute_mandel_py, "Python")
mandel_hy = timer(compute_mandel_hy, "Hybrid")
mandel_cy = timer(compute_mandel_cy, "Cython")

Nx = 320
Ny = 240
steps = 255

1. 用我们之前设置的常量运行每个装饰的例程。将最终调用（Cython 版本）的输出记录在vals变量中：

mandel_py(Nx, Ny, steps)
mandel_hy(Nx, Ny, steps)
vals = mandel_cy(Nx, Ny, steps)

1. 最后，绘制 Cython 版本的输出，以检查例程是否正确计算了 Mandelbrot 集：

fig, ax = plt.subplots()
ax.imshow(vals.T, extent=(-2.5, 0.5, -1.2, 1.2))
plt.show()

运行run.py文件将在终端打印每个例程的执行时间，如下所示：

Time taken for Python: 6.276328802108765
Time taken for Hybrid: 5.816391468048096
Time taken for Cython: 0.03116750717163086

Mandelbrot 集的绘图可以在以下图像中看到：

图 10.4：使用 Cython 代码计算的 Mandelbrot 集的图像

这是我们对 Mandelbrot 集的期望。

它是如何工作的...

在这个示例中发生了很多事情，所以让我们从解释整个过程开始。Cython 接受用 Python 语言的扩展编写的代码，并将其编译成 C 代码，然后用于生成可以导入 Python 会话的 C 扩展库。实际上，您甚至可以使用 Cython 直接将普通 Python 代码编译为扩展，尽管结果不如使用修改后的语言好。在这个示例中的前几个步骤中，我们在修改后的语言中定义了 Python 代码的新版本（保存为.pyx文件），其中包括类型信息以及常规 Python 代码。为了使用 Cython 构建 C 扩展，我们需要定义一个设置文件，然后创建一个文件来生成结果。

Cython 代码的最终编译版本比其 Python 等效代码运行速度快得多。Cython 编译的 Python 代码（在本示例中称为混合代码）的性能略优于纯 Python 代码。这是因为生成的 Cython 代码仍然必须处理带有所有注意事项的 Python 对象。通过在.pyx文件中向 Python 代码添加类型信息，我们开始看到性能的重大改进。这是因为in_mandel函数现在有效地被定义为一个 C 级别函数，它不与 Python 对象交互，而是操作原始数据类型。

Cython 代码和 Python 等效代码之间存在一些小但非常重要的区别。在步骤 1中，您可以看到我们像往常一样导入了 NumPy 包，但我们还使用了cimport关键字将一些 C 级别的定义引入了作用域。在步骤 2中，我们在定义in_mandel例程时使用了cdef关键字而不是def关键字。这意味着in_mandel例程被定义为一个 C 级别函数，不能从 Python 级别使用，这在调用这个函数时（这经常发生）节省了大量开销。

关于这个函数定义的唯一其他真正的区别是在签名和函数的前几行中包含了一些类型声明。我们在这里应用的两个装饰器禁用了访问列表（数组）元素时的边界检查。boundscheck装饰器禁用了检查索引是否有效（在 0 和数组大小之间），而wraparound装饰器禁用了负索引。尽管它们禁用了 Python 内置的一些安全功能，但这两个装饰器在执行过程中都会对速度产生适度的改进。在这个示例中，禁用这些检查是可以的，因为我们正在使用循环遍历数组的有效索引。

设置文件是我们告诉 Python（因此也是 Cython）如何构建 C 扩展的地方。Cython 中的cythonize例程在这里起着关键作用，因为它触发了 Cython 构建过程。在步骤 9和10中，我们使用setuptools中的Extension类定义了扩展模块，以便我们可以为构建定义一些额外的细节；具体来说，我们为 NumPy 编译设置了一个环境变量，并添加了 NumPy C 头文件的include文件。这是通过Extension类的define_macros关键字参数完成的。我们在步骤 13中使用setuptools命令来构建 Cython 扩展，并且添加了--inplace选项，这意味着编译后的库将被添加到当前目录，而不是放在一个集中的位置。这对开发来说是很好的。

运行脚本相当简单：从每个定义的模块中导入例程 - 其中两个实际上是 C 扩展模块 - 并计算它们的执行时间。我们必须在导入别名和例程名称上有一些创造性，以避免冲突。

还有更多...

Cython 是改进代码性能的强大工具。然而，在优化代码时，您必须始终谨慎地花费时间。使用像 Python 标准库中提供的 cProfiler 这样的性能分析工具可以用来找到代码中性能瓶颈出现的地方。在这个示例中，性能瓶颈出现的地方是相当明显的。在这种情况下，Cython 是解决问题的良药，因为它涉及对（双重）for循环内的函数进行重复调用。然而，它并不是解决性能问题的通用方法，往往情况下，通过重构代码以利用高性能库，可以大大提高代码的性能。

Cython 与 Jupyter 笔记本集成良好，并且可以无缝地在笔记本的代码块中使用。Cython 也包含在 Python 的 Anaconda 发行版中，因此在使用 Anaconda 发行版安装了 Cython 后，就无需额外设置即可在 Jupyter 笔记本中使用 Cython。

在从 Python 生成编译代码时，Cython 并不是唯一的选择。例如，NumBa 包（numba.pydata.org/）提供了一个即时（JIT）编译器，通过简单地在特定函数上放置装饰器来优化 Python 代码的运行时。NumBa 旨在与 NumPy 和其他科学 Python 库一起使用，并且还可以用于利用 GPU 加速代码。

使用 Dask 进行分布式计算

Dask 是一个用于在多个线程、进程或甚至计算机之间进行分布式计算的库，以有效地进行大规模计算。即使您只是在一台笔记本电脑上工作，这也可以极大地提高性能和吞吐量。Dask 提供了 Python 科学堆栈中大多数数据结构的替代品，如 NumPy 数组和 Pandas DataFrames。这些替代品具有非常相似的接口，但在内部，它们是为分布式计算而构建的，以便它们可以在多个线程、进程或计算机之间共享。在许多情况下，切换到 Dask 就像改变import语句一样简单，可能还需要添加一些额外的方法调用来启动并发计算。

在这个示例中，我们将学习如何使用 Dask 对 DataFrame 进行一些简单的计算。

准备工作

对于这个示例，我们需要从 Dask 包中导入dataframe模块。按照 Dask 文档中的约定，我们将使用别名dd导入此模块：

import dask.dataframe as dd

我们还需要这一章的代码库中的sample.csv文件。

如何做...

按照以下步骤使用 Dask 对 DataFrame 对象执行一些计算：

1. 首先，我们需要将数据从sample.csv加载到 Dask 的DataFrame中：

data = dd.read_csv("sample.csv")

1. 接下来，我们对 DataFrame 的列执行标准计算：

sum_data = data.lower + data.upper
print(sum_data)

与 Pandas DataFrames 不同，结果不是一个新的 DataFrame。print语句给了我们以下信息：

Dask Series Structure:
npartitions=1
    float64
        ...
dtype: float64
Dask Name: add, 6 tasks

1. 要实际获得结果，我们需要使用compute方法：

result = sum_data.compute()
print(result.head())

结果现在如预期所示：

0 -0.911811
1 0.947240
2 -0.552153
3 -0.429914
4 1.229118
dtype: float64

1. 我们计算最后两列的均值的方式与 Pandas DataFrame 完全相同，但我们需要添加一个调用compute方法来执行计算：

means = data.loc[:, ("lower", "upper")].mean().compute()
print(means)

打印的结果与我们的预期完全一致：

lower -0.060393
upper -0.035192
dtype: float64

它是如何工作的...

Dask 为计算构建了一个任务图，描述了需要对数据集合执行的各种操作和计算之间的关系。这样可以将计算步骤分解，以便可以按正确的顺序在不同的工作器之间进行计算。然后将此任务图传递给调度程序，调度程序将实际任务发送给工作器执行。Dask 配备了几种不同的调度程序：同步、线程、多进程和分布式。可以在compute方法的调用中选择调度程序的类型，或者全局设置。如果没有给出一个合理的默认值，Dask 会选择一个合理的默认值。

同步、线程和多进程调度程序在单台机器上工作，而分布式调度程序用于与集群一起工作。Dask 允许您以相对透明的方式在调度程序之间切换，尽管对于小任务，您可能不会因为设置更复杂的调度程序而获得任何性能优势。

compute方法是这个示例的关键。通常会在 Pandas DataFrames 上执行计算的方法现在只是设置了一个通过 Dask 调度程序执行的计算。直到调用compute方法之前，计算才会开始。这类似于Future作为异步函数调用结果的代理返回，直到计算完成才会实现。

还有更多...

Dask 提供了 NumPy 数组的接口，以及本示例中显示的 DataFrames。还有一个名为dask_ml的机器学习接口，它提供了类似于 scikit-learn 包的功能。一些外部包，如xarray，也有 Dask 接口。Dask 还可以与 GPU 一起工作，以进一步加速计算并从远程源加载数据，这在计算分布在集群中时非常有用。