IPython-笔记本精要-二-

IPython 笔记本精要（二）

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

译者：飞龙

协议：CC BY-NC-SA 4.0

第五章。使用 SciPy、Numba 和 NumbaPro 的高级计算

本章中，用户将学习如何使用SciPy执行科学计算。然后，将介绍Numba包作为加速计算的方法。最后，将介绍NumbaPro的 GPU 并行执行能力。

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

• SciPy概述

• 使用SciPy的高级数学算法

• 使用Numba和NumbaPro加速计算

在运行本章中的示例之前，请通过在计算单元中运行以下命令来加载pylab：

%pylab inline

SciPy概述

SciPy是一个广泛应用于数学和科学计算的库。以下是库中所有可用模块的完整列表：

模块	功能
cluster	聚类算法
constants	物理和数学常数
fftpack	快速傅里叶变换
integrate	积分与常微分方程
interpolate	插值与样条
io	输入与输出
linalg	线性代数
ndimage	图像处理
odr	正交距离回归
optimize	优化与根求解
signal	信号处理
sparse	稀疏矩阵
spatial	空间数据结构
special	特殊函数
stats	统计分布
weave	C/C++ 集成

在脚本中导入SciPy模块的标准方式是使用以下命令行：

from scipy import signal

然后，可以使用通常的模块引用语法来调用各个函数，如下所示：

signal.correlate(…)

然而，许多最常用的函数可以在SciPy层次结构的顶层找到。此外，我们使用 IPython 的交互模式，并使用（本书中假设的）魔法命令，如下所示：

%pylab inline

许多函数将无需显式引用模块即可使用。

在下一节中，我们将展示SciPy中可用的一些函数示例。读者不需要深入了解示例中将使用的数学技术和算法。

使用 SciPy 的高级数学算法

在本节中，我们将涵盖SciPy中一些可用的算法。以下每个子节都包含一个来自应用科学重要领域的代表性示例。这些示例的选择不需要广泛的领域知识，但仍然具有现实性。我们呈现的主题和示例如下：

• 解方程和寻找最优值：我们将研究一个市场模型，它需要求解一个非线性系统，以及一个需要非标准优化的设施选址问题。

• 微积分和常微分方程：我们将展示一个利用积分微积分进行的体积计算，并介绍牛顿的炮，这是由艾萨克·牛顿提出的一个思想实验，我们将用常微分方程系统对其建模。最后，我们将介绍一个三维系统——著名的洛伦兹方程，这是早期表现出混沌行为的例子。

求解方程和寻找最优值

为了说明这个主题，我们使用了经济学中一个标准的供需模型。在这个模型中，供需与价格之间有函数关系，市场均衡通过确定供需曲线的交点来找到。我们在示例中使用的数学公式有些是任意的（因此可能不现实），但这些公式超出了教科书中通常假设供需关系为线性关系的范围。

定义供需曲线的公式如下：

我们将使用函数工厂模式。在一个单元格中运行以下代码：

def make_supply(A, B, C):
 def supply_func(q):
 return A * q / (C  - B * q)
 return supply_func
def make_demand(K, L):
 def demand_func(q):
 return K / (1 + L * q)
 return demand_func

上述代码并没有直接定义供需曲线，而是指定了函数工厂。这样做可以更方便地处理参数，这在应用问题中非常常见，因为我们希望相同的模型能适用于多种情况。

接下来，我们设置参数值，并调用函数工厂来定义实际评估供需曲线的函数，如下所示：

A, B, C = 23.3, 9.2, 82.4
K, L = 1.2, 0.54
supply = make_supply(A, B, C)
demand = make_demand(K, L)

以下代码行生成了曲线的图像：

q = linspace(0.01,5,200)
plot(q, supply(q), lw = 2)
plot(q, demand(q), lw = 2)
title('Supply and demand curves')
xlabel('Quantity (thousands of units)')
ylabel('Price ($)')
legend(['Supply', 'Demand'], loc='upper left')

以下是上述代码行输出的图像：

选择的供需曲线反映了合理的假设：随着价格上涨，供应增加而需求减少。即使价格为零，需求也是有限的（反映了市场中对产品感兴趣的人口有限）。另一方面，供给曲线具有垂直渐近线（在图中未显示），表明存在生产限制（即使价格趋近于无限，也只能提供有限的市场供应量）。

市场的均衡点是供需曲线的交点。为了找到均衡点，我们使用optimize模块，该模块除了提供优化函数外，还具有求解数值方程的功能。推荐用来求解一元函数的函数是brentq()，如下代码所示：

from scipy import optimize
def opt_func(q):
 return supply(q) - demand(q)
q_eq = optimize.brentq(opt_func, 1.0, 2.0)
print q_eq, supply(q_eq), demand(q_eq)

brentq()函数假设我们要求解的方程的右边是0。因此，我们首先定义opt_func()函数，该函数计算供给和需求之间的差值。这个函数是brentq()的第一个参数。接下来的两个数值参数给出了包含解的区间。选择一个确切包含方程解的区间非常重要。在我们的例子中，通过查看图形很容易做到这一点，因为可以明显看到曲线在 1 和 2 之间交叉。

运行前面的代码会产生以下输出：

1.75322153719 0.616415252177 0.616415252177

第一个值是均衡点，即在最优价格下能够销售的单位数（以千为单位）。最优价格是通过供给曲线和需求曲线共同计算得出的（以检查这些值是否确实相同）。

为了说明一个有两个变量的优化问题，我们考虑一个最优设施位置问题。假设一座工厂有若干个制造站点，需要从一个单一的供给站点分配材料。工厂的车间是矩形的，分配轨道必须与工厂的墙壁平行。这最后一个要求使问题变得有趣。需要最小化的函数与所谓的出租车****距离相关，具体如以下图所示：

第一步是定义给定制造站点的位置，如下所示：

points = array([[10.3,15.4],[6.5,8.8],[15.6,10.3],[4.7,12.8]])

这些位置存储为一个 4 x 2 的NumPy数组，名为points，每行代表一个点。以下命令会绘制之前命令行中提到的这些点的图：

plot(points[:,0],points[:,1], 'o', ms=12, mfc='LightSkyBlue')

这些点使用圆形标记显示，标记由参数o指定，该参数还关闭了连接这些点的线段。ms和mfc选项分别指定标记的大小（以像素为单位）和颜色。生成的输出为以下图像：

下一步是定义要最小化的函数。我们再次倾向于定义一个函数工厂的方法，如下所示：

def make_txmin(points):
 def txmin_func(p):
 return sum(abs(points - p))
 return txmin_func

这段代码的关键在于计算出租车距离的方式，充分利用了NumPy数组操作的灵活性。此操作在以下代码行中完成：

return sum(abs(points - p))

这段代码首先计算了向量差值points-p。请注意，points是一个 4 x 2 的数组，而p是一个 1 x 2 的数组。NumPy意识到数组的维度不同，并利用其广播规则。结果是，数组p被从points数组的每一行中减去，这正是我们想要的。接着，abs()函数计算结果数组中每个条目的绝对值，最后sum()函数将所有条目相加。这一切都在一行代码中完成，工作量非常大！

然后我们需要使用函数工厂来定义实际计算出租车距离的函数。

txmin = make_txmin(points)

函数工厂只是简单地用包含实际位置的数组作为参数进行调用。此时，问题已经完全设置好，我们准备计算最优解，以下代码可以完成这一过程：

from scipy import optimize
x0 = array([0.,0.])
res = optimize.minimize(
 txmin, x0,
 method='nelder-mead',
 options={'xtol':1e-5, 'disp':True})

最小化通过调用minimize()函数来计算。该函数的前两个参数是之前单元格中定义的目标函数txmin()和初始猜测值x0。我们选择原点作为初始猜测值，但在实际问题中，我们会利用所收集的信息来选择一个接近实际最小值的猜测值。该函数提供了几种优化方法，适用于不同类型的目标函数。我们使用Nelder-Mead 方法，这是一种启发式算法，不要求目标函数具有光滑性，非常适合当前问题。最后，我们为该方法指定了两个选项：所需的容差和用于打印诊断信息的显示选项。这将生成以下输出：

Optimization terminated successfully.
 Current function value: 23.800000
 Iterations: 87
 Function evaluations: 194

上述输出表示已成功找到最小值，并给出了其值。请注意，与任何数值优化方法一样，通常只能保证找到局部最小值。在这种情况下，由于目标函数是凸的，因此最小值必定是全局最小值。该函数的结果存储在SciPy数据结构中，这种数据结构是optimize模块中定义的OptimizeResult类型。要获取设施的最优位置，我们可以使用以下命令：

print res.x

上述命令的输出如下：

[  8.37782286  11.36247412]

为了完成这个示例，我们展示了显示最优解的代码：

plot(points[:,0],points[:,1], 'o', ms=12, mfc='LightSkyBlue')
plot(res.x[0], res.x[1],'o', ms=12, mfc='LightYellow')
locstr = 'x={:5.2f}, y={:5.2f}'.format(res.x[0], res.x[1])
title('Optimal facility location: {}'.format(locstr))

对plot()函数的调用类似于之前示例中的调用。为了给图表添加格式化良好的标题，我们首先定义locstr字符串，用于显示最优位置坐标。这个字符串是一个 Python 格式化字符串，格式规范为{:5.2f}，即宽度为5且精度为2位的小数。结果如下图所示：

微积分和微分方程

作为微积分计算的一个例子，我们将展示如何计算旋转体的体积。该旋转体是通过将下图中显示的曲线绕y轴旋转得到的：

这条曲线是通过以下代码绘制的：

def make_gen(a,b):
 def gen_func(y):
 return a/pi * arccos(2.0 * y / b - 1.0)
 return gen_func
a = 5.
b = 4.
gen = make_gen(a,b)
x = linspace(0,b,200)
y = gen(x)
subplot(111, aspect='equal')
plot(x,y,lw=2)

该曲线本质上是一个拉伸和转置的反余弦函数，定义在make_gen()函数中。它依赖于两个参数，a和b，分别指定其高度和长度。make_gen()函数是一个函数工厂，返回一个实际计算曲线值的函数。实际定义曲线的函数称为gen()（代表生成器），因此这是要绘制的函数。

当这条曲线围绕垂直轴旋转时，我们得到了如下绘制的实体：

当然，前面的图是使用 IPython 生成的，使用了以下代码：

from mpl_toolkits.mplot3d import Axes3D
na = 50
nr = 50
avalues = linspace(0, 2*pi, na, endpoint=False)
rvalues = linspace(b/nr, b, nr)
avalues_r = repeat(avalues[...,newaxis], nr, axis=1)
xvalues = append(0, (rvalues*cos(avalues_r)).flatten())
yvalues = append(0, (rvalues*sin(avalues_r)).flatten())
zvalues = gen(sqrt(xvalues*xvalues+yvalues*yvalues))
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(xvalues, yvalues, zvalues, 
 color='Cyan',alpha=0.65,linewidth=0.)

该代码中的关键函数是最后一行中对plot_trisurf()的调用。该函数接受三个NumPy数组，xvalues、yvalues和zvalues，指定了表面上点的坐标。数组xvalues和yvalues定义了一系列同心圆上的点，如下图所示：

z坐标的值是通过在每个点计算gen(sqrt(x*x+y*y))得到的，这样就将同心圆上的所有点在 3D 图中分配了相同的高度。

要计算实体的体积，我们使用圆柱壳法。该方法的工作原理超出了本书的范围，但归根结底是计算一个积分，如下公式所示：

在这个公式中，*f(x)*函数代表围绕y轴旋转的曲线。为了计算这个积分，我们使用scipy.integrate包。我们使用quad()函数，适用于没有奇点的函数的通用积分。以下是这个公式的代码：

import scipy.integrate as integrate
int_func  = lambda x: 2 * pi * x * gen(x)
integrate.quad(int_func, 0, b)

导入integrate模块后，我们定义要进行积分的函数。请注意，由于这是一行计算，我们使用了lambda语法。最后，我们调用quad()来执行积分。调用的参数是被积函数和积分的边界（在这种情况下是从0到b）。以下是前面几行代码的输出：

(94.24777961000055, 1.440860870616234e-07)

第一个数字是积分值，第二个是误差估计。

在下一个示例中，我们考虑牛顿的大炮，这是现代物理和微积分根源的一个思维实验。以下图像展示了这种情况，这是艾萨克·牛顿爵士的著作《世界体系论》中的一幅雕刻：

牛顿让我们想象有一门大炮坐落在一座非常高的山顶上。如果大炮发射一颗抛射体，它会飞行一段时间，最终落到地面。抛射体的初速度越大，它落地的距离越远。让我们假设我们可以任意快地发射抛射体，并且没有空气阻力。这样，随着初速度的增加，最终抛射体将绕地球飞行，如果大炮足够迅速地移除，抛射体将永远继续绕地球轨道运行。牛顿用这个例子来解释月球是如何绕地球运行的，而不单单依靠重力作用下掉落。

要建模这个情况，我们需要使用牛顿的引力定律作为一个微分方程组：

我们不会尝试解释这些公式是如何得出的，唯一对我们重要的是有四个状态变量，前两个表示抛射体的位置，后两个表示其速度向量。由于运动发生在通过地球中心的平面内，因此只需要两个位置变量。G和M分别是表示牛顿万有引力常数和地球质量的常量。抛射体的质量不出现，因为引力质量与惯性质量恰好抵消。

使用SciPy来解决这个问题的第一步是定义这一组微分方程，具体实现如下代码：

M = 5.9726E24
G = 6.67384E-11
C = G * M
def ode_func(xvec, t):
 x1, x2, v1, v2 = xvec
 d = (x1 * x1 + x2 * x2) ** 1.5
 return array([v1, v2, -C * x1 / d, -C * x2 / d ])

我们只需要做的是定义一个函数，计算微分方程组的右侧。我们首先定义常量M和G（使用国际单位制），以及辅助常量C，因为G和M仅通过它们的乘积出现在方程中。系统由ode_func()函数表示。这个函数必须至少接受两个参数：一个NumPy数组xvec和一个浮动值t。在我们的例子中，xvec是一个四维向量，因为我们的系统有四个状态变量。t变量在系统中没有被使用，因为没有外力（如果我们发射的是火箭而不是抛射体，那么就会有外力）。然而，它仍然必须作为输入参数列出。

在ode_func()内部，我们首先通过赋值提取xvec向量的各个元素，如下所示：

x1, x2, v1, v2 = xvec

这不是严格必要的，但有助于提高可读性。然后我们计算辅助量d（这是最后两个方程的分母）。最后，输出数组根据系统中的公式进行计算。请注意，没有计算任何导数，因为求解器所需的所有信息都包含在方程的右侧。

我们现在准备使用以下代码行来求解微分方程组：

import scipy.integrate as integrate
earth_radius = 6.371E6
v0 = 8500
h0 = 5E5
ic = array([0.0, earth_radius + h0, v0, 0.0])
tmax = 8700.0
dt = 10.0
tvalues = arange(0.0, tmax, dt)
xsol = integrate.odeint(ode_func, ic, tvalues)

上述代码的第一行导入了integrate模块，该模块用于求解微分方程。然后，我们需要指定抛射物的初始位置和速度。我们假设火炮位于北极，位于 50,000 米高的塔楼上（虽然这显然不现实，我们选择如此大的值是为了增强轨道的可见性）。由于地球不是完美的球体，我们使用半径的平均值。初始速度设置为8500米/秒。

初始条件存储在一个NumPy数组中，通过以下赋值完成：

ic = array([0.0, earth_radius + h0, v0, 0.0])

下一步是定义初始时间（在我们的例子中为零）和解所需的时间数组。这通过以下三行代码完成：

tmax = 8650.0
dt = 60.0
tvalues = arange(0.0, tmax, dt)

我们首先定义tmax为模拟的持续时间（以秒为单位）。变量dt存储我们希望记录解的时间间隔。在上述代码中，解将在 8,650 秒内每 60 秒记录一次。最终时间通过反复试验选择，以大致对应抛射物的一个轨道。

我们现在准备计算解，这通过调用odeint()函数来完成。解存储在向量xsol中，该向量为每个计算解的时间提供一行。要查看该向量的前几行，我们可以运行以下命令：

xsol[:10]

上述命令产生了以下输出：

array([[  0.00000000e+00,   6.87100000e+06,   8.50000000e+03,
 0.00000000e+00],
 [  5.09624253e+05,   6.85581217e+06,   8.48122162e+03,
 -5.05935282e+02],
 [  1.01700042e+06,   6.81036510e+06,   8.42515172e+03,
 -1.00800330e+03],
 [  1.51991202e+06,   6.73500470e+06,   8.33257580e+03,
 -1.50243025e+03],
 [  2.01620463e+06,   6.63029830e+06,   8.20477026e+03,
 -1.98562415e+03],
 [  2.50381401e+06,   6.49702131e+06,   8.04345585e+03,
 -2.45425372e+03],
 [  2.98079103e+06,   6.33613950e+06,   7.85073707e+03,
 -2.90531389e+03],
 [  3.44532256e+06,   6.14878788e+06,   7.62903174e+03,
 -3.33617549e+03],
 [  3.89574810e+06,   5.93624710e+06,   7.38099497e+03,
 -3.74461812e+03],
 [  4.33057158e+06,   5.69991830e+06,   7.10944180e+03,
 -4.12884566e+03]])

这些值是从0秒到360秒的抛射物的位置和速度向量，时间间隔为60秒。

我们确实希望绘制轨道图。这可以通过在一个单元格中运行以下代码来完成：

subplot(111, aspect='equal')
axis(earth_radius * array([-1.5, 1.5, -1.8, 1.2]))
earth = Circle((0.,0.), 
 earth_radius, 
 ec='Black', fc='Brown', lw=3)
gca().add_artist(earth)
plot(xsol[:,0], xsol[:,1], lw=2, color='DarkBlue')
title('Newton\'s Canon, $v_0={}$ m/s'.format(v0))

我们希望两个轴使用相同的尺度，因为两个轴都表示以米为单位的空间坐标。这在第一行代码中完成。第二行设置轴的限制，使轨道的图表能够舒适地适应图像。

然后，我们使用以下代码绘制一个圆形来表示地球：

earth = Circle((0.,0.), 
 earth_radius, 
 ec='Black', fc='Brown', lw=3)
gca().add_artist(earth)

我们在图表中并未强调使用Artist对象，因为这些对象的级别较低，通常在科学绘图中不需要。这里，我们通过指定圆心、半径和外观选项来构建一个Circle对象：黑色边缘颜色、棕色填充颜色和线宽为 3。第二行代码展示了如何将Circle对象添加到图表中。

绘制地球后，我们使用以下代码绘制轨道：

plot(xsol[:,0], xsol[:,1], lw=2, color='DarkBlue')

这是对plot()函数的标准调用。请注意，我们只绘制xsol数组的前两列，因为这表示抛射物的位置（回想一下，其他两列表示速度）。以下图像是我们得到的输出：

微分方程的数值解是一个复杂的话题，完整的处理超出了本书的范围，但我们将呈现odeint()函数的完整形式并评论其中的一些选项。odeint()函数是 Python 对ODEPACK的lsoda求解器的封装，ODEPACK是一个 Fortran 库。有关求解器的详细信息，请访问people.sc.fsu.edu/~jburkardt/f77_src/odepack/odepack.html

以下代码行是odeint()的完整签名：

odeint(ode_func, x0, tvalues, args=(), Dfun=None, col_deriv=0,
 full_output=0, ml=None, mu=None, rtol=None, atol=None,
 tcrit=None, h0=0.0, hmax=0.0, hmin=0.0,ixpr=0, mxstep=0,
 mxhnil=0, mxordn=12, mxords=5, printmessg=0)

参数ode_func、x0和tvalues已经讨论过。参数args允许我们向所求解的方程传递额外的参数。这是一种非常常见的情况，下面的示例将对此进行说明。在这种情况下，定义系统的函数必须具有以下签名：

ode_func(x, t, p1, p2, … , pn)

这里，p1、p2和pn是额外的参数。这些参数对于单个解是固定的，但可以在不同的解之间发生变化（通常用于表示环境）。传递给args的元组长度必须与ode_func()所需的参数数量完全相等。

以下是常见选项意义的部分列表：

• Dfun是一个计算系统雅可比矩阵的函数。这可能提高解的精度。

• col_deriv指定了雅可比矩阵的导数是沿着列（True，更快）还是行（False）排列。

• 如果full_output设置为True，输出将包含关于解算过程的诊断信息。如果错误积累并且解算过程未成功完成，这可能会很有用。

在本节的最后一个示例中，我们展示了洛伦兹振荡器，这是大气对流的简化模型，也是一个在某些参数值下表现出混沌行为的著名方程。我们还将利用这个示例演示如何在三维空间中绘制解。

洛伦兹系统由以下方程定义：

我们通过定义一个表示该系统的 Python 函数开始，如下所示：

def ode_func(xvec, t, sigma, rho, beta):
 x, y, z = xvec
 return array([sigma * (y - x),
 x * (rho - z) - y,
 x * y - beta * z ])

这个系统与前一个系统的唯一区别在于有了sigma、rho和beta参数。注意，它们只是作为额外的参数添加到ode_func()中。求解这个方程几乎与求解前一个示例一样：

tmax = 50
tdelta = 0.005
tvalues = arange(0, tmax, tdelta) 
ic = array([0.0, 1.0, 0.0])
sol = integrate.odeint(ode_func, ic, tvalues, 
 args=(10., 28., 8./3.))

我们定义了时间数组和初始条件，就像在之前的示例中做的那样。请注意，由于这是一个三维问题，因此初始条件是一个包含三个分量的数组。接下来是对odeint()的调用。此时调用增加了一个额外的参数：

args=(10., 28., 8./3.)

这将sigma、rho和beta分别设置为10、28和8/3。这些值已知对应于混沌解。

然后可以通过以下代码绘制解决方案：

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(8,8))
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')
x, y, z = sol.transpose() 
ax.plot(x, y, z, lw=0.5, color='DarkBlue')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')

前三行代码设置了三维绘图的坐标轴。下一行提取了适合绘图的数据格式：

x, y, z = sol.transpose()

这段代码展示了一个常见的模式。数组sol包含了解决方案的坐标（以列的形式），因此我们将数组转置，使得数据沿数组的行排列，然后将每一行分配给变量x、y和z。

其他几行代码也很简单：我们调用了plot()函数，然后为坐标轴添加标签。以下是我们得到的输出图形：

上面的图像被称为经典的洛伦兹蝴蝶，是一个奇异吸引子的 striking 示例。

使用 Numba 和 NumbaPro 加速计算

本节我们将讨论Numba和NumbaPro，这两个非常激动人心的库，用于加速NumPy代码。Numba和NumbaPro由Continuum Analytics公司创建，后者也生产 Anaconda 发行版。Numba是标准的 Anaconda 发行版的一部分，但NumbaPro是一个商业产品，必须单独购买作为Accelerate包的一部分。不过，NumbaPro可以免费试用一段时间。

这些库的独特之处在于它们允许通过增加几行代码来加速代码执行。作为第一个示例，我们来看看以下几行代码来乘法两个矩阵：

def matrix_multiply(A, B):
 m, n = A.shape
 n, r = B.shape
 C = zeros((m, r), float64)
 for i in range(m):
 for j in range(r):
 acc = 0
 for k in range(n):
 acc += A[i, k] * B[k, j]
 C[i, j] = acc
 return C

上面的代码使用了矩阵乘法的简单定义，看起来就像我们在用 C 语言实现算法时写的代码。这不是 Python 风格的代码，显然也没有经过优化。（在实际情况中，通常会直接使用NumPy的内置矩阵乘法。）特别需要注意的是，矩阵的维度并未进行检查：假设A的列数等于B的行数。

首先让我们尝试使用小矩阵进行计算，如下所示：

A = array([[1,2,0],[1,-3,4],[0,-2,1],[3,7,-4]], dtype=float64)
B = array([[3,4],[-2,0],[2,4]], dtype = float64)
C = matrix_multiply(A, B)
print A
print B
print C

我们首先定义矩阵A和B（注意它们的维度是兼容相乘的）。与本节中的所有示例一样，我们小心地指定了数据类型（这可能有助于优化）。然后，我们简单地调用matrix_multiply，将结果存储在数组C中，并打印这三个矩阵。结果如下所示：

[[  1\.   2\.   0.]
 [  1\.  -3\.   4.]
 [  0\.  -2\.   1.]
 [  3\.   7\.  -4.]]
[[  3    4]
 [ -2    0]
 [  2    4]]
[[ -1\.   4.]
 [  17\. 20.]
 [  6\.   4.]
 [ -13\. -4.]]

你可以通过手动检查几个条目来验证算法的正确性。或者，我们也可以检查结果是否与内置的矩阵乘法一致，如下所示：

C - A.dot(B)
array([[ 0\. ,  0.],
 [ 0\. ,  0.],
 [ 0\. ,  0.],
 [ 0\. ,  0.]])

一切看起来都很好。现在，我们希望定义一些更大的随机矩阵，如下所示：

n = 100
A = rand(n, n)
B = rand(n, n)

在 64 位架构下，前面的代码行会自动生成 64 位浮点数矩阵。接下来，我们将矩阵相乘并计时，如下所示：

%%timeit
C = matrix_multiply(A, B)

上述计算的输出结果如下：

1 loops, best of 3: 472 ms per loop

计时结果当然会因运行代码的机器而异。此示例在运行 Microsoft Windows 7 64 位操作系统的 Intel Core i7 处理器，主频为 3.5 GHz，内存为 16 GB 的计算机上运行。

现在让我们看看如何快速优化这个函数。首先，从Numba模块加载jit函数，如下所示：

from numba import jit

然后，定义带有@jit装饰器的函数，如下所示：

@jit
def matrix_multiply_jit(A, B):
 m, n = A.shape
 n, r = B.shape
 C = zeros((m, r), float64)
 for i in range(m):
 for j in range(r):
 acc = 0
 for k in range(n):
 acc += A[i, k] * B[k, j]
 C[i, j] = acc
 return C

注意，代码中唯一的更改是添加装饰器。（我们还更改了函数的名称以避免混淆，但这并非必要。）装饰器是一个高级的 Python 主题，但我们不需要深入讨论它们的工作原理。有关装饰器的更多信息，请参考 Simeon Franklin 在simeonfranklin.com/blog/2012/jul/1/python-decorators-in-12-steps/上的优秀博客文章。

现在让我们计时我们的代码，如下所示：

%%timeit
C = matrix_multiply_jit(A, B)

以下是结果输出：

1000 loops, best of 3: 1.8 ms per loop

这是一行代码实现的 260 倍改进！在这里应该保持适度的观点，因为这种加速不能期望适用于通用代码。请记住，我们故意编写代码，不使用NumPy中已经优化过的函数。为了比较和充分披露，让我们将其与内置的dot()方法进行比较：

%%timeit
C = A.dot(B)

结果输出如下：

10000 loops, best of 3: 28.6 µs per loop

因此，即使加速了，我们的函数也无法与内置的NumPy竞争。我们再次强调，本节的目标是介绍加速技术的概述，而不是深入探讨复杂的优化方法。

了解@jit装饰器的工作原理是值得的。当调用由@jit装饰的函数时，库会尝试推断参数和返回值的数据类型，并动态生成函数的编译版本，然后调用它。结果是一个类似于用 C 编写的代码的函数调用。

可以指定数据类型，而不是让参数和返回值的类型被推断，这可能会提高性能。以下表格列出了Numba支持的数据类型及其缩写：

数据类型	缩写
boolean	b1
bool_	b1
byte	u1
uint8	u1
uint16	u2
uint32	u4
uint64	u8
char	i1
int8	i1
int16	i2
int32	i4
int64	i8
float_	f4
float32	f4
double	f8
float64	f8
complex64	c8
complex128	c16

所有这些名称都在Numba模块中定义。例如，要定义一个将两个浮点值相加的函数，我们使用以下代码：

from numba import jit, f8
@jit (f8(f8,f8))
def my_sum(a, b):
 return a + b

注意装饰器语法，如下所示：

@jit (f8(f8,f8))

这指定了一个接受两个 float64 参数并返回一个 float64 值的函数。然后，函数按如下方式调用：

my_sum(3.5, 6.9)

这将产生预期的结果。然而，如果我们尝试类似以下代码，我们会遇到错误：

a = array([1.4, 2.0])
b = array([2.3, 5,2])
my_sum(a,b)

然而，使用 @jit 装饰器可以使用数组。要定义一个将两个一维数组相加的函数，可以使用以下代码行：

@jit (f8:)
def vector_sum(a, b):
 return a + b

注意如何指定向量。二维数组表示为 f8[:,:]，三维数组表示为 f8[:,:,:]，以此类推。

NumbaPro 是 Numba 的商业版，并增加了多个增强功能。我们将重点讨论使用图形处理单元（GPU）作为示例，进行并行处理，这是一个激动人心的新技术，现已轻松应用于笔记本中。

要运行以下示例，读者必须拥有 NumbaPro、CUDA 兼容的 GPU（以下简称“设备”）以及最新的 CUDA 兼容驱动程序。

CUDA 兼容设备的列表可以在developer.nvidia.com/cuda-gpus找到。验证您的设备兼容性后，请从developer.nvidia.com/cuda-downloads下载并安装适合平台的最新版本 CUDA SDK。CUDA 工具包附带了一些示例，您可以使用这些示例来测试安装情况。

NumbaPro 下载可以通过store.continuum.io/cshop/accelerate/获得。下载并安装 Accelerate 库。

要测试设置，启动一个 IPython 笔记本，并在单元格中运行以下内容：

import numbapro
numbapro.check_cuda()

如果一切正常，系统将打印出由 Anaconda 安装的 CUDA 库的列表以及系统中可用的 CUDA 兼容设备的列表。您还会看到显示末尾的 PASSED 消息。

尽管 CUDA 编程是通向大规模并行计算的相对简单路径，但在我们开始第一个 CUDA 程序之前，仍然需要掌握一些概念。我们将在此简要概述架构的基础内容，只讨论足够运行后续示例的内容。有关完整规格，请参见CUDA 编程指南，该指南可通过docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#programming-model获得。

GPU 最初是为比计算机的 CPU 更快速地处理渲染操作而设计的。这种处理加速在很大程度上是通过大规模并行化渲染管道所需的图形操作来实现的。

一个CUDA兼容的 GPU 由一组流多处理器（SMs）组成。每个 SM 本身无法与当前的 CPU 在速度上竞争。然而，多个 SM 能够合作解决问题，这一点弥补了这一不足。SM 还可以访问 GPU 中的内存，这部分内存称为设备内存。

在CUDA中，运行在 CPU 中的代码和运行在设备（GPU）中的代码有严格的分离。一个特别的限制是，CPU 代码只能访问常规计算机内存，而设备代码只能访问设备内存。运行在设备中的代码被指定在一个名为内核的函数中。内核被编译成 SM 可以理解的低级语言，并异步地运行在每个 SM 中（意味着每个 SM 按照自己的节奏运行，除非遇到特殊的同步指令）。因此，CUDA中的简单计算通常需要以下三个步骤：

1. 将输入数据从计算机内存传输到设备内存。

2. 在设备中启动内核。

3. 将设备内存中的输出数据传输到计算机内存，以便 CPU 能够再次访问它。

正如你将看到的，Python CUDA使内存传输变得透明。（如果需要，仍然可以通过编程方式控制它。）

内核同时在一组 SM 中启动，每个线程独立进行计算。每个 SM 可以并行运行多个线程，并且可以访问设备的所有内存。（架构更为复杂，还有其他类型的内存可用，这里不作讨论。）在最简单的情况下，每个线程只访问几个内存区域，每个区域包含一个 64 位浮动值，而且线程访问的内存不会被其他线程访问。因此，无需同步。在更复杂的问题中，同步可能成为一个主要问题。

正在运行的线程集合具有二级数组结构：

• 线程按块组织。每个块是一个最多具有3个维度的线程数组。块的维度存储在一个名为blockDim的变量中。块中的线程由变量threadIdx标识。这是一个具有三个整数字段的结构：threadIdx.x、threadIdx.y和threadIdx.z。这些字段唯一标识块中的每个线程。

• 块按网格组织。网格是一个最多具有3个维度的块数组。网格的维度存储在一个名为gridDim的变量中。网格中的块由变量gridIdx标识。这是一个具有三个整数字段的结构：gridIdx.x、gridIdx.y和gridIdx.z。这些字段唯一标识网格中的每个块。

下面的图给出了这种组织结构的一个示例：

在前面的示例中，gridDim是(2, 3, 1)，因为有两行三列的块（以及一个单独的空间维度）。网格中的所有块都是一维的，所以blockDim是(4, 1, 1)。例如，底部行中第一个块中的第三个线程由以下代码行标识：

blockIdx.x=0, blockIdx.y=1, blockIdx.z=1
threadIdx.x=2, threadIdx.y=1, threadIdx.z=1

在运行时，每个线程都可以访问这些标识信息。

CUDA架构的一个关键点如下：

• 同一块中的所有线程始终在一个 SM 中并行执行，直到块中的所有线程都执行完毕

• 不同的块可以根据 SM 的可用性并行或串行执行计算

我们现在准备使用 Python CUDA定义内核。我们将编写一个函数，在 GPU 中计算两个向量的和。在一个单元格中运行以下代码：

from numbapro import cuda
@cuda.jit('void(float64[:], float64[:], float64[:])')
def sum(a, b, result):
 i = cuda.threadIdx.x 
 result[i] = a[i] + b[i]

我们假设只有一个线程块，每个线程负责在单个位置添加数组的元素。线程负责的数组位置由threadIdx.x的值确定。请注意，内核没有返回值。我们需要指定一个数组result，用于存储计算的返回值。

现在让我们看看这个函数是如何被调用的。请注意，网格和块的几何形状并没有在内核中定义。（如果需要，内核可以获取几何信息；稍后会详细讲解。）这在启动内核时完成：

a = array([1,2,3,4,5], dtype=float64)
b = array([5,4,3,2,1], dtype=float64)
result = array([0,0,0,0,0], dtype=float64)
sum1,5
print result

上面的代码行给出了以下输出：

[ 6\.  6\.  6\.  6\.  6.]

这段代码的关键点在于以下这一行：

sum1,5

上面的代码行在一个包含1个块、每个块包含5个线程的网格中启动内核。网格和块都是一维的。现在让我们添加更大的向量：

n = 64
a = arange(0,n,dtype=float64)
b = arange(n,0,-1, dtype=float64)
result = zeros(n,dtype=float64)
sum1,n
print result[:5]

上面的代码行本质上与之前的代码相同，只是稍微通用了一些，因为数组的大小可以变化。我们想要做的是增加n的大小。如果你尝试使用n=10000这样的值，将会出现CUDA_ERROR_INVALID_VALUE类型的错误。问题在于单个 SM 可以运行的线程数是有限的，也就是说，单个块中可以执行的线程数有限。为了处理更大的向量，我们需要修改代码，使其能够处理多个块。为此，请以以下方式修改sum()函数的定义：

from numbapro import cuda
@cuda.jit('void(float64[:], float64[:], float64[:], int32)')
def sum(a, b, result, n):
 tx = cuda.threadIdx.x
 bx = cuda.blockIdx.x
 bsz = cuda.blockDim.x
 i = tx + bx * bsz
 if i < n:
 result[i] = a[i] + b[i]

首先需要注意的是，我们包含了一个类型为int32的参数，用于保存正在加法的数组的大小。现在的重点是，位于不同块中的线程必须访问不同的内存区域，因此计算与线程相关联的索引i变得更加复杂。本质上，我们必须知道当前块之前有多少个块，然后将其乘以块的维度，并加上当前线程的索引。然后，在访问相关的内存位置之前，我们会检查索引是否有效。这可以防止线程访问不属于输入/输出数组的区域，并且是更复杂代码中一个至关重要的检查。要测试代码，请运行以下内容：

n = 100000
a = arange(0,n,dtype=float64)
b = arange(n,0,-1, dtype=float64)
result = zeros(n,dtype=float64)
sum1000,64
print result[:5]

上述代码应该可以顺利运行。请注意，我们指定了一个包含1000个块且每个块有64个线程的网格。网格中的块数是无限制的，设备负责以最佳方式分配 SM。请注意，块数必须足够大，以覆盖输入/输出数组。在我们的例子中，这意味着blockDim.x * gridDim.x >= n。

现在我们可以开始使用大向量进行计算。试试下面的代码：

n = 100000
a = rand(n)
b = rand(n)
result = zeros(n, dtype=float64)
bd = 10000
gd = 64
if bd * gd < n:
 print 'Block/grid dimensions too small'
else:
 sumbd,gd
print result[:10]

读者应尝试不同的n、bd和gd值。请记住，gd的最大值取决于你电脑中的设备。一个有趣的实验是检查计算在更大n值下的扩展性。

概述

本章我们介绍了在SciPy中使用高级数学算法，包括解方程、寻找最优值、积分和微分方程。本章最后讨论了如何利用 GPU 并行化来加速计算。

附录 A. IPython 笔记本参考卡

启动笔记本

要启动笔记本，请打开终端窗口并运行以下命令：

ipython notebook

在运行上述命令时，确保你处于包含笔记本的目录中。

键盘快捷键

一些重要的键盘快捷键如下：

• 要进入编辑模式，按下Enter键或点击单元格

• 要进入命令模式，按Esc

编辑模式中的快捷键

编辑模式中使用的一些重要快捷键如下：

• 要运行一个单元格，使用以下快捷键：

  • 要运行一个单元格并跳转到下一个单元格，按Shift + Enter

  • 要运行一个单元格但保持在同一单元格中，按Ctrl + Enter

  • 要运行一个单元格并在其下方插入一个新单元格，按Alt + Enter

  • 要在当前单元格中创建一个新行，按Enter

• 要缩进内容，按Tab

• 要启动代码补全，开始在单元格中输入然后按Tab

• 要全选，按Ctrl + A

• 要撤销一个操作，按Ctrl + Z

• 要重做一个操作，按Ctrl + Y或Ctrl + Shift + Z

• 要跳转到单元格的开头，按Ctrl + Home

• 要跳转到单元格的末尾，按Ctrl + End

• 要拆分一个单元格，按Ctrl + Shift + -

命令模式中的快捷键

• 要列出键盘快捷键，按H

• 要将单元格模式更改为以下之一，快捷键如下：

  • 代码：按Y

  • Markdown：按M

  • 标题：按 1 到 6 之间的数字，选择标题的级别

  • 原始 NBConvert：按R

• 要选择当前单元格上方的单元格，按上箭头或K

• 要选择当前单元格下方的单元格，按下箭头或J

• 要将一个单元格向上移动一个位置，按Ctrl + K

• 要将一个单元格向下移动一个位置，按Ctrl + J

• 要在当前单元格上方插入一个新单元格，按A

• 要在当前单元格下方插入一个新单元格，按B

• 要剪切一个单元格，按X

• 要复制一个单元格，按C

• 要将一个单元格粘贴到当前单元格下方，按V

• 要将一个单元格粘贴到当前单元格上方，按Shift + V

• 要删除一个单元格，按D

• 要撤销删除操作，按Z

• 要将当前单元格与下方的单元格合并，按Shift + M

• 要切换行号显示，按L

导入模块

加载一些重要模块的步骤如下：

• 要加载NumPy和 matplotlib 以进行交互式工作，并支持内联图形，请运行以下命令：

  pylab inline
  
  

• 要加载NumPy和 matplotlib 而不将名称导入当前命名空间，并且支持内联图形，请运行以下命令行：

  pylab -–no-import-all inline
  
  

• 要加载SciPy模块，使用以下任何标准的 Python 导入命令：

  import scipy.<module>
  import scipy.<module> as <local-name>
  from scipy.<module> import <function>
  
  

如果使用了–no-import-all选项，函数和对象必须加上适当的模块名称前缀，如下所示：

• 对于NumPy函数和对象，使用numpy或np。

• 对于交互式图形，使用pyplot或plt。

系统中安装的库模块以及用户创建的.py扩展名的模块，可以通过标准的 Python 机制进行导入。

获取帮助

获取帮助的方式有很多种：

• 要启动交互式帮助，可以运行以下命令：

  help()
  
  

• 要获取函数或对象的帮助，可以运行以下命令：

  help(object)
  help(function)
  <object>?
  <function>?
  
  

• 对于自动补全，开始输入函数/对象/方法的名称，然后按Tab。

• 要获取工具提示，开始输入函数/对象/方法的名称，然后按Shift + Tab。

附录 B. Python 简要回顾

介绍

本附录将简要介绍 Python 语法。这并不是一个 Python 编程课程，而是供不熟悉该语言的读者作为快速入门使用。以下主题将在本附录中涵盖：

• 基本类型、表达式、变量及其赋值

• 序列类型

• 字典

• 控制结构

• 函数、对象和方法

基本类型、表达式、变量及其赋值

任何可以在 Python 代码中引用的数据都被视为对象。对象用于表示从原子数据（如数字）到非常复杂的数据结构（如多维数组、数据库连接和各种格式的文档）。

在对象层次结构的根部是数值数据类型。这些包括以下内容：

• 整数：Python 中有三种类型的整数。

  • 普通整数：它们在本地架构中表示，大多数系统中通常是 32 位或 64 位的有符号值。

  • 长整型：它们是具有无限范围的整数，取决于可用内存。大多数情况下，程序员不需要关心普通整数和长整型之间的区别。Python 会透明地处理这两种类型之间的转换。

  • 布尔值：它们表示False和True两个值。在大多数情况下，它们分别等价于0和1。

• 浮点数：它们表示本地的双精度浮点数。

• 复数：它们表示复数，表示为一对双精度浮点数。

下表展示了每种数据类型的字面量（即常量）示例：

数据类型	字面量
整数	0, 2, 4, …, 43882838388``5L, 5l（长整型）0xFE4（十六进制）03241（八进制）
实数（float）	5.34, 1.2, 3., 0``1.4e-32（科学记数法）
复数	1.0+3.4j, 1+2j, 1j, 0j, complex(4.3, 2.5)

虚数单位由j表示，但只有在它跟随数字字面量时（否则它表示名为j的变量）。因此，要表示虚数单位，必须使用1j，复数零为0j。复数的实部和虚部总是以双精度浮点值存储。

注意

请注意，NumPy大大扩展了数值类型的集合，以便进行高效的数值计算。

赋值语句用于将值存储到变量中，如下所示：

a = 3
b = 2.28
c = 12
d = 1+2j

Python 支持多重同时赋值，因此前四行代码可以等效地写成一行，如下所示：

a, b, c, d = 3, 2.28, 12, 1+2j

在多重赋值中，右侧的所有表达式都会在赋值之前被求值。例如，交换两个变量的值的常见惯用法如下所示：

v, w = w, v

作为练习，读者可以尝试预测以下语句的结果，前提是已知前面的变量赋值：

a, b, c = a + b, c + d, a * b * c * d
print a, b, c, d

以下示例展示了如何计算二次方程的两个解：

a, b, c = 2., -1., -4.
x1, x2 = .5 * (-b - (b ** 2 - 4 * a * c) ** 0.5), .5 * (-b + (b ** 2 - 4 * a * c) ** 0.5)
print x1, x2

请注意，我们通过使用小数点强制将变量 a、b 和 c 转换为浮点值。在进行数值计算时，这是一种良好的实践。下表包含了 Python 运算符的部分列表：

运算符	Python 运算符
算术	+（加法）-（减法，取负）*（乘法）/（除法，见表格下方的注释）//（整数除法）%（余数）
比较	==（等于）>（大于）<（小于）>=（大于或等于）<=（小于或等于）!=（不等于）
布尔	and or not
位运算布尔	&（与）&#124;（或）^（异或）~（非）
位移	<<（左移）>>（右移）

注释

使用除法运算符（/）时需要小心。如果操作数是整数，则此操作的结果是整数商。例如，34/12 结果为 2。要获得浮点结果，我们必须输入浮点操作数，例如 34./12.，或者添加以下语句：

from __future__ import division

// 运算符始终表示整数除法。

算术运算符遵循运算顺序规则，括号可以改变该顺序。比较运算符的优先级低于算术运算符，而 or、and 和 not 运算符的优先级更低。因此，像以下的表达式会得到预期的结果：

2 + 3 < 5 ** 2 and 4 * 3 != 13

换句话说，前面的命令行被解析如下：

(((2 + 3) < (5 ** 2)) and ((4 * 3) != 13))

逻辑运算符 and 和 or 采用短路求值，因此，例如，命令中的第二个比较不会被评估：

2 < 3 or 4 > 5

位运算符和移位运算符的优先级规则可能不像直观那样清晰，因此建议始终使用括号来指定运算顺序，这也有助于代码的清晰度。

Python 还支持增强赋值。例如，以下命令行首先将值 5 赋给 a，然后将 a 的值加一：

a = 5
a += 1

注释

Python 不支持递增/递减运算符，如 C 语言中的 a++ 和 ++a。

所有 Python 运算符都有对应的增强赋值语句。任何运算符 $ 的一般语义如下所示：

v $= <expression>

前面的语句等价于以下：

v = v $ (<expression>)

注释

请注意，$ 不是一个有效的 Python 运算符，它仅作为一个占位符表示通用运算符。

序列类型

Python 序列类型用于表示有序的对象集合。它们分为可变序列类型和不可变序列类型。在这里，我们只讨论lists（可变）和tuples与strings（均为不可变）。其他序列类型将在本节末尾提到。

列表

以下示例演示如何在 Python 中构造一个列表并将其赋值给变量：

numbers = [0, 1.2, 234259399992, 4+3j]

列表中的个别条目可以通过索引表示法访问，如下所示：

numbers[2]

请注意，索引总是从0开始。允许使用负索引，它们表示从列表末尾开始的位置。例如，numbers[-1]是最后一个元素，numbers[-2]是倒数第二个元素，以此类推。

由于列表是可变的序列类型，我们可以就地修改其中的条目：

numbers[0] = -3
numbers[2] += numbers[0]
print numbers

另一种重要的引用 Python 序列类型元素的方法是切片，它允许从列表中提取子列表。由于这个话题对NumPy数组非常重要，我们将讨论延迟到附录 C，NumPy 数组。

Python 列表具有一组非常实用的功能，其中一些在以下代码示例中得到了体现：

• 要查找列表的长度，请使用以下命令：

  len(numbers)
  
  

• 要就地反转列表，请使用以下命令：

  numbers.reverse()
  print numbers
  
  

• 要添加一个新元素，请使用以下命令：

  numbers.append(35)
  print numbers
  
  

• 要就地对列表进行排序，请使用以下命令：

  values = [1.2, 0.5, -3.4, 12.6, 3.5]
  values.sort()
  print values
  values.sort(reverse=True)
  print values
  
  

• 要在指定位置插入一个值，请使用以下命令：

  values.insert(3, 6.8)
  print values
  
  

• 要扩展列表，请使用以下命令：

  values.extend([7,8,9])
  print values
  
  

Python 有一些便捷的方式来构造常用的列表。range()函数返回一个等间隔的整数列表。最简单的形式返回一个从0开始的连续整数列表：

range(10)

上述命令返回以下列表：

[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

请注意，最后一个元素比函数调用中给出的参数少 1。通常规则是，range(n)返回一个包含n个元素的列表，从零开始，因此最后一个元素是n-1。要从非零值开始，请使用以下带有两个参数的版本：

range(3, 17)

第三个参数指定一个增量。以下命令行会生成一个包含所有小于 100 的正 6 的倍数的列表：

range(6,100,6)

也可以使用负增量：

range(20, 2, -3)

列表支持连接，表示为+运算符：

l1 = range(1, 10)
l2 = range(10, 0, -1)
l3 = l1 + l2
print l3

注意

请注意，对于NumPy数组，+运算符被重新定义为表示向量/矩阵加法。

乘法运算符（*）可用于通过重复给定列表的元素来构造一个新列表，如下所示：

l4 = 3*[4,-1,5]
print l4

在 Python 中构造列表的最灵活方法是使用列表推导式。完整的讨论超出了本附录的范围，但以下示例说明了其中的一些可能性：

• 要显示从0到10（包括10）的整数的平方列表，请使用以下命令行：

  [n ** 2 for n in range(11)]
  
  

• 要显示一个整数的所有除数列表，请使用以下命令行：

  k = 60
  [d for d in range(1, k+1) if k % d == 0]
  
  

• 要显示小于 100 的所有素数列表，请使用以下命令行（效率较低）：

  [k for k in range(2,101) if len([d for d in range(1, k+1) if k % d == 0])==2]
  
  

• 要显示具有整数坐标的点的元组列表及其与原点的距离，请使用以下命令行：

  [(i,j,(i*i+j*j)**0.5) for i in range(5) for j in range(6)]
  
  

元组

元组与列表类似，但它们是不可变的——一旦创建，其元素就不能修改。

会发生变化。以下命令行将导致错误信息：

t1 = (2,3,5,7)
t1[2] = -4

元组在 Python 中有一些特殊用途。它们可以作为字典中的索引（因为它们是不可变的）。它们还构成了 Python 用来从函数返回多个值的机制。例如，内建函数 divmod() 返回一个元组，其中包含整数商和余数：

divmod(213, 43)

元组支持与列表相同的序列接口，除了会修改元组的那些方法。例如，没有名为 sort() 的方法可以就地对元组进行排序。

字符串

Python 字符串表示一个不可变的字符序列。有两种字符串类型：str，表示 ASCII 字符串，和 unicode，表示 Unicode 字符串。

字符串文字是由单引号或双引号括起来的字符序列，如下所示：

s1 = 'I am a string'
s2 = "I am a string"
print s1
print s2

单引号和双引号之间没有语义上的区别，唯一的不同是单引号字符串可以包含双引号，而双引号字符串可以包含单引号。例如，以下命令行是正确的：

s3 = "I'm a string"
print s3

字符串有两个主要用途：作为字典的索引和打印信息。当打印信息时，字符串有一个 format() 方法，可以轻松地显示信息。我们经常使用这个功能为图形添加注释。这里有一个例子：

n = 3
message = 'The square root of {:d} is approximately {:8.5f}.'.format(n, n ** 0.5)
print message

在前面的例子中，有两个格式说明符：

• {:d}：这指定了整数值的十进制格式

• {:8.5f}：这指定了一个宽度为 8，小数点后有 5 位的浮点数格式

格式说明符会按顺序与参数匹配，在这个例子中是 n 和 n ** 0.5。

字符串有丰富的接口。如果你需要编写与字符串相关的代码，很可能有一个内建函数可以实现，且只需要很少的修改。所有可用字符串方法的列表以及格式化功能，可以参考 docs.python.org/2/library/stdtypes.html#string-methods。

字典

Python 字典是一种包含键值对的数据结构。键必须是不可变类型，通常是字符串或元组。以下是一个构建字典的示例：

grades = {'Pete':87, 'Annie':92, 'Jodi':78}

要访问某个项，我们提供键作为索引，如下所示：

print grades['Annie']

字典是可变的，所以我们可以使用它们来更改条目的值。如果 Jodi 做了额外的工作来提高她的成绩，我们可以这样更改：

grades['Jodi'] += 10
print grades['Jodi']

要向字典添加条目，只需为新键赋值：

grades['Ivan']=94

然而，尝试访问一个不存在的键会导致错误。

一个需要注意的重要点是，字典是无序的。以下代码是迭代字典的标准惯用法：

for key, item in grades.iteritems():
 print "{:s}'s grade in the test is {:d}".format(key, item)

这里的关键点是，输出与字典中条目的添加顺序完全无关。

如需了解更多字典接口的细节，可以参考docs.python.org/2/library/stdtypes.html#mapping-types-dict。

控制结构

控制结构允许改变代码执行的流程。我们关注的有两种结构：分支和循环。

分支结构根据测试结果执行不同的代码。以下示例展示了一个改进版的代码，用于求解二次方程。使用了if-then-else结构来处理实数解和虚数解的情况，如下所示：

a, b, c = 2., -4., 5.
discr = b ** 2 - 4 * a * c
if discr >= 0:
 sqroot = discr ** 0.5
 x1 = 0.5 * (-b + sqroot)
 x2 = 0.5 * (-b - sqroot)
else:
 sqroot = (-discr) ** 0.5
 x1 = 0.5 * (-b + sqroot * 1j)
 x2 = 0.5 * (-b - sqroot * 1j)
print x1, x2

上面的代码首先计算二次方程的判别式。然后，使用if-then-else语句来判断根是实数还是虚数，这取决于判别式的符号。注意代码的缩进，Python 通过缩进来定义语句块的边界。if-then-else结构的一般形式如下：

if <condition>:
 <statement block T>
else:
 <statement block F>

首先，条件<condition>被评估。如果为True，则执行语句<statement block T>；否则，执行语句<statement block F>。else:子句可以省略。

Python 中最常见的循环结构是for语句。下面是一个示例：

numbers = [2, 4, 3, 6, 2, 1, 5, 10]
for n in numbers:
 r = n % 2
 if r == 0:
 print 'The integer {:d} is even'.format(n)
 else:
 print 'The integer {:d} is odd'.format(n)

我们首先定义一个整数列表。for语句使得变量n依次取列表中的每个值，并对每个值执行缩进块的代码。注意，for循环内部有一个if-then-else结构。此外，print语句是双重缩进的。

for循环常用于执行简单的搜索。一个常见的场景是，当满足某个条件时需要跳出循环。以下代码查找整数范围内的第一个完全平方数：

for n in range(30, 90):
 if int(n ** 0.5) ** 2 == n:
 print n
 break

对于给定范围内的每个n值，我们先计算n的平方根，取其整数部分，再计算平方。如果结果等于n，则进入if块，打印n，然后跳出循环。

如果在范围内没有完全平方数怎么办？将前面的函数range(30, 60)改为range(125, 140)。运行命令行时不会输出任何内容，因为在125到140之间没有完全平方数。现在，将命令行改为如下：

for n in range(125, 140):
 if int(n ** 0.5) ** 2 == n:
 print n
 break
else:
 print 'There are no perfect squares in the range'

else子句只有在执行没有跳出循环时才会被执行，此时会打印出相应的消息。

另一个常见的情况是某些值在迭代过程中需要跳过。以下示例中，我们打印一系列-1到1之间随机数的平方根，但只有当这些数为正时才打印：

import random
numbers = [-1 + 2 * rand() for _ in range(20)]
for n in numbers:
 if n < 0:
 continue
 print 'The square root of {:8.6} is {:8.6}'.format(n, n ** 0.5)

当 Python 遇到continue语句时，它会跳过当前执行块的剩余部分，并继续执行控制变量的下一个值。

另一个经常使用的控制结构是while循环。该结构在条件为真时执行一组命令。例如，假设我们想计算一组随机生成的值的累加和，但仅当和超过某个值时才停止。可以使用以下代码来实现：

import random
bound = 10.
acc = 0.
n = 0
while acc < bound:
 v = random.random()
 acc += v
 print 'v={:5.4}, acc={:6.4}'.format(v, acc)

另一种比预期更常见的情况需要一种叫做永远循环的模式。这种情况发生在循环开始时需要检查的条件不可用时。例如，下面的代码实现了著名的3n+1游戏：

n = 7
while True:
 if n % 2 == 0:
 n /= 2
 else:
 n = 3 * n + 1
 print n
 if n == 1:
 break

游戏从一个任意整数开始，这里是7。然后，在每次迭代中，我们测试n是否为偶数。如果是，我们将其除以2；否则，将其乘以3并加上1。然后，我们检查是否达到了1。如果是，我们退出循环。由于我们不知道是否需要在循环结束前退出，因此我们使用一个永远循环，如下所示：

while True:
 <statements>
 if <condition>:
 break
 <possibly more statements>

一些程序员避免使用这种结构，因为如果不小心，很容易导致无限循环。然而，它在某些情况下确实非常有用。顺便说一下，3n+1问题中的循环是否对所有初始值都会停止还是一个开放性问题！读者可以尝试使用初始值n=27来感受一下。

函数、对象和方法

现在我们来介绍使 Python 如此灵活和强大的构造，它的面向对象特性。我们在前面的章节中已经看到了一些面向对象的代码示例（面向对象的范式是 Python 的核心，很难写出不使用它的代码），但现在我们将对这些特性进行更具体的讨论。

函数

我们已经看到许多使用函数的例子。例如，len()函数用于计算列表的长度：

lst = range(1000)
print len(lst)

调用函数的最基本语法如下：

function_name(arg1, arg2, …, argn)

在这种情况下，arg1、arg2、…、argn被称为位置参数，因为它们是根据出现的位置进行匹配的。例如，我们考虑内置函数pow()。这个函数最多接受三个参数：

pow(b, n, m)

在这种形式下，前面的函数使用优化算法来计算b的n次方模m。（如果你想知道，这是公共密钥密码学中的一个重要操作。）参数b、n和m是根据它们的位置来关联的。例如，要计算12的十次方模15，我们使用以下命令：

pow(12, 10, 15)

Python 还支持任意大小的参数序列。例如，max()函数计算任意序列中的最大值：

max(2,6,8,-3,3,4)

前面的命令返回值8。

传递参数给函数的第三种方式是使用关键字参数。这非常有用，因为通常很难记住位置参数的准确顺序。（例如，我不太愿意编写超过三个或四个位置参数的函数。）

例如，内置的int()函数可以用于将字符串转换为整数。可选的关键字参数base让我们可以指定转换的进制。例如，以下命令行将给定的2进制整数赋值给n：

n = int('10111010100001', base=2)
print n

关键字参数总是有默认值。在我们的示例中，如果没有指定基数，默认假设为10。

我们经常需要编写自己的函数。这是通过关键字def来实现的。作为示例，假设我们要编写代码来实现著名的bisection方法以数值求解方程。一个可能的解决方案如下：

def bisection(f, a, b, tol=1e-5, itermax=1000):
 fa = f(a)
 fb = f(b)
 if fa * fb > 0:
 raise ValueError('f(a) and f(b) must have opposite signs')
 niter = 0
 while abs(a-b) > tol and niter < itermax:
 m = 0.5 * (a + b)
 fm = f(m)
 if fm * fa < 0:
 b, fb = m, fm
 else:
 a, fa = m, fm
 return min(a, b), max(a, b)

上述函数需要三个重要且必要的参数：

• f函数接受一个浮点值作为输入，并返回一个浮点值作为输出。

• 浮点值a和b，它们指定一个包含函数零点的区间。

另外两个参数是可选的。参数tol指定结果所需的容差，itermax指定最大迭代次数。要使用bisection()函数，必须先定义函数f。我们将借此机会展示另一种定义 Python 函数的方式，如下：

from math import cos, pi
f = lambda x: cos(x) - x

我们现在准备通过以下命令来调用函数：

bisection(f, 0, pi/2)

上述函数返回以下输出：

(0.7390851262506977, 0.7390911183631504)

请注意，我们设计了该函数返回一个包含零的区间。该区间的长度小于tol，除非达到了最大迭代次数。如果我们希望更小的容差，可以使用以下函数：

bisection(f, 0, pi/2, tol=1E-10)

现在，假设我们关注计算所花费的时间。我们可以通过以下方式限制最大次数：

bisection(f, 0, pi/2, itermax=10, tol=1E-20)

请注意，关键字参数的顺序是无关紧要的，并且在前面的示例中没有达到所需的容差。

对象和方法

对象是 Python 中最通用的数据抽象。实际上，从程序员的角度来看，Python 中的一切都是对象。

对象不过是结构化数据的集合，并且有一个操作这些数据的接口。对象是通过class构造定义的，但我们这里的目标并不是展示如何定义类。尽管设计新类是一个高级话题，但使用现有类却相当简单。

作为示例，让我们探索内置类型str。首先，我们定义一个可以操作的str对象，如下所示：

message = 'Mathematics is the queen of science'

首先，让我们将消息转换为大写，如下所示：

message.upper()

我们说前述语句调用了message对象的upper()方法。方法仅仅是与对象关联的函数。以下是str对象的其他一些方法：

• 要查找子字符串的第一次出现（如果未找到该字符串，返回-1），请使用以下命令行：

  message.find('queen')
  
  

• 要将字符串拆分为单词，请使用以下命令行：

  words = message.split()
  print words
  
  

• 要计算s子字符串出现的次数，请使用以下命令行：

  message.count('e')
  
  

• 要用其他内容替换子字符串，请使用以下命令行：

  message.replace('Mathematics', 'Mme. Curie')
  
  

注意

请注意，前述方法不会改变原始字符串对象，而是返回新的修改过的字符串。字符串是不可变的。对于可变对象，方法可以自由地改变对象中的数据。

总结

在本附录中，我们概述了 Python 语法和特性，涵盖了基本类型、表达式、变量和赋值、基本数据结构、函数、对象和方法。

附录 C. NumPy 数组

介绍

数组是NumPy引入的基本数据结构，它们是我们在本书中讨论的所有科学计算和数据分析库的基础。本附录将简要概述以下数组特性：

• 数组创建和成员访问

• 索引和切片

数组创建和成员访问

NumPy数组是ndarray类的对象，该类表示一个固定大小的多维同质数据集合。

在这里，我们假设已经通过以下命令行导入了NumPy库：

import numpy as np

一旦完成这些步骤，我们可以从一个列表的列表创建ndarray（以后简称为数组对象或简单称为数组），如以下命令行所示：

a = np.array([[-2,3,-4,0],[2,-7,0,0],[3,-4,2,1]],dtype=np.float64)
print a

与 Python 列表和元组不同，数组对象的所有元素必须是相同类型的。这些类型由NumPy对象表示，并称为数组的dtype（数据类型）。在前面的示例中，我们明确指定dtype为float64，它表示一个 64 位浮动小数值。

数组具有几个属性，用于提供关于数据布局的信息。最常用的属性如下：

• 数组的形状可以通过以下命令计算：

  a.shape
  
  

  前面的命令返回元组(3, 4)，因为这是一个具有三行四列的二维数组。有些人可能会感到惊讶，shape属性并非只读，我们可以利用它来重塑数组：

  a.shape = (6,2)
  print a
  
  

  运行前面的示例后，运行a.shape(3,4)可以返回原始维度。

• 数组的维度数量可以通过以下命令获取：

  a.ndim
  
  

  当然，这将返回2。在NumPy中，一个重要的概念是数组的轴。二维数组有两个轴，编号为 0 和 1。如果我们把数组看作是一个数学矩阵，那么轴 0 是垂直的，指向下方，轴 1 是水平的，指向右侧。某些数组方法有一个可选的axis关键字参数，允许用户指定在哪个轴上执行操作。

• 要获取数组中元素的个数，可以使用以下命令：

  a.size
  
  

  在前面的示例中，返回的输出是12，正如预期的那样。

• 数组的一个最终属性是计算数组的转置。这可以通过以下命令完成：

  b = a.T
  print b
  
  

  这创建的一个重要内容是数组a的视图。NumPy包被设计成能高效处理非常大的数组，并且在大多数情况下，除非绝对必要，或明确指示，否则避免复制数据。

• 运行以下代码行：

  print a
  b[1,2] = 11
  print a
  
  

  请注意，数组a的条目2, 1已被更改，这表明变量a和b都指向内存中的同一位置。

• 可以通过empty()函数如下创建一个包含未初始化数据的数组：

  c = np.empty(shape=(3,2), dtype=np.float64)
  print c
  
  

• 使用未初始化的数据是不推荐的，因此最好使用 zeros() 或 ones() 函数，如下所示：

  • 要使用 zeros() 函数，执行以下命令：

    d = np.zeros(shape=(3,2), dtype=np.float64)
    print d
    
    

  • 要使用 ones() 函数，执行以下命令：

    e = np.ones(shape=(3,2), dtype=np.float64)
    print e
    
    

    a_like = np.zeros_like(a)
    print a_like
    
    

  也有一些函数可以创建具有与现有数组相同形状和数据类型的新数组：

• ones_like() 和 empty_like() 函数生成与给定数组相同形状的全 1 数组和未初始化的数据数组。

• NumPy 还有一个 eye() 函数，可以返回给定维度和 dtype 的单位矩阵：

  f = np.eye(5, dtype=np.float64)
  print f
  
  

  行和列的数量不必相同。在这种情况下，生成的矩阵将仅是左单位矩阵或右单位矩阵，具体取决于情况：

  g = np.eye(5, 3, dtype=np.float64)
  print g
  
  

• 数组也可以从现有数据创建。copy() 函数可以如下面这样克隆数组：

  aa = np.copy(a)
  print a
  print aa
  
  

• frombuffer() 函数从暴露（单维）缓冲区接口的对象创建数组。以下是一个示例：

  ar = np.arange(0.0, 1.0, 0.1, dtype=np.float64)
  v = np.frombuffer(ar)
  v.shape = (2, 5)
  print v
  
  

  arange() 函数是 NumPy 对 Python range 的扩展。它的语法类似，但允许浮动范围值。

• loadtxt() 函数从文本文件中读取数组。假设文本文件 matrix.txt 包含以下数据：

   1.3  4.6  7.8
  -3.6  0.4  3.54
   2.4  1.7  4.5
  
  

  然后，我们可以使用以下命令读取数据：

  h = np.loadtxt('matrix.txt', dtype=np.float64)
  print h
  
  

• 数组也可以以 .npy 格式保存和加载：

  np.save('matrix.npy',h)
  hh = np.load('matrix.npy')
  print hh
  
  

索引和切片

为了说明索引，我们先使用以下命令创建一个包含随机数据的数组：

import numpy.random
a = np.random.rand(6,5)
print a

这会创建一个形状为 (6,5) 的数组，包含随机数据。数组的各个元素可以通过常规的索引表示法来访问，例如，a[2,4]。

操作 NumPy 数据的一个重要技巧是使用 切片。切片可以被认为是数组的一个子数组。例如，假设我们想要提取数组 a 中的中间两行和前两列的子数组。考虑以下命令：

b = a[2:4,0:2]
print b

现在，让我们做一个非常重要的观察。切片仅仅是数组的一种视图，并没有实际复制数据。可以通过运行以下命令来验证：

b[0,0]=0
print a

所以，b 的更改会影响数组 a！如果我们确实需要一个副本，我们需要明确表示我们想要一个副本。可以使用以下命令行来实现：

c = np.copy(a[2:4,0:2])
c[0,0] = -1
print a

在切片表示法 i:j 中，我们可以省略 i 或 j，在这种情况下，切片将表示对应轴的开始或结束：

print a[:4,3:]

省略 i 和 j 表示整个轴：

print a[:,2:4]

最后，我们可以使用 i:j:k 的表示法来指定切片中的步幅 k。在以下示例中，我们首先创建一个更大的随机数组来说明这一点：

a = np.random.rand(10,6)
print a
print
print a[1:7:2,5:0:-3]

现在，让我们考虑更高维度数组的切片。我们将通过创建一个非常大的三维数组来开始：

d1, d2, d3 = 4, 5, 3
a = np.random.rand(d1, d2, d3)
print a

假设我们想提取最后一个轴中索引为 1 的所有元素。这可以通过使用省略号对象轻松实现，示例如下：

print a[...,1]

上述命令行等同于以下命令：

print a[:,:,1]

在切片时，也可以沿着某一轴扩展矩阵，如下所示：

print a[0, :, np.newaxis, 0]

将前一个命令行的输出与以下输出进行比较：

print a[0, :, 0]