Python3-科学计算教程-二-

Python3 科学计算教程（二）

原文：Scientific Computing with Python 3

协议：CC BY-NC-SA 4.0

五、高级数组概念

在本章中，我们将解释数组的一些更高级的方面。首先，我们将介绍数组视图的概念，然后是布尔数组以及如何比较数组。我们简要描述索引和矢量化，解释稀疏数组，以及广播等一些特殊主题。

数组视图和副本

为了精确控制内存的使用方式，NumPy 提供了数组视图的概念。视图是较小的数组，与较大的数组共享相同的数据。这就像引用一个单独的对象一样(参见第一章、入门一节基本类型)。

数组视图

一个最简单的视图示例由一个数组的切片给出:

M = array([[1.,2.],[3.,4.]])
v = M[0,:] # first row of M

前面的切片是M的视图。它与M共享相同的数据。修改v也会修改M:

v[-1] = 0.
v # array([[1.,0.]])
M # array([[1.,0.],[3.,4.]]) # M is modified as well

可以使用数组属性base访问拥有数据的对象:

v.base # array([[1.,0.],[3.,4.]])
v.base is M # True

如果数组拥有其数据，则属性基为 none:

M.base # None

切片作为视图

哪些切片将返回视图，哪些将返回副本，都有精确的规则。只有基本切片(主要是带有:的索引表达式)返回视图，而任何高级选择(例如用布尔值切片)都将返回数据的副本。例如，可以通过列表(或数组)索引来创建新的矩阵:

a = arange(4) # array([0.,1.,2.,3.])
b = a[[2,3]] # the index is a list [2,3]
b # array([2.,3.])
b.base is None # True, the data was copied
c = a[1:3]
c.base is None # False, this is just a view

在前面的例子中，数组b不是视图，而用更简单的切片获得的数组c是视图。

有一个特别简单的数组切片，它返回整个数组的视图:

N = M[:] # this is a view of the whole array M

作为视图进行置换和重塑

其他一些重要的操作返回视图。例如，转置返回一个视图:

M = random.random_sample((3,3))
N = M.T
N.base is M # True

这同样适用于所有整形操作:

v = arange(10)
C = v.reshape(-1,1) # column matrix
C.base is v # True

数组副本

有时需要明确请求复制数据。这可以简单地通过名为array的 NumPy 函数来实现:

M = array([[1.,2.],[3.,4.]])
N = array(M.T) # copy of M.T

我们可能会验证数据确实已被复制:

N.base is None # True

比较数组

比较两个数组并不像看起来那么简单。考虑下面的代码，它旨在检查两个矩阵是否彼此接近:

A = array([0.,0.])
B = array([0.,0.])
if abs(B-A) < 1e-10: # an exception is raised here
    print("The two arrays are close enough")

当执行if语句时，该代码引发异常:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

在本节中，我们将解释为什么会这样，以及如何补救这种情况。

布尔数组

布尔数组对于高级数组索引很有用(参见用布尔数组索引一节)。布尔数组只是一个条目类型为bool的数组:

A = array([True,False]) # Boolean array
A.dtype # dtype('bool')

任何作用于数组的比较运算符都将创建一个布尔数组，而不是简单的布尔:

M = array([[2, 3],
           [1, 4]])
M > 2 # array([[False, True],
             # [False, True]])
M == 0 # array([[False, False],
             # [False, False]])
N = array([[2, 3],
           [0, 0]])
M == N # array([[True, True],
              # [False, False]])
...

请注意，因为数组比较创建布尔数组，所以不能在条件语句中直接使用数组比较，例如，if语句。解决方法是使用方法all和any:

A = array([[1,2],[3,4]])
B = array([[1,2],[3,3]])
A == B # creates array([[True, True], [True, False]]) 
(A == B).all() # False
(A != B).any() # True
if (abs(B-A) < 1e-10).all():
    print("The two arrays are close enough")

检查是否相等

检查两个浮点数组的相等性并不是直截了当的，因为两个浮点可能非常接近而不相等。在 NumPy 中，可以检查与allclose是否相等。该函数检查给定精度下两个数组的相等性:

data = random.rand(2)*1e-3
small_error = random.rand(2)*1e-16
data == data + small_error # False
allclose(data, data + small_error, rtol=1.e-5, atol=1.e-8)   # True

公差以相对公差界限rtol和绝对误差界限atol给出。命令allclose是:(abs(A-B) < atol+rtol*abs(B)).all()的缩写。

注意allclose也可以应用于标量:

data = 1e-3
error = 1e-16
data == data + error # False
allclose(data, data + error, rtol=1.e-5, atol=1.e-8)  #True

数组上的布尔运算

布尔数组中不能使用and、or和not。事实上，这些操作符强制从数组转换为布尔值，这是不允许的。相反，我们可以使用下表中给出的运算符(表 5.1 )对布尔数组进行组件逻辑运算:

| **逻辑运算符** | **布尔数组的替换** | | `A and B` | `A & B` | | `A or B` | `A | B` | | `not A` | `~ A` |

表 5.1 逻辑运算符和、或和不适用于数组。

A = array([True, True, False, False])
B = array([True, False, True, False])
A and B # error!
A & B # array([True, False, False, False])
A | B # array([True, True, True, False])
~A # array([False, False, True, True])

以下是逻辑运算符在布尔数组中的用法示例:

假设我们有一系列的数据被一些测量误差破坏了。进一步假设我们运行一个回归，它给出了每个值的偏差。我们希望获得所有异常值和所有偏差小于给定阈值的值:

data = linspace(1,100,100) # data
deviation = random.normal(size=100) # the deviations 
           #don't forget the parentheses in next statement!
exceptional = data[(deviation<-0.5)|(deviation>0.5)] 
exceptional = data[abs(deviation)>0.5] # same result 
small = data[(abs(deviation)<0.1)&(data<5.)] # small deviation and data

数组索引

我们已经看到，可以通过切片和整数的组合来索引数组，这是基本的切片技术。然而，还有很多可能性，允许多种方式来访问和修改数组元素。

用布尔数组索引

根据数组的值，只访问和修改数组的一部分通常很有用。例如，人们可能想要访问数组的所有正元素。事实证明，使用布尔数组是可能的，布尔数组就像掩码一样，只选择数组的某些元素。这种索引的结果是总是一个向量。例如，考虑以下示例:

B = array([[True, False],
           [False, True]])
M = array([[2, 3],
           [1, 4]])
M[B] # array([2,4]), a vector

其实M[B]通话相当于M.flatten()[B]。然后，可以用另一个向量替换结果向量。例如，可以用零替换所有元素(更多信息请参考广播部分):

M[B] = 0
M # [[0, 3], [1, 0]]

或者可以用其他值替换所有选定的值:

M[B] = 10, 20
M # [[10, 3], [1, 20]]

通过结合创建布尔数组(M > 2)、智能索引(用布尔数组索引)和广播，可以使用以下优雅的语法:

M[M>2] = 0    # all the elements > 2 are replaced by 0

这里的表达广播指的是标量 0 到适当形状的向量的默认转换。

使用位置

命令where给出了一个有用的构造，它可以将布尔数组作为一个条件，并返回满足该条件的数组元素的索引，或者根据布尔数组中的值返回不同的值。

基本结构是:

where(condition, a, b)

这将在条件为True时返回a的值，在条件为False时返回b的值。

例如，考虑一个重量函数:

下面的代码实现了一个 Heaviside 函数:

def H(x):
    return where(x < 0, 0, 1)
x = linspace(-1,1,11)  # [-1\. -0.8 -0.6 -0.4 -0.2 0\. 0.2 0.4 0.6 0.8 1\. ]
print(H(x))            # [0 0 0 0 0 1 1 1 1 1 1]

第二个和第三个参数可以是与条件大小相同的数组(布尔数组)，也可以是标量。我们再举两个例子来演示如何根据条件操作数组或标量中的元素:

x = linspace(-4,4,5)
# [-4\. -2.  0.  2.  4.]

print(where(x > 0, sqrt(x), 0))
# [ 0.+0.j 0.+0.j 0.+0.j 1.41421356+0.j  2.+0.j ]
print(where(x > 0, 1, -1)) # [-1 -1 -1  1  1]

如果省略第二个和第三个参数，则返回包含满足条件的元素索引的元组。

例如，考虑以下代码中只有一个参数的where的使用:

a = arange(9)
b = a.reshape((3,3))

print(where(a > 5))   # (array([6, 7, 8]),)

print(where(b > 5))   # (array([2, 2, 2]), array([0, 1, 2]))

性能和矢量化

说到 Python 代码的性能，通常归结为解释代码和编译代码之间的区别。Python 是一种解释的编程语言，基本的 Python 代码是直接执行的，不需要对机器代码进行任何中间编译。使用编译语言，代码需要在执行前被翻译成机器指令。

解释语言的好处很多，但是解释代码在速度上无法与编译代码竞争。为了让你的代码更快，你可以用像 FORTRAN、C 或 C++这样的编译语言编写一些部分。这就是 NumPy 和 SciPy 所做的。

因此，最好尽可能在解释版本上使用 NumPy 和 SciPy 中的函数。矩阵乘法、矩阵向量乘法、矩阵分解、标量积等 NumPy 数组运算的速度比任何纯 Python 等价类都要快得多。考虑标量积的简单情况。标量积比编译的 NumPy 函数dot(a,b)慢得多(对于大约有 100 个元素的数组，慢 100 多倍):

def my_prod(a,b):
    val = 0
    for aa,bb in zip(a,b):
        val += aa*bb
    return val

测量函数的速度是科学计算的一个重要方面。测量执行时间详见第十三章测试中测量执行时间一节。

矢量化

为了提高性能，必须经常向量化代码。用 NumPy 切片、操作和函数替换for循环和其他较慢的代码部分可以带来显著的改进。例如，通过迭代元素将标量简单地添加到向量中是非常慢的:

for i in range(len(v)):
    w[i] = v[i] + 5

其中使用 NumPy 的加法要快得多:

w = v + 5

使用 NumPy 切片也可以在迭代for循环时显著提高速度。为了演示这一点，让我们考虑在二维数组中形成邻居的平均值:

def my_avg(A):
    m,n = A.shape
    B = A.copy()
    for i in range(1,m-1):
        for j in range(1,n-1):
            B[i,j] = (A[i-1,j] + A[i+1,j] + A[i,j-1] + A[i,j+1])/4
    return B

def slicing_avg(A):
    A[1:-1,1:-1] = (A[:-2,1:-1] + A[2:,1:-1] +
    A[1:-1,:-2] + A[1:-1,2:])/4
    return A

这些函数都给每个元素分配了四个相邻元素的平均值。使用切片的第二个版本要快得多。

除了用 NumPy 函数代替for循环和其他较慢的构造，还有一个有用的函数叫做vectorize，参考第四章、线性代数-数组中函数作用于数组一节。这将采用一个函数，并创建一个矢量化版本，尽可能使用函数将该函数应用于数组的所有元素。

考虑以下向量化函数的示例:

def my_func(x):
    y = x**3 - 2*x + 5
    if y>0.5:
        return y-0.5
    else:
        return 0

通过迭代数组来应用这一点非常慢:

for i in range(len(v)):
    v[i] = my_func(v[i])

相反，使用vectorize创建一个新函数，如下所示:

my_vecfunc = vectorize(my_func)

然后，该函数可以直接应用于数组:

v = my_vecfunc(v)

矢量化选项要快得多(长度为 100 的数组快 10 倍左右)。

广播

NumPy 中的广播表示猜测两个数组之间共同的兼容形状的能力。例如，当添加一个向量(一维数组)和一个标量(零维数组)时，标量被扩展为向量，以便允许添加。一般的机制叫做广播。我们将首先从数学的角度回顾这个机制，然后继续给出 NumPy 中广播的精确规则。

数学观

广播经常在数学中进行，主要是隐含的。例如 f(x) + C 或 f(x) + g(y) 等表达式。我们将在本节中对该技术进行明确的描述。

我们牢记函数和 NumPy 数组之间非常密切的关系，如第 4 章线性代数-数组的数学预备知识一节所述。

常量函数

广播最常见的例子之一是增加一个函数和一个常数；如果 C 是一个标量，人们通常会写道:

这是对符号的滥用，因为人们不应该能够添加函数和常数。然而，常量被隐式地传播给函数。常数 C 的广播版本是功能定义的:

现在将两个函数加在一起是有意义的:

我们并不是为了学究而学究，而是因为类似的情况可能会出现在数组中，如下面的代码所示:

vector = arange(4) # array([0.,1.,2.,3.])
vector + 1\.        # array([1.,2.,3.,4.])

在这个例子中，发生的一切就好像标量1.已经被转换成与vector长度相同的数组，即array([1.,1.,1.,1.])，然后被添加到vector中。

这个例子非常简单，所以我们继续展示不太明显的情况。

几个变量的函数

一个更复杂的广播例子出现在构建几个变量的函数时。例如，假设给我们一个变量的两个函数， f 和 g ，我们想要根据公式构造一个新的函数 F :

这显然是一个有效的数学定义。我们想把这个定义表达为两个变量中两个函数的和，这两个变量定义为

现在我们可以简单地写道:

这种情况类似于添加列矩阵和行矩阵时出现的情况:

C = arange(2).reshape(-1,1) # column
R = arange(2).reshape(1,-1) # row
C + R                       # valid addition: array([[0.,1.],[1.,2.]])

这在对两个变量的函数进行采样时特别有用，如部分典型示例所示。

一般机制

我们已经看到了如何添加一个函数和一个标量，以及如何从一个变量的两个函数构建一个两个变量的函数。现在让我们把重点放在使这成为可能的一般机制上。一般的机制包括两个步骤:重塑和延伸。

首先，函数 g 被重塑为一个函数，该函数接受两个参数。其中一个论点是一个伪论点，我们认为它是零，作为一个惯例:

数学上，的定义域现在是然后函数 f 以类似于以下的方式被重塑:

现在和都取两个参数，虽然其中一个总是零。我们继续下一步，扩展。将常数转换为常数函数的步骤相同(参考常数函数示例)。

功能扩展为:

功能扩展为:

现在两个变量 F 的函数，由 F(x，y) = f(x) + g(y) 草率地定义，可以不用参考它的参数来定义:

例如，让我们描述一下前面的常数机制。常数是标量，也就是零参数的函数。因此，整形步骤是定义一个(空的)变量的函数:

现在，扩展步骤简单地通过以下方式进行:

惯例

最后一个要素是关于如何向函数添加额外参数的约定，即如何自动执行整形。按照惯例，一个函数通过在左边加零来自动重塑。

例如，如果两个参数的函数 g 必须被重新整形为三个参数，那么新的函数将被定义为:

广播数组

我们现在重复观察数组仅仅是几个变量的函数(参考第四章、线性代数-数组中的数学预备知识一节)。因此，数组广播遵循与上述数学函数完全相同的过程。广播在 NumPy 中自动完成。

在下图(图 5.1 )中，我们展示了将形状矩阵(4，3)添加到大小矩阵(1，3)时会发生什么。第二矩阵具有形状(4，3):

图 5.1:矩阵和向量之间的广播。

广播问题

当 NumPy 被赋予两个不同形状的数组，并被要求执行要求这两个形状相同的操作时，这两个数组被广播到一个公共形状。

假设两个数组具有形状s₁T3】和s₂T7】。该广播分两步进行:**

1. 如果形状s₁T3 比形状s₂T7】短，则在形状 s ₁ 的左侧添加 1。这是一次重塑。**
2. 当形状具有相同的长度时，数组被扩展以匹配形状 s ₂ (如果可能的话)。

假设我们想要将形状(3)的向量添加到形状(4，3)的矩阵中。矢量需要广播。第一次手术是整形；向量的形状从(3)转换为(1，3)。第二个操作是扩展；形状从(1，3)转换为(4，3)。

例如，假设大小为 n 的向量将被广播到形状( m，n ):

1. v 自动重塑为(1， n )。
2. v 延伸至( m 、 n )。

为了证明这一点，我们考虑由下式定义的矩阵:

M = array([[11, 12, 13, 14],
           [21, 22, 23, 24],
           [31, 32, 33, 34]])

和向量由下式给出:

v = array([100, 200, 300, 400])

现在我们可以直接添加M和v:

M + v # works directly

结果是这个矩阵:

形状不匹配

不可能自动将长度为n的向量v传播到形状(n,m)。下图说明了这一点:

播放会失败，因为形状(n,)可能不会自动播放到形状(m, n)。解决办法是手动将v重塑为(n,1)的形状。广播现在将照常工作(仅扩展):

M + v.reshape(-1,1)

下面是另一个示例，通过以下方式定义矩阵:

M = array([[11, 12, 13, 14],
           [21, 22, 23, 24],
           [31, 32, 33, 34]])

和一个矢量:

v = array([100, 200, 300])

现在自动广播将失败，因为自动整形不起作用:

M + v # shape mismatch error

因此，解决方案是手动处理整形。在那种情况下，我们想要的是在右边加 1，也就是把向量转换成列矩阵。然后，广播直接工作:

M + v.reshape(-1,1)

形状参数-1，参见第四章、线性代数-数组的访问和更改形状*部分。*结果就是这个矩阵:

典型例子

让我们研究一些广播可能派上用场的典型例子。

重新缩放行

假设M是一个 n × m 矩阵，我们要将每行乘以一个系数。系数存储在一个向量coeff中，该向量具有 n 个分量。在这种情况下，自动重塑将不起作用，我们必须执行:

rescaled = M*coeff.reshape(-1,1)

重新缩放列

这里的设置是相同的，但是我们想要用存储在长度为 m 的向量coeff中的系数来重新缩放每一列。在这种情况下，自动整形将起作用:

rescaled = M*coeff

显然，我们也可以手动进行整形，并通过以下方式获得相同的结果:

rescaled = M*coeff.reshape(1,-1)

两个变量的函数

假设 u 和 v 是向量，我们想用元素W_ij= u_I+v_j组成矩阵。这对应于函数 F(x，y) = x + y 。矩阵 W 仅由下式定义:

W=u.reshape(-1,1) + v

如果向量 u 和 v 分别为[0，1]和[0，1，2]，则结果为:

更一般地说，假设我们要对函数 w ( x，y):=cos(x)+sin(2y)。假设向量 x 和 y 被定义，采样值的矩阵 w 由下式得到:

w = cos(x).reshape(-1,1) + sin(2*y)

请注意，这与ogrid结合使用非常频繁。从ogrid获得的向量已经被方便地成形用于广播。这允许对函数 cos(x)+sin(2y)进行以下优雅的采样:

x,y = ogrid[0:1:3j,0:1:3j] 
# x,y are vectors with the contents of linspace(0,1,3)
w = cos(x) + sin(2*y)

ogrid的语法需要一些解释。第一，ogrid是没有功能的。这是一个具有__getitem__方法的类的实例(参见第 8 章、类中的属性一节)。这就是为什么它用括号代替圆括号的原因。

这两个命令是等效的:

x,y = ogrid[0:1:3j, 0:1:3j]
x,y = ogrid.__getitem__((slice(0, 1, 3j),slice(0, 1, 3j)))

前面例子中的步幅参数是一个复数。这表示是步数而不是步长。乍一看，步幅参数的规则可能会令人困惑:

• 如果步幅是实数，那么它定义了开始和停止之间的步长，并且停止不包括在列表中。
• 如果步幅是复数s，那么s.imag的整数部分定义了开始和停止之间的步数，停止包含在列表中。

ogrid输出的另一个例子是具有两个数组的元组，可以用于广播:

x,y = ogrid[0:1:3j, 0:1:3j]

给出:

array([[ 0\. ],
       [ 0.5],
       [ 1\. ]])
array([[ 0\. ,  0.5,  1\. ]])

这相当于:

x,y = ogrid[0:1.5:.5, 0:1.5:.5]

稀疏矩阵

具有少量非零条目的矩阵称为稀疏矩阵。例如，当在数值求解偏微分方程的上下文中描述离散微分算子时，在科学计算中会出现稀疏矩阵。

稀疏矩阵通常有很大的维数，有时大到整个矩阵(零条目)甚至无法容纳在可用的内存中。这是稀疏矩阵特殊类型的一个动机。另一个动机是可以避免零矩阵条目的更好的操作性能。

线性代数中一般的非结构化稀疏矩阵的算法数量非常有限。其中大多数本质上是迭代的，并且基于稀疏矩阵的矩阵向量乘法的有效实现。

稀疏矩阵的例子是对角矩阵或带状矩阵。这些矩阵的简单模式允许简单的存储策略；主对角线以及次对角线和超对角线存储在 1D 数组中。从稀疏表示到经典数组类型的转换可以通过命令diag完成，反之亦然。

一般来说，没有这么简单的结构，稀疏矩阵的描述需要特殊的技术和标准。在这里，我们为稀疏矩阵提供了一个面向行和一个面向列的类型，这两个类型都可以通过模块scipy.sparse获得。

图 5.2:弹性板有限元模型的刚度矩阵。像素表示 1250 × 1250 矩阵中的非零条目

稀疏矩阵格式

scipy.sparse模块从稀疏矩阵中提供许多不同的存储格式。这里我们只描述最重要的几个:中国南车、中国南车和 LIL。LIL 格式应该用于生成和修改稀疏矩阵；CSR 和 CSC 是矩阵-矩阵和矩阵-向量运算的有效格式。

压缩稀疏行

压缩稀疏行格式(CSR)使用三个数组:data、indptr和indices:

• 1D 数组data按顺序存储所有非零值。它的元素和非零元素一样多，通常用变量nnz表示。
• 1D 数组indptr包含整数，使得indptr[i]是data中元素的索引，该元素是行 i 的第一个非零元素。如果整行 i 为零，则indptr[i]==indptr[i+1]。如果原矩阵有 m 行，那么len(indptr)==m+1。
• 1D 数组indices包含列索引信息的方式是indices[indptr[i]:indptr[i+1]]是一个整数数组，行 i 中非零元素的列索引。显然，len(indices)==len(data)==nnz。

让我们看一个例子:矩阵的企业社会责任格式:

由三个数组给出:

data = (1\. 2\. 3\. 4.)
indptr = (0 2 2 3 5)
indices = (0 2 0 0 3)

模块scipy.sparse提供了一个类型csr_matrix，带有一个构造函数，可以通过以下方式使用:

• 以 2D 数组作为参数
• 在scipy.sparse中使用其他稀疏格式之一的矩阵
• 使用形状参数(m,n)，生成 CSR 格式的零矩阵
• 通过用于data的 1D 数组和具有形状(2,len(data))的整数数组ij，使得ij[0,k]是矩阵的行索引，ij[1,k]是矩阵的data[k]的列索引
• 三个参数data、indptr和indices可以直接给构造函数

前两个选项用于转换目的，而后两个选项直接定义稀疏矩阵。

考虑 python 中的上述示例，如下所示:

import scipy.sparse as sp
A = array([[1,0,2,0],[0,0,0,0],[3.,0.,0.,0.],[1.,0.,0.,4.]])
AS = sp.csr_matrix(A)

其中，提供了以下属性:

AS.data      # returns array([ 1., 2., 3., 1., 4.]) 
AS.indptr    # returns array([0, 2, 2, 3, 5])
AS.indices   # returns array([0, 2, 0, 0, 3])
AS.nnz       # returns 5

压缩稀疏列

CSR 格式有一个面向列的孪生格式——压缩稀疏列(CSC)格式。与 CSR 格式相比，它唯一的不同是indptr和indices数组的定义，这两个数组现在是与列相关的。CSC 格式的类型是csc_matrix，它的使用对应于csr_matrix，在本节前面已经解释过了。

继续 CSC 格式的相同示例:

import scipy.sparse as sp
A = array([[1,0,2,0],[0,0,0,0],[3.,0.,0.,0.],[1.,0.,0.,4.]])
AS = sp.csc_matrix(A)
AS.data         # returns array([ 1., 3., 1., 2., 4.]) 
AS.indptr       # returns array([0, 3, 3, 4, 5])
AS.indices      # returns array([0, 2, 3, 0, 3])
AS.nnz          # returns 5

基于行的链表格式

链表稀疏格式将非零矩阵条目按行存储在列表data中，使得data[k]是行 k 中非零条目的列表。如果该行中的所有条目都为 0，则它包含一个空列表。

第二个列表rows在位置 k 包含行 k 中非零元素的列索引列表。下面是基于行的链表格式(LIL) 格式的一个例子:

import scipy.sparse as sp
A = array([[1,0,2,0],[0,0,0,0], [3.,0.,0.,0.], [1.,0.,0.,4.]]) 
AS = sp.lil_matrix(A)
AS.data     # returns array([[1.0, 2.0], [], [3.0], [1.0, 4.0]], dtype=object)
AS.rows     # returns array([[0, 2], [], [0], [0, 3]], dtype=object)
AS.nnz      # returns 5

以 LIL 格式更改和切片矩阵

LIL 格式最适合切片，即提取 LIL 格式的子矩阵，并通过插入非零元素来改变稀疏模式。下一个示例演示了切片:

BS = AS[1:3,0:2]
BS.data     # returns array([[], [3.0]], dtype=object)
BS.rows     # returns array([[], [0]], dtype=object)

插入新的非零元素会自动更新属性:

AS[0,1] = 17 
AS.data # returns array([[1.0, 17.0, 2.0], [], [3.0], [1.0, 4.0]])
AS.rows              # returns array([[0, 1, 2], [], [0], [0, 3]])
AS.nnz               # returns 6

其他稀疏矩阵格式不鼓励这些操作，因为它们效率极低。

生成稀疏矩阵

NumPy 命令eye、identity、diag和rand有它们稀疏的对应物。他们提出了一个额外的论点；它指定结果矩阵的稀疏矩阵格式。

以下命令生成单位矩阵，但采用不同的稀疏矩阵格式:

import scipy.sparse as sp
sp.eye(20,20,format = 'lil') 
sp.spdiags(ones((20,)),0,20,20, format = 'csr') 
sp.identity(20,format ='csc')

sp.rand命令接受描述生成的随机矩阵密度的附加参数。密集矩阵的密度为 1，而零矩阵的密度为 0:

import scipy.sparse as sp 
AS=sp.rand(20,200,density=0.1,format=’csr’)
AS.nnz # returns 400

与 NumPy 命令zeroes没有直接对应关系。完全用零填充的矩阵是通过实例化相应的类型来生成的，其中形状参数作为构造函数参数:

import scipy.sparse as sp
Z=sp.csr_matrix((20,200))
Z.nnz    # returns 0

稀疏矩阵方法

有几种方法可以将一种稀疏类型转换为另一种类型或转换为数组:

AS.toarray # converts sparse formats to a numpy array 
AS.tocsr
AS.tocsc
AS.tolil

稀疏矩阵的类型可以通过方法issparse、isspmatrix_lil、isspmatrix_csr和isspmatrix_csc来检查。

稀疏矩阵上的元素操作+、*、/和**被定义为用于 NumPy 数组。不管操作数的稀疏矩阵格式如何，结果总是一个csr_matrix。将元素操作函数应用于稀疏矩阵需要首先将它们转换为 CSR 或 CSC 格式，并将函数应用于它们的data属性，如下一个示例所示。

稀疏矩阵的元素正弦可以通过对其data属性的操作来定义:

import scipy.sparse as sp
def sparse_sin(A):
    if not (sp.isspmatrix_csr(A) or sp.isspmatrix_csc(A)):
        A = A.tocsr()
A.data = sin(A.data)
return A

对于矩阵-矩阵或矩阵-向量乘法，有一种稀疏矩阵方法，dot。它返回一个csr_matrix或一个array1D 号码:

import scipy.sparse as sp
A = array([[1,0,2,0],[0,0,0,0],[3.,0.,0.,0.],[1.,0.,0.,4.]])
AS = sp.csr_matrix(A)
b = array([1,2,3,4])
c = AS.dot(b)      # returns array([ 7., 0., 3., 17.]) 
C = AS.dot(AS)     # returns  csr_matrix
d = dot(AS,b)      # does not return the expected result! 

型式

避免在稀疏矩阵上使用 NumPy 的命令dot，因为这可能会导致意想不到的结果。请使用来自scipy.sparse的命令dot。

其他线性代数运算，如系统求解、最小二乘、特征值和奇异值由scipy.sparse.linalg模块提供。

总结

观点的概念是你应该从本章中学到的重要话题之一。错过这个主题会让你在调试代码时很困难。布尔数组出现在本书的不同地方。在处理数组时，它们是避免冗长的if构造和循环的方便而紧凑的工具。在几乎所有大型计算项目中，稀疏矩阵都成为一个问题。您看到了这些是如何处理的，以及哪些相关方法可用。

六、绘图

Python 中的绘图可以通过 matplotlib 模块的pyplot部分来完成。使用 matplotlib，您可以创建高质量的图形和图形，还可以绘制和可视化您的结果。Matplotlib 是开源的免费软件，【21】。matplotlib 网站还包含优秀的文档和示例， [35] 。在本节中，我们将向您展示如何使用最常见的功能。接下来几节中的示例假设您已经将模块导入为:

from matplotlib.pyplot import *

如果想使用 IPython 中的绘图命令，建议在启动 IPython shell 后直接运行魔法命令%matplotlib。这为交互式绘图做好了准备。

基本标绘

标准绘图功能为plot。调用plot(x,y)创建一个图形窗口，其图形为 y 作为 x 的函数。输入参数是长度相等的数组(或列表)。也可以使用plot(y)，在这种情况下 y 中的值将相对于它们的指数绘制，即plot(y)是plot(range(len(y)),y)的缩写。

下面是一个示例，展示了如何使用 200 个采样点为xϵ【-2π，2π】绘制 sin( x )并在每第四个点设置标记:

# plot sin(x) for some interval
x = linspace(-2*pi,2*pi,200)
plot(x,sin(x))

# plot marker for every 4th point
samples = x[::4]
plot(samples,sin(samples),'r*')

# add title and grid lines
title('Function sin(x) and some points plotted')
grid()

结果如下图所示(图 6.1 ):

图 6.1:显示了带有网格线的函数 sin(x)的曲线图。

如您所见，标准图是一条纯蓝色曲线。每个轴会自动缩放以适应这些值，但也可以手动设置。颜色和绘图选项可以在前两个输入参数之后给出。这里，r*表示红色的星形标记。下一节将更详细地介绍格式化。title命令在绘图区域上方放置一个标题文本字符串。

多次调用plot会将图叠加在同一个窗口中。要获得新的干净图形窗口，请使用figure()。figure命令可能包含一个整数，例如figure(2)，可用于在图形窗口之间切换。如果没有具有该编号的图形窗口，则会创建一个新的图形窗口，否则，该窗口将被激活以进行打印，并且所有后续的打印命令都将应用于该窗口。

可以使用legend功能解释多个绘图，并为每个绘图调用添加标签。以下示例使用命令polyfit和polyval将多项式拟合到一组点，并用图例绘制结果:

# —Polyfit example—
x = range(5)
y = [1,2,1,3,5]
p2 = polyfit(x,y,2)
p4 = polyfit(x,y,4)

# plot the polynomials and points
xx = linspace(-1,5,200) 
plot(xx, polyval(p2, xx), label='fitting polynomial of degree 2')
plot(xx, polyval(p4, xx),
                label='interpolating polynomial of degree 4') 
plot(x,y,'*')

# set the axis and legend
axis([-1,5,0,6])
legend(loc='upper left', fontsize='small')

这里还可以看到如何使用axis([xmin,xmax,ymin,ymax])手动设置轴的范围。legend命令接受关于位置和格式的可选参数；这种情况下图例放在左上角，用小字号排版，如下图所示(图 6.2 )。

图 6.2:拟合到相同点的两个多项式。

作为基本绘图的最后一个示例，我们演示了如何在二维中绘制散点图和对数图。

2D 点散点图示例:

# create random 2D points
import numpy
x1 = 2*numpy.random.standard_normal((2,100))
x2 = 0.8*numpy.random.standard_normal((2,100)) + array([[6],[2]])
plot(x1[0],x1[1],'*')
plot(x2[0],x2[1],'r*')
title('2D scatter plot')

图 6.3(a):散点图的一个例子

以下代码是使用loglog的对数图的示例:

# log both x and y axis 
x = linspace(0,10,200) 
loglog(x,2*x**2, label = 'quadratic polynomial',
                            linestyle = '-', linewidth = 3)
loglog(x,4*x**4, label = '4th degree polynomial',
                            linestyle = '-.', linewidth = 3)
loglog(x,5*exp(x), label = 'exponential function', linewidth = 3)
title('Logarithmic plots')
legend(loc = 'best')

图 6.3(b):带有对数 x 轴和 y 轴的曲线图示例

上图所示的例子(图 6.3(a) 和图 6.3(b) )使用了plot和loglog的一些允许特殊格式化的参数。在下一节中，我们将更详细地解释参数。

格式化

图形和情节的外观可以根据您的需要进行设计和定制。一些重要的变量是linewidth，控制地块线的粗细；设置轴标签的xlabel、ylabel、color为地块颜色，transparent为透明度。本节将告诉您如何使用其中的一些。下面是一个包含更多关键字的示例:

k = 0.2
x = [sin(2*n*k) for n in range(20)]
plot(x, color='green', linestyle='dashed', marker='o', 
                       markerfacecolor='blue', markersize=12, linewidth=6)

如果您只需要基本的样式更改，例如设置颜色和线条样式，可以使用一些简短的命令。下表(表 6.1 )显示了这些格式化命令的一些示例。您可以使用短字符串语法plot(...,'ro-')，或者更明确的语法plot(..., marker='o', color='r', linestyle='-')。

表 6.1:一些常见的绘图格式参数

要用我们写的'o'标记将颜色设置为绿色:

plot(x,'go')

要绘制直方图而不是常规图，使用hist命令:

# random vector with normal distribution
sigma, mu = 2, 10
x = sigma*numpy.random.standard_normal(10000)+mu 
hist(x,50,normed=1)
z = linspace(0,20,200)
plot(z, (1/sqrt(2*pi*sigma**2))*exp(-(z-mu)**2/(2*sigma**2)),'g')
# title with LaTeX formatting 
title('Histogram with '.format(mu,sigma))

图 6.4 具有 50 个箱的正态分布和指示真实分布的绿色曲线

结果图与上图相似(图 6.4 )。标题和任何其他文本都可以使用 LaTeX 进行格式化，以显示数学公式。LaTeX 格式包含在一对$符号中。另外，请注意使用format方法进行的字符串格式化，请参考第 2 章、变量和基本类型中的字符串一节。

有时字符串格式的括号会干扰 LaTeX 括号环境。如果出现这种情况，用双支架替换 LaTeX 支架，例如x_{1}应替换为x_{{1}}。文本可能包含与字符串转义序列重叠的序列，例如，\tau将被解释为制表符\t。一个简单的解决方法是在字符串前添加r，例如r'\tau'；这使它成为一个原始字符串。

使用subplot命令可以在一个图形窗口中放置几个图。考虑下面的例子，它迭代地平均正弦曲线上的噪声。

def avg(x):
    """ simple running average """
    return (roll(x,1) + x + roll(x,-1)) / 3
# sine function with noise
x = linspace(-2*pi, 2*pi,200)
y = sin(x) + 0.4*rand(200)

# make successive subplots
for iteration in range(3):
    subplot(3, 1, iteration + 1)
    plot(x,y, label = '{:d} average{}'.format(iteration, 's' if iteration > 1 else ''))
    yticks([])
    legend(loc = 'lower left', frameon = False)
    y = avg(y) #apply running average 
subplots_adjust(hspace = 0.7)

图 6.5:在同一图形窗口中绘制多次的示例。

函数avg使用roll调用来移动数组的所有值。subplot采用三个参数:垂直图的数量、水平图的数量和指示绘制位置的索引(按行计数)。请注意，我们使用subplots_adjust命令添加额外的空间来调整两个支线剧情之间的距离。

一个有用的命令是savefig，它允许你将一个图形保存为一个图像(这也可以在图形窗口中完成)。该命令支持许多图像和文件格式，它们由文件扩展名指定为:

savefig('test.pdf')  # save to pdf

或者

savefig('test.svg')  # save to svg (editable format)

您可以将图像放在非白色背景下，例如网页。为此，可以设置transparent参数使图形背景透明:

savefig('test.pdf', transparent=True)

如果您打算将图形嵌入到 LaTeX 文档中，建议您通过将图形的边界框紧密设置在绘图周围来减少周围的空白，如下所示:

savefig('test.pdf', bbox_inches='tight')

网格和等高线

常见的任务是矩形上标量函数的图形表示:

为此，首先我们必须在矩形[ a 、 b ] x [ c 、 d 上生成一个网格。这是使用meshgrid命令完成的:

n = ... # number of discretization points along the x-axis
m = ... # number of discretization points along the x-axis 
X,Y = meshgrid(linspace(a,b,n), linspace(c,d,m))

X和Y为(n,m)形状的数组，使得包含网格点的坐标，如下图*(图 6.6)* :

图 6.6:由网格离散的矩形

由meshgrid离散的矩形将用于可视化迭代的行为。但是首先我们将使用它来绘制函数的水平曲线。这是通过命令contour完成的。

例如，我们选择罗森布鲁克的香蕉函数:

它用于挑战优化方法。函数值向香蕉形谷下降，该谷本身向函数的全局最小值(1，1)缓慢下降。

首先我们使用contour显示水平曲线。

rosenbrockfunction = lambda x,y: (1-x)**2+100*(y-x**2)**2 
X,Y = meshgrid(linspace(-.5,2.,100), linspace(-1.5,4.,100))
Z = rosenbrockfunction(X,Y) 
contour(X,Y,Z,logspace(-0.5,3.5,20,base=10),cmap='gray') 
title('Rosenbrock Function: ')
xlabel('x')
ylabel('y')

这将在第四个参数给出的级别绘制级别曲线，并使用颜色图gray。此外，我们使用从 10 ^0.5 到 10 ³ 的对数间隔步长，使用函数logscale来定义级别，如下图所示。

图 6.7:罗森布鲁克函数的等高线图

在前面的例子中，由关键字lambda指示的匿名函数用于保持代码紧凑。匿名函数在第 7 章、函数、匿名函数的匿名函数 lambda 关键字一节中进行了解释。如果等级没有作为参数提供给contour，函数会自己选择合适的等级。

contourf功能执行与contour相同的功能，但根据不同的级别用颜色填充绘图。等高线图是可视化数值方法行为的理想选择。我们在这里通过展示优化方法的迭代来说明这一点。

我们继续前面的例子，描述了由鲍威尔方法【27】生成的罗森布罗克函数的最小值的步骤，我们将应用它来寻找罗森布罗克函数的最小值:

import scipy.optimize as so
rosenbrockfunction = lambda x,y: (1-x)**2+100*(y-x**2)**2
X,Y=meshgrid(linspace(-.5,2.,100),linspace(-1.5,4.,100))
Z=rosenbrockfunction(X,Y)
cs=contour(X,Y,Z,logspace(0,3.5,7,base=10),cmap='gray')
rosen=lambda x: rosenbrockfunction(x[0],x[1])
solution, iterates = so.fmin_powell(rosen,x0=array([0,-0.7]),retall=True)
x,y=zip(*iterates)
plot(x,y,'ko') # plot black bullets
plot(x,y,'k:',linewidth=1) # plot black dotted lines
title("Steps of Powell's method to compute a  minimum")
clabel(cs)

迭代法fmin_powell应用鲍威尔法求最小值。它由给定的开始值 x ₀ 启动，并在给出选项retall=True时报告所有迭代。经过十六次迭代，找到了解 x= 0 ，y= 0。迭代在下面的等高线图中被描述为项目符号(图 6.8 )。

图 6.8:罗森布鲁克函数的等高线图，带有优化方法的搜索路径

contour还创建了一个轮廓集对象，我们将其分配给变量cs。然后clabel用它来标注相应函数值的级别，如上图所示(图 6.8 )。

图像和轮廓

让我们看一些将数组可视化为图像的例子。下面的函数将为曼德勃罗分形**创建一个颜色值矩阵。**这里我们考虑一个不动点迭代，它依赖于一个复杂的参数 c :

根据该参数的选择，它可能会也可能不会创建复数值的有界序列 z _n 。

对于 c 的每个值，我们检查 z _n 是否超过规定的界限。如果在maxit次迭代中，它保持在界限之下，我们假设序列是有界限的。

请注意，在下面这段代码中，meshgrid是如何用来生成复杂参数值矩阵的 c:

def mandelbrot(h,w, maxit=20):
    X,Y = meshgrid(linspace(-2, 0.8, w), linspace(-1.4, 1.4, h))
    c = X + Y*1j
    z = c
    exceeds = zeros(z.shape, dtype=bool)

    for iteration in range(maxit):
        z  = z**2 + c
        exceeded = abs(z) > 4
        exceeds_now = exceeded & (logical_not(exceeds))  
        exceeds[exceeds_now] = True        
        z[exceeded] = 2  # limit the values to avoid overflow
    return exceeds

imshow(mandelbrot(400,400),cmap='gray')
axis('off')

命令imshow将矩阵显示为图像。所选的颜色图用白色显示了序列显示为无界的区域，用黑色显示了其他区域。这里我们使用axis('off')来关闭轴，因为这对于图像可能不是很有用。

图 6.9:使用 imshow 将矩阵可视化为图像的示例。

默认情况下，imshow使用插值使图像看起来更好。当矩阵很小时，这一点很明显。下图显示了使用:

imshow(mandelbrot(40,40),cmap='gray')

和

imshow(mandelbrot(40,40), interpolation='nearest', cmap='gray')

在第二个例子中，像素值被复制。

图 6.10:使用 imshow 的线性插值与使用最近邻插值的区别

有关使用 Python 处理和绘制图像的更多详细信息，请参考【30】。

Matplotlib 对象

到目前为止，我们已经使用了 matplotlib 的pyplot模块。这个模块让我们可以轻松地直接使用最重要的绘图命令。大多数情况下，我们感兴趣的是创建一个图形并立即显示它。但是，有时我们想生成一个图形，稍后应该通过改变它的一些属性来修改它。这要求我们以面向对象的方式处理图形对象。在本节中，我们将介绍一些修改图形的基本步骤。对于 Python 中更复杂的面向对象的绘图方法，您必须离开pyplot并直接进入matplotlib及其大量的文档。

轴对象

当创建一个以后应该修改的图时，我们需要引用一个图形和一个轴对象。为此，我们必须先创建一个图形，然后定义一些轴及其在图形中的位置。我们不应该忘记将这些对象分配给一个变量:

fig = figure()
ax = subplot(111)

根据subplot的使用，一个图形可以有多个轴对象。在第二步中，绘图与给定的轴对象相关联:

fig = figure(1)
ax = subplot(111)
x = linspace(0,2*pi,100) 
# We set up a function that modulates the amplitude of the sin function
amod_sin = lambda x: (1.-0.1*sin(25*x))*sin(x)
# and plot both...
ax.plot(x,sin(x),label = 'sin') 
ax.plot(x, amod_sin(x), label = 'modsin')

这里我们使用了一个由lambda关键字表示的匿名函数。我们将在后面的章节匿名函数中解释这个构造-在第 7 章、函数中的 lambda 关键字。事实上，这两个绘图命令用两个Lines2D对象填充列表ax.lines:

ax.lines #[<matplotlib.lines.Line2D at ...>, <matplotlib.lines.Line2D at ...>]

使用标签是一个很好的做法，这样我们以后就可以用一种简单的方式识别对象:

for il,line in enumerate(ax.lines):
    if line.get_label() == 'sin':
       break

我们现在以一种允许进一步修改的方式设置事情。我们目前得到的图如前图所示(图 6.11，左)。

修改线条属性

我们只是通过标签来识别特定的线对象。它是带有索引il的列表ax.lines列表的一个元素。它的所有属性都收集在字典中

dict_keys(['marker', 'markeredgewidth', 'data', 'clip_box', 'solid_capstyle', 'clip_on', 'rasterized', 'dash_capstyle', 'path', 'ydata', 'markeredgecolor', 'xdata', 'label', 'alpha', 'linestyle', 'antialiased', 'snap', 'transform', 'url', 'transformed_clip_path_and_affine', 'clip_path', 'path_effects', 'animated', 'contains', 'fillstyle', 'sketch_params', 'xydata', 'drawstyle', 'markersize', 'linewidth', 'figure', 'markerfacecolor', 'pickradius', 'agg_filter', 'dash_joinstyle', 'color', 'solid_joinstyle', 'picker', 'markevery', 'axes', 'children', 'gid', 'zorder', 'visible', 'markerfacecoloralt'])

这可以通过以下命令获得:

ax.lines[il].properties()

它们可以通过相应的 setter 方法进行更改。让我们更改正弦曲线的线条样式:

ax.lines[il].set_linestyle('-.')
ax.lines[il].set_linewidth(2)

我们甚至可以修改数据，如图所示:

ydata=ax.lines[il].get_ydata()
ydata[-1]=-0.5
ax.lines[il].set_ydata(ydata)

结果如下图*(图 6.11，右)*:

图 6.11:调幅正弦函数(左)和最后一个数据点损坏的曲线(右)。

注释

一个有用的轴方法是annotate。它在给定位置设置注释，并用箭头指向图形中的另一个位置。箭头可以在字典中被赋予属性:

annot1=ax.annotate('amplitude modulated\n curve', (2.1,1.0),(3.2,0.5),
       arrowprops={'width':2,'color':'k', 'connectionstyle':'arc3,rad=+0.5', 
                   'shrink':0.05},
       verticalalignment='bottom', horizontalalignment='left',fontsize=15, 
                   bbox={'facecolor':'gray', 'alpha':0.1, 'pad':10})
annot2=ax.annotate('corrupted data', (6.3,-0.5),(6.1,-1.1),
       arrowprops={'width':0.5,'color':'k','shrink':0.1},
       horizontalalignment='center', fontsize=12)

在上面的第一个注释示例中，箭头指向坐标为( 2.1，1.0 )的点，文本的左下方坐标为( 3.2，0.5 )。如果没有另外指定，坐标是在方便的数据坐标系中给出的，这是指用于生成图的数据。

此外，我们演示了由arrowprop字典指定的几个箭头属性。您可以通过shrink键缩放箭头。设置'shrink':0.05将箭头尺寸减小 5%，以与它所指向的曲线保持一段距离。您可以使用connectionstyle键让箭头跟随样条弧或赋予它其他形状。

文本属性甚至文本周围的边界框都可以通过额外的关键字参数来标注方法，参见下图(图 6.12，左侧):

实验注释有时需要移除我们想要拒绝的尝试。因此，我们将注释对象分配给一个变量，这允许我们通过其remove方法移除注释:

annot1.remove()

填充曲线之间的区域

填充是突出曲线之间差异的理想工具，例如预期数据之上的噪声、近似与精确函数等。

填充是通过轴方法完成的

ax.fill_between(x,y1,y2)

对于我们使用的下一个数字:

axf = ax.fill_between(x, sin(x), amod_sin(x), facecolor='gray')

where是一个非常方便的参数，需要一个布尔数组来指定额外的填充条件。

axf = ax.fill_between(x, sin(x), amod_sin(x),where=amod_sin(x)-sin(x) > 0, facecolor=’gray’)

选择要填充的区域的布尔数组是amod_sin(x)-sin(x) > 0。

下图显示了带有两种填充区域的曲线:

图 6.12:带有注释和填充区域的调幅正弦函数(左)和使用 where 参数仅部分填充区域的修改图形(右)。

如果您自己测试这些命令，在尝试部分填充之前，不要忘记删除完整的填充，否则您将看不到任何变化:

axf.remove()

相关填充命令为fill、fill_betweenx。

刻度和标签

如果谈话、海报和出版物中的人物没有过多的不必要信息，他们看起来会好得多。你想把观众引向那些包含信息的部分。在我们的示例中，我们通过从 x 轴和 y 轴移除刻度并引入与问题相关的刻度标签来清理图片:

图 6.13:振幅调制正弦函数的完整示例，带有注释和填充区域以及修改的刻度和刻度标签。

ax.set_xticks(array([0,pi/2,pi,3/2*pi,2*pi]))
ax.set_xticklabels(('$0$','$\pi/2$','$\pi$','$3/2 \pi$','$2 \pi$'),fontsize=18)
ax.set_yticks(array([-1.,0.,1]))
ax.set_yticklabels(('$-1$','$0$','$1$'),fontsize=18)

请注意，我们在字符串中使用 LaTeX 格式来表示希腊字母，正确设置公式，并使用 LaTeX 字体。增加字体大小也是一种很好的做法，这样可以将生成的图形缩小到文本文档中，而不会影响坐标轴的可读性。本指导例最终结果见上图(图 6.13 )。

制作三维图

有一些有用的matplotlib工具包和模块可以用于各种特殊用途。在本节中，我们描述了一种生成三维图的方法。

mplot3d工具包提供点、线、轮廓、表面和所有其他基本组件的三维绘图，以及三维旋转和缩放。如下例所示，通过将关键字projection='3d'添加到轴对象来制作三维图:

from mpl_toolkits.mplot3d import axes3d

fig = figure()
ax = fig.gca(projection='3d')
# plot points in 3D
class1 = 0.6 * random.standard_normal((200,3))
ax.plot(class1[:,0],class1[:,1],class1[:,2],'o')
class2 = 1.2 * random.standard_normal((200,3)) + array([5,4,0])
ax.plot(class2[:,0],class2[:,1],class2[:,2],'o')
class3 = 0.3 * random.standard_normal((200,3)) + array([0,3,2])
ax.plot(class3[:,0],class3[:,1],class3[:,2],'o')

可以看到，需要从mplot3d导入axes3D类型。生成的图显示了分散的三维数据，如下图所示(图 6.14 )

图 6.14:使用 mplot3d 工具包绘制三维数据

绘制曲面同样简单。以下示例使用内置函数get_test_data创建用于绘制曲面的样本数据。考虑以下带有透明度的曲面图示例。

X,Y,Z = axes3d.get_test_data(0.05)

fig = figure()
ax = fig.gca(projection='3d')
# surface plot with transparency 0.5 
ax.plot_surface(X,Y,Z,alpha=0.5)

α值设置透明度。曲面图如下图所示(图 6.15 )。

图 6.15:用三个 2D 投影绘制表面网格的示例。

您也可以在任何坐标投影中绘制等高线，如下例所示。

fig = figure()
ax = fig.gca(projection = '3d')
ax.plot_wireframe(X,Y,Z,rstride = 5,cstride = 5)

# plot contour projection on each axis plane
ax.contour(X,Y,Z, zdir='z',offset = -100)
ax.contour(X,Y,Z, zdir='x',offset = -40)
ax.contour(X,Y,Z, zdir='y',offset = 40)

# set axis limits
ax.set_xlim3d(-40,40)
ax.set_ylim3d(-40,40)
ax.set_zlim3d(-100,100)

# set labels
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z axis')

注意设置轴限制的命令。用于设置轴限制的标准matplotlib命令是axis([-40, 40, -40, 40])，这对于 2D 地块很有效。然而，axis([-40,40,-40,40,-40,40])不起作用。对于三维绘图，您需要使用面向对象版本的命令，ax.set_xlim3d(-40,40)和类似的命令。标注轴也是如此；请注意设置标签的命令。对于 2D 地块你可以做xlabel(’X axis’)和ylabel(’Y axis’)但是没有zlabel命令。相反，在三维绘图中，您需要使用ax.set_xlabel(’X axis’)和类似的工具，如前例所示。

这段代码的结果如下

有许多选项可用于设置地块外观的格式，包括曲面的颜色和透明度。mplot3d文档网站【23】有详细信息。

用剧情拍电影

如果您有演变的数据，除了在图形窗口中显示之外，您可能还想将其保存为电影，类似于savefig命令。一种方法是使用 visvisvis 提供的visvis模块(更多信息请参考【37】)。

这里有一个简单的例子，用隐式表示法来演化一个圆。让圆用函数的零电平来表示。或者，考虑零点集合内的圆盘。如果 f 的值以 v 的速率减小，则圆将以的速率向外移动。

这可以实现为:

import visvis.vvmovie as vv

# create initial function values
x = linspace(-255,255,511)
X,Y = meshgrid(x,x)
f = sqrt(X*X+Y*Y) - 40 #radius 40

# evolve and store in a list
imlist = []
for iteration in range(200):
    imlist.append((f>0)*255)
    f -= 1 # move outwards one pixel
vv.images2swf.writeSwf('circle_evolution.swf',imlist)

结果是一部 Flash 电影(*。swf 文件)的一个不断增长的黑色圆圈，如下图所示(图 6.16) :

图 6.16:一个演化圆的例子

在本例中，使用数组列表来创建电影。visvis模块还可以保存 GIF 动画，在某些平台上还可以保存 AVI 动画(。gif 和。avi 文件 ) ，也有可能直接从人物窗口捕捉电影画面。然而，这些选项需要在您的系统上安装更多的软件包(例如，PyOpenGL和PIL，Python 图像库)。详见visvis网页上的文档。

另一种选择是使用savefig创建图像，每帧一个。

# create initial function values
x = linspace(-255,255,511)
X,Y = meshgrid(x,x)
f = sqrt(X*X+Y*Y) - 40 #radius 40
for iteration in range(200):
    imshow((f>0)*255)
    gray()
    axis('off')
    savefig('circle_evolution_{:d}.png'.format(iteration))
    f -= 1

然后，可以使用标准的视频编辑软件(例如美柯德或 ImageMagick)来组合这些图像。这种方法的优点是，您可以通过保存高分辨率图像来制作高分辨率视频。

总结

图形表示是呈现数学结果或算法行为的最简洁形式。本章为您提供了绘图的基本工具，并向您介绍了以面向对象的方式处理图形对象(如图形、轴和线)的更复杂的方法。

在本章中，您学习了如何制作图，不仅是经典的 x/y 图，还有 3D 图和直方图。我们给了你一个拍电影的开胃菜。您还看到了如何修改绘图，将它们视为具有相关方法和属性的图形对象，这些方法和属性可以设置、删除或修改。

练习

Ex。1 →编写一个函数，给定椭圆的中心坐标( x，y )、半轴 a 和 b 旋转角度θ绘制椭圆。

Ex。2 →写一个短程序，取一个 2D 数组，比如前面的 Mandelbrot 轮廓图像，用相邻值的平均值迭代替换每个值。在图形窗口中更新数组的等高线图，以动画显示等高线的演变。解释行为。

Ex。3 →考虑一个带有整数值的 N × N 矩阵或图像。地图

是点的环形正方形网格映射到自身的一个示例。这有一个有趣的特性，它通过剪切扭曲图像，然后使用 modulu 函数mod将图像之外的部分移回。迭代应用，这导致随机图像的方式，最终返回原来的。执行以下顺序:

并将前 N 步保存到文件中或在图形窗口中绘制出来。

作为示例图像，您可以使用来自scipy.misc的经典 512 × 512 Lena 测试图像。

from scipy.misc import lena
I = lena()

结果应该如下所示:

|  |  | … |  | … |  | … |  |  | | Zero | one | | One hundred and twenty-eight | | Two hundred and fifty-six | | Five hundred and eleven | Five hundred and twelve |

型式

计算 x 和 y 映射，并使用数组索引(参考第 5 章、高级数组概念中的数组索引部分)来复制像素值。

Ex。4 →图像读取和绘制。SciPy 自带imread功能(在scipy.misc模块)读取图像，(参考第 12 章、输入输出中读写图像一节)。写一个简短的程序，从文件中读取一个图像，并在给定的灰度值下绘制覆盖在原始图像上的图像轮廓。

型式

您可以通过对颜色通道进行平均来获得图像的灰度版本，如下所示:mean(im,axis=2)

Ex。5 →图像边缘。2D 拉普拉斯的过零点是图像边缘的良好指示。修改上一练习中的程序，使用scipy.ndimage模块中的gaussian_laplace或laplace功能计算 2D 拉普拉斯算子，并将边缘覆盖在图像顶部。

Ex。6 →使用orgid代替meshgrid,重新表述曼德罗德分形示例(参见图像和轮廓部分)另请参见第 5 章、高级数组概念中双变量函数中的解释ogrid。orgid、mgrid和meshgrid有什么区别？

七、函数

本章介绍函数，这是编程中的一个基本构件。我们展示了如何定义它们，如何处理输入和输出，如何正确使用它们，以及如何将它们视为对象。

基础

在数学中，一个函数被写成一个映射，它唯一地将一个元素 y 从范围 R 分配给域 D 中的每个元素 x 。

这由 f : D → R 表示

或者，当考虑特定元素 x 和 y 时，可以写 f : x → y

在这里， f 被称为函数的名称， f(x) 是其应用于 x 时的值。在这里， x 有时被称为 f. 我们先看一个例子，然后再考虑 Python 中的函数。

例如， D = ℝ x ℝ和 y = f(x ₁ ，x₂)= x₁-x₂。这个函数把两个实数映射到它们的差上。

在数学中，函数可以有数字、向量、矩阵，甚至其他函数作为自变量。下面是一个带有混合参数的函数示例:

。

在这种情况下，会返回一个数字。使用函数时，我们必须区分两个不同的步骤:

• 函数的定义
• 函数的求值，即对于给定值 x 的 f(x) 的计算

第一步完成一次，而第二步可以针对各种参数执行多次。编程语言中的函数遵循相同的概念，并将其应用于各种类型的输入参数，例如字符串、列表或任何对象。我们通过再次考虑给定的例子来演示函数的定义:

def subtract(x1, x2):
    return x1 - x2

关键字def表示我们要定义一个函数。subtract是函数的名称，x1和 *x2* 是它的参数。冒号表示我们正在使用一个 block 命令，函数返回的值跟在return关键字后面。现在，我们可以评估这个函数。调用此函数时，其参数由输入参数代替:

r = subtract(5.0, 4.3)

计算结果 0.7 并分配给r变量。

参数和参数

定义函数时，其输入变量称为函数的参数。执行函数时使用的输入称为其参数。

传递参数-按位置和关键字

我们将再次考虑前面的例子，其中函数采用两个参数，即x1和x2。

它们的名字用来区分这两个数字，在这种情况下，不改变结果就不能互换。第一个参数定义从中减去第二个参数的数字。调用subtract时，每个参数都被一个参数代替。只有争论的顺序很重要；参数可以是任何对象。例如，我们可以称之为:

z = 3 
e = subtract(5,z)

除了这种调用函数的标准方式(即通过位置传递参数)之外，有时使用关键字传递参数可能会更方便。参数的名称是关键字；考虑以下实例:

z = 3 
e = subtract(x2 = z, x1 = 5)

这里，参数是按名称而不是按调用中的位置分配给参数的。调用函数的两种方式可以结合起来，这样由 position 给出的参数排在第一位，由 keyword 给出的参数排在最后。我们使用功能plot来显示这一点，该功能在第 6 章、标绘中有描述:

plot(xp, yp, linewidth = 2,label = 'y-values')

改变论点

参数的目的是为函数提供必要的输入数据。更改函数内部的参数值通常不会影响函数外部的参数值:

def subtract(x1, x2):
    z = x1 - x2
    x2 = 50.
    return z
a = 20.
b = subtract(10, a)    # returns -10
a    # still has the value 20

这适用于所有不可变的参数，如字符串、数字和元组。如果可变参数(如列表或字典)被更改，情况就不同了。

例如，将可变输入参数传递给函数，并在函数内部更改它们，也可以在函数外部更改它们:

def subtract(x):
    z = x[0] - x[1]
    x[1] = 50.
    return z
a = [10,20]
b = subtract(a)    # returns -10
a    # is now [10, 50.0]

这样的函数滥用它的参数来返回结果。我们强烈建议您不要这样构造，并且建议您不要更改函数内部的输入参数(更多信息请参考默认参数部分)。

访问在本地命名空间之外定义的变量

Python 允许函数访问在其任何封闭程序单元中定义的变量。与局部变量相反，这些变量被称为全局变量。后者只能在函数中访问。例如，考虑以下代码:

import numpy as np # here the variable np is defined
def sqrt(x):
    return np.sqrt(x) # we use np inside the function

这个功能不应该被滥用。下面的代码是一个不应该做什么的例子:

a = 3
def multiply(x):
    return a * x # bad style: access to the variable a defined outside

当改变变量a时，函数multiply默认改变其行为:

a=3
multiply(4)  # returns 12
a=4
multiply(4)  # returns 16

在这种情况下，最好通过参数列表提供变量作为参数:

def multiply(x, a):
    return a * x

使用闭包时，全局变量可能很有用。名称空间和范围在第 11 章、名称空间、范围和模块中有更广泛的讨论。

默认参数

有些函数可以有许多参数，其中一些可能只在非标准情况下有用。如果参数可以自动设置为标准(默认)值，这将是切实可行的。我们通过查看scipy.linalg模块中的命令norm来演示默认参数的使用。它计算矩阵和向量的各种范数。

以下片段要求计算 3 × 3 单位矩阵的弗罗贝纽斯范数是等价的(关于矩阵范数的更多信息可以在【10】中找到):

import scipy.linalg as sl
sl.norm(identity(3))
sl.norm(identity(3), ord = 'fro')
sl.norm(identity(3), 'fro')

注意，在第一次调用中，没有给出关于ord关键字的信息。Python 如何知道它应该计算 Frobenius 范数而不是另一个范数，例如欧几里德 2-范数？

上一个问题的答案是使用默认值。默认值是函数定义已经给出的值。如果调用函数时没有提供这个参数，Python 会使用程序员在定义函数时提供的值。

假设我们调用的函数subtract只有一个参数；我们会收到一条错误消息:

TypeError: subtract() takes exactly 2 arguments (1 given)

为了允许省略参数x2，函数的定义必须提供默认值，例如:

def subtract(x1, x2 = 0): 
    return x1 - x2

总而言之，参数可以作为位置参数和关键字参数给出。必须首先给出所有的位置参数。只要省略的参数在函数定义中有默认值，就不需要提供所有关键字参数。

小心可变默认参数

默认参数是根据函数定义设置的。使用默认值时，更改函数内部的可变参数会产生副作用，例如:

def my_list(x1, x2 = []):
    x2.append(x1)
    return x2
my_list(1)  # returns [1]
my_list(2)  # returns [1,2]

可变参数数

列表和字典可用于定义或调用具有可变数量参数的函数。让我们定义一个列表和一个字典如下:

data = [[1,2],[3,4]]    
style = dict({'linewidth':3,'marker':'o','color':'green'})

然后我们可以使用带星号(*)的参数调用plot函数:

plot(*data,**style)

以*为前缀的变量名，如前面例子中的*data，意味着提供了一个在函数调用中被解包的列表。这样，列表会生成位置参数。类似地，以**为前缀的变量名，如示例中的**style，将字典解包为关键字参数。参考下图(图 7.1 ):

图 7.1:函数调用中的星号参数

您可能还想使用相反的过程，其中所有给定的位置参数都打包成一个列表，所有关键字参数在传递给函数时都打包成一个字典。

在函数定义中，这由分别以*和**为前缀的参数表示。您会经常在代码文档中找到*args和**kwargs参数，参见图 7.2 。

图 7.2:函数定义中的星形参数

返回值

Python 中的函数总是返回单个对象。如果一个函数必须返回多个对象，这些对象将被打包并作为单元组对象返回。

例如，下面的函数取一个复数 z ，并根据欧拉公式将其极坐标表示返回为幅度 r 和角度:

Python 的对应物是这样的:

def complex_to_polar(z):
    r = sqrt(z.real ** 2 + z.imag ** 2)
    phi = arctan2(z.imag, z.real)
    return (r,phi)  # here the return object is formed

这里，我们使用sqrt(x) NumPy 函数作为数字x的平方根，使用arctan2(x,y)作为表达式 tan ^-1 ( x/y )。

让我们试试我们的功能:

z = 3 + 5j  # here we define a complex number
a = complex_to_polar(z)
r = a[0]
phi = a[1]

最后三个语句可以用一行写得更优雅:

r,phi = complex_to_polar(z)

我们可以通过调用polar_to_comp来测试我们的功能；参考练习 1 。

如果一个函数没有return语句，则返回值None。在许多情况下，函数不需要返回值。这可能是因为传递给函数的变量可能会被修改。例如，考虑以下函数:

def append_to_list(L, x):
    L.append(x)

前面的函数不返回任何内容，因为它修改了作为参数给出的一个对象。我们在参数和参数部分提到了为什么这是有用的。有许多方法的行为方式是相同的。仅提及列表方法，append、extend、reverse和sort方法不返回任何内容(即返回None)。当一个对象被一个方法以这种方式修改时，这种修改被就地调用。除非通过查看代码或文档，否则很难知道一个方法是否改变了一个对象。

函数或方法不返回任何内容的另一个原因是当它打印出消息或写入文件时。

执行在第一个出现的return语句处停止。该语句后的行是死代码，永远不会执行:

def function_with_dead_code(x):
    return 2 * x
    y = x ** 2 # these two lines ...
    return y   # ... are never executed!

递归函数

在数学中，许多函数都是递归定义的。在这一节中，我们将展示这个概念如何在编程函数时使用。这使得程序与其数学对应物的关系非常清楚，这可能会降低程序的可读性。

然而，我们建议您谨慎使用这种编程技术，尤其是在科学计算领域。在大多数应用中，更直接的迭代方法更有效。从下面的例子中，这一点将立即变得清楚。

切比雪夫多项式由三项递归定义:

这样的递归需要初始化，即 T ₀ ( x ) =1、TT8】1(x)=x .

在 Python 中，这三项递归可以通过以下函数定义来实现:

def chebyshev(n, x):
    if n == 0:
        return 1.
    elif n == 1:
        return x
    else:
        return 2\. * x * chebyshev(n - 1, x) 
                      - chebyshev(n - 2 ,x)

该函数的调用方式如下:

chebyshev(5, 0.52) # returns 0.39616645119999994

这个例子也说明了大幅浪费计算时间的风险。函数求值的数量随着递归级别呈指数增长，并且这些求值中的大多数只是以前计算的重复。虽然使用递归程序来演示代码和数学定义之间的强关系可能很有诱惑力，但是生产代码将避免这种编程技术(参考练习 6)。我们还提到了一种叫做记忆化的技术(更多细节请参考【22】)，它将递归编程与缓存技术相结合，以保存复制的函数求值。

递归函数通常有一个级别参数。在前面的例子中，它是 n. 它用于控制功能的两个主要部分:

• 基本情况，这里，前两个if分支
• 递归体，其中使用较小的级别参数调用函数本身一次或多次。

递归函数的执行传递的级别数称为递归深度。这个数量不能太大；否则计算可能不再有效，最终会出现以下错误:

RuntimeError: maximum recursion depth exceeded

最大递归深度取决于您使用的计算机的内存。当函数定义中缺少初始化步骤时，也会出现此错误。我们鼓励对非常小的递归深度使用递归程序(有关更多信息，请参考第 9 章、迭代的无限迭代一节。

功能文档

您应该在开始时使用字符串来记录您的函数。这称为文档字符串:

def newton(f, x0):
    """
    Newton's method for computing a zero of a function
    on input:
    f  (function) given function f(x)
    x0 (float) initial guess 
    on return:
    y  (float) the approximated zero of f
    """
     ...

调用help(newton)时，这个文档字符串会和这个函数的调用一起显示出来:

Help on function newton in module __main__:

newton(f, x0)
     Newton's method for computing a zero of a function
     on input:
     f  (function) given function f(x)
     x0 (float) initial guess
     on return:
     y  (float) the approximated zero of f

文档字符串在内部保存为给定函数的属性__doc__。在示例中，它是newton.__doc__。您应该在文档字符串中提供的最小信息是函数的目的以及输入和输出对象的描述。有工具可以通过收集程序中的所有文档字符串来自动生成完整的代码文档(更多信息请参考 [32] )。

函数是对象

函数是对象，就像 Python 中的其他东西一样。人们可以传递函数作为参数，更改它们的名称，或者删除它们。例如:

def square(x):
    """
    Return the square of x
    """
    return x ** 2
square(4) # 16
sq = square # now sq is the same as square
sq(4) # 16
del square # square doesn't exist anymore
print(newton(sq, .2)) # passing as argument

在科学计算中应用算法时，将函数作为参数传递是非常常见的。计算给定函数的零点的scipy.optimize中的函数fsolve或计算积分的scipy.integrate中的quad就是典型的例子。

函数本身可以有不同类型的不同数量的参数。因此，当将您的函数f作为参数传递给另一个函数g时，请确保f与g的文档字符串中描述的形式完全相同。

fsolve的文档字符串给出了其func参数的信息:

func -- A Python function or method which takes at least one
               (possibly vector) argument.

部分应用

让我们从一个双变量函数的例子开始。

函数可以看作是两个变量的函数。人们通常认为ω不是自由变量，而是函数族的固定参数:

这种解释将两个变量中的一个函数简化为一个变量中的一个函数t，给定一个固定的参数值 ω 。通过固定(冻结)函数的一个或几个参数来定义新函数的过程称为部分应用。

使用 Python 模块functools可以轻松创建部分应用，该模块正是为此目的提供了一个名为partial的函数。我们通过构造一个函数来说明这一点，该函数返回给定频率的正弦值:

import functools
def sin_omega(t, freq):
    return sin(2 * pi * freq * t)

def make_sine(frequency):
    return functools.partial(sin_omega, freq = frequency)

使用闭包

使用函数是对象的观点，部分应用可以通过编写一个函数来实现，这个函数本身返回一个新的函数，减少了输入参数的数量。例如，该函数可以定义如下:

def make_sine(freq):
    "Make a sine function with frequency freq"
    def mysine(t):
        return sin_omega(t, freq)
    return mysine

在本例中，内部函数mysine可以访问变量freq；它既不是这个函数的局部变量，也不是通过参数列表传递给它的。Python 允许这样的构造(参见第 11 章、命名空间、范围和模块中的命名空间部分)。

匿名函数 lambda 关键字

Python 中使用关键字 lambda 定义匿名函数，即；函数没有名字，用一个表达式描述。您可能只想对一个可以用简单表达式表示的函数执行操作，而不需要命名该函数，也不需要用冗长的def块定义该函数。

注

名称λ源于微积分和数理逻辑的一个特殊分支，-微积分。

例如，为了计算下面的表达式，我们可以使用 SciPy 的函数quad，它要求函数被积分作为它的第一个参数，积分界限作为接下来的两个参数:

这里，要集成的函数只是一个简单的单行函数，我们使用lambda关键字来定义它:

import scipy.integrate as si
si.quad(lambda x: x ** 2 + 5, 0, 1)

语法如下:

lambda parameter_list: expression

lambda函数的定义只能由一个表达式组成，特别是不能包含循环。lambda与其他函数一样，函数是对象，可以分配给变量:

parabola = lambda x: x ** 2 + 5
parabola(3) # gives 14

λ结构总是可替换的

需要注意的是，lambda 构造只是 Python 中的语法糖。任何 lambda 结构都可以用显式函数定义来代替:

parabola = lambda x: x**2+5 
# the following code is equivalent
def parabola(x):
    return x ** 2 + 5

使用构造的主要原因是为了非常简单的函数，而完整的函数定义会过于繁琐。

lambda函数提供了第三种构造闭包的方法，正如我们继续前面的例子所展示的:

我们使用sin_omega函数计算不同频率下正弦函数的积分:

import scipy.integrate as si
for iteration in range(3):
    print(si.quad(lambda x: sin_omega(x, iteration*pi), 0, pi/2.) )

作为装饰者的功能

在部分应用部分，我们看到了如何使用一个函数来修改另一个函数。装饰器是 Python 中的语法元素，它方便地允许我们改变函数的行为，而不改变函数本身的定义。让我们从以下情况开始:

假设我们有一个决定矩阵稀疏程度的函数:

def how_sparse(A):
    return len(A.reshape(-1).nonzero()[0])

如果该函数不是以数组对象作为输入调用的，它将返回一个错误。更准确地说，它不适用于不实现reshape方法的对象。例如，how_sparse函数不适用于列表，因为列表没有reshape方法。下面的 helper 函数修改任何带有一个输入参数的函数，以便尝试对数组进行类型转换:

def cast2array(f):
    def new_function(obj):
        fA = f(array(obj))
        return fA
    return new_function

因此，修改后的函数how_sparse = cast2array(how_sparse)可以应用于任何可以转换为数组的对象。如果用类型转换功能修饰how_sparse的定义，也可以实现同样的功能。还建议考虑functools.wraps(详见【8】):

@cast2array
def how_sparse(A):
    return len(A.reshape(-1).nonzero()[0])

要定义装饰器，您需要一个可调用的对象，例如修改要装饰的函数定义的函数。主要目的是:

• 通过从不直接为其功能服务的函数中分离出部分来提高代码的可读性(例如，记忆)
• 将一系列相似函数的公共序言和结尾部分放在一个公共位置(例如，类型检查)
• 能够轻松地关闭和打开功能的附加功能(例如，测试打印、跟踪)

总结

函数不仅是使程序模块化的理想工具，而且也反映了数学思维。您学习了函数定义的语法，以及如何区分定义和调用函数。

我们认为函数是可以被其他函数修改的对象。使用函数时，熟悉变量范围的概念以及信息如何通过参数传递到函数中是很重要的。

有时，用所谓的匿名函数动态定义函数是很方便的。为此，我们引入了关键字 lambda。

练习

例 1 →编写一个函数polar_to_comp，该函数接受两个参数 r 和并返回复数指数函数使用 NumPy 函数exp。

Ex 2 →在 Python 模块functools的描述中(参见【8】了解更多关于 functools 的详细信息)您可以找到以下 Python 函数:

def partial(func, *args, **keywords):
    def newfunc(*fargs, **fkeywords):
        newkeywords = keywords.copy()
        newkeywords.update(fkeywords)
        return func(*(args + fargs), **newkeywords)
    newfunc.func = func
    newfunc.args = args
    newfunc.keywords = keywords
    return newfunc

解释并测试该功能。

Ex 3 →为函数how_sparse编写装饰器，通过将小于 1.e-16 的元素设置为零来清理输入矩阵 A (考虑第节中作为装饰器的示例)。

Ex 4 → A 连续函数 f 与f(A)f(b)<0 在区间【 a，b 中改变其符号，并且在该区间中至少有一个根(零)。这样的根可以用等分法找到。此方法从给定的时间间隔开始。然后研究子区间中的符号变化，

如果符号在第一个子区间 b 中发生变化，则重新定义为:

否则，它会以相同的方式重新定义为:

并且重复该过程直到 b-a 小于给定的公差。

• 将此方法实现为将作为参数的函数:
  • –功能 f
  • –初始间隔[ a，b ]
  • –公差
• 该功能bisec应返回最终间隔及其中点。
• 在区间[1.1，1.4] *、*以及在[1.3，1.4]中交替使用函数arctan和多项式 f(x) = 3 x ² -5 测试该方法。

Ex。5 →两个整数的最大公约数可以用如下递归描述的欧几里德算法计算:

写一个计算两个整数的最大公约数的函数。编写另一个函数，使用以下关系计算这些数字的最小公倍数:

Ex。6 →研究切比雪夫多项式的递归实现，考虑递归函数部分的例子。以非递归方式重写程序，研究计算时间与多项式次数的关系(另见timeit模块)。

八、类

在数学中，当我们写罪的时候，我们指的是一个数学对象，我们从初等微积分中知道许多方法。例如:

• 我们可能希望在 x= 0.5 处计算 sin x ，即计算 sin(0.5)，它返回一个实数
• 我们可能想计算它的导数，这给了我们另一个数学对象，cos
• 我们可能需要计算泰勒多项式的前三个系数

这些方法不仅可以应用于 sin，还可以应用于其他足够平滑的函数。然而，还有其他的数学对象(例如数字 5 )这些方法对它们没有意义。具有相同方法的对象被组合在抽象类中，例如函数。每一个可以应用于函数的语句和方法都特别适用于 s in 或 cos。这类的其他例子可能是有理数，它有分母和分子方法；一个区间，它有左右边界法；一个无穷序列，我们可以问它是否有极限，等等。

在这种情况下，sin 被称为类的实例。数学短语*让 g 是一个函数...*在本文中称为实例化。这里， g 是函数的名称；可以分配给它的许多属性之一。另一个属性可能是它的域。

数学对象 p(x) = 2x ² - 5 就像正弦函数一样。每种功能方法都适用于 p ，但是我们也可以为 p 定义特殊的方法。例如，我们可能会要求 p 的系数。这些方法可以用来定义多项式的种类。由于多项式是函数，它们还继承了函数类的所有方法。

在数学中，我们经常对完全不同的运算使用同一个运算符符号。例如，在 5+4 和 sin + cos 中，运算符符号+有不同的含义。通过使用相同的符号，人们试图表达数学运算的相似性。我们通过将这些术语应用于数学示例，从面向对象编程中引入了这些术语:

• 类
• 实例和实例化
• 遗产
• 方法
• 属性
• 操作员超载

在本章中，我们将展示如何在 Python 中使用这些概念。

课程介绍

我们将用有理数的例子来说明类的概念，即形式为q = q_Nq_D的数，其中qT10】N 和qT14】D 都是整数。

图 8.1:一个类声明的例子

我们这里使用有理数只是作为类概念的一个例子。关于 Python 中有理数的未来工作，请使用分数模块(请参考【6】)。

类语法

一个类的定义是由一个带有class关键字的块命令、类的名称和块中的一些语句组成的(参见图 8.1 ):

class RationalNumber: 
      pass

这个类的一个实例(或者换句话说，类型为RationalNumber的对象)是由

r = RationalNumber()

查询type(a)返回答案<class'__main__.RationalNumber'>。如果我们想调查一个对象是否是这个类的实例，我们可以使用这个:

if isinstance(a, RationalNumber):
    print('Indeed it belongs to the class RationalNumber')  

到目前为止，我们已经生成了一个RationalNumber类型的对象，它还没有数据。此外，没有定义对这些对象执行操作的方法。这将是下一部分的主题。

_ _ init _ _ 方法

现在我们为示例类提供一些属性；也就是说，我们给它定义数据。在我们的例子中，这个数据是分母和分子的值。为此，我们必须定义一个方法，__init__，用于用这些值初始化类:

class RationalNumber:
    def __init__(self, numerator, denominator):
        self.numerator = numerator
        self.denominator = denominator

在解释我们添加到类中的特殊__init__函数之前，我们演示一个RationalNumber对象的实例化:

q = RationalNumber(10, 20)    # Defines a new object
q.numerator    # returns 10
q.denominator    # returns 20

使用类名创建一个类型为RationalNumber的新对象，就像它是一个函数一样。这种说法做了两件事:

• 它首先创建一个空对象q。
• 然后对其应用__init__功能；即执行q.__init__(10, 20)。

__init__的第一个参数是指新对象本身。在函数调用时，第一个参数被对象的实例替换。这适用于该类的所有方法，而不仅仅适用于特殊方法__init__。第一个参数的特殊作用体现在命名为self的惯例上。在前面的例子中，__init__函数定义了新对象的两个属性，numerator和denominator。

属性和方法

使用类的主要原因之一是对象可以组合在一起并绑定到一个公共对象。我们在看有理数时已经看到了这一点；分母和分子是我们绑定到RationalNumber类实例的两个对象。它们被称为实例的属性。一个对象是一个类实例的一个属性的事实从它们被引用的方式中变得显而易见，我们之前已经默认使用了这一点:

<object>.attribute

以下是实例化和属性引用的一些示例:

q = RationalNumber(3, 5) # instantiation
q.numerator     # attribute access
q.denominator

a = array([1, 2])    # instantiation
a.shape

z = 5 + 4j    # instantiation
z.imag

一旦定义了一个实例，我们就可以设置、更改或删除该特定实例的属性。语法与常规变量相同:

q = RationalNumber(3, 5) 
r = RationalNumber(7, 3)
q.numerator = 17
del r.denominator

更改或删除属性可能会产生不希望的副作用，甚至会使对象变得无用。我们将在相互依赖的属性一节中了解更多信息。由于函数也是对象，我们也可以将函数用作属性；它们被称为实例的方法:

<object>.method(<arguments...>)

例如，让我们向类RationalNumber添加一个方法，将数字转换为浮点数:

class RationalNumber:
...
    def convert2float(self):
        return float(self.numerator) / float(self.denominator)

同样，该方法将对象本身的引用self作为其第一个(也是唯一一个)参数。我们将此方法用于常规函数调用:

q = RationalNumber(10, 20)    # Defines a new object
q.convert2float() # returns 0.5   

这相当于以下调用:

RationalNumber.convert2float(q)

再次注意，对象实例是作为函数的第一个参数插入的。第一个参数的使用解释了如果该特定方法与其他参数一起使用时会出现的错误消息:

q.convert2float(15)调用会引发以下错误信息:

TypeError: convert2float() takes exactly 1 argument (2 given)

这不起作用的原因是q.convert2float(15)精确地等同于RationalNumber.convert2float(q,15)，这是失败的，因为RationalNumber.convert2float只接受一个参数。

特殊方法

特殊的方法__repr__让我们能够定义对象在 Python 解释器中的表示方式。对于有理数，该方法的可能定义如下:

class RationalNumber:
 ...
    def __repr__(self):
        return '{} / {}'.format(self.numerator,self.denominator)

定义了这个方法，只需输入q返回 10 / 20。

我们希望有一种方法可以将两个有理数相加。第一次尝试可能会产生这样的方法:

class RationalNumber:
...
    def add(self, other):
        p1, q1 = self.numerator, self.denominator
        if isinstance(other, int):
            p2, q2 = other, 1
        else:
            p2, q2 = other.numerator, other.denominator
        return RationalNumber(p1 * q2 + p2 * q1, q1 * q2)

对此方法的调用采用以下形式:

q = RationalNumber(1, 2)
p = RationalNumber(1, 3)
q.add(p)   # returns the RationalNumber for 5/6

如果我们能写q + p就更好了。但到目前为止，RationalNumber类型没有定义加号。这是通过使用__add__特殊方法完成的。因此，仅仅将add重命名为__add__就允许使用有理数的加号:

q = RationalNumber(1, 2)
p = RationalNumber(1, 3)
q + p # RationalNumber(5, 6)

表达式q + p实际上是表达式q.__add__(p)的别名。在表(表 8.1) 中，你会发现二元运算符的特殊方法，如+、-或*。

| **操作员** | **方法** | **操作员** | **方法** | | `+` | `__add__` | `+=` | `__iadd__` | | `*` | `__mul__` | `*=` | `__imul__` | | `-` | `__sub__` | `-=` | `__isub__` | | `/` | `__truediv__` | `/=` | `__itruediv__` | | `//` | `__floordiv__` | `//=` | `__ifloordiv__` | | `**` | `__pow__` | | | | `==` | `__eq__` | `!=` | `__ne__` | | `<=` | `__le__` | `<` | `__lt__` | | `>=` | `__ge__` | `>` | `__gt__` | | `()` | `__call__` | `[]` | `__getitem__` |

表 8.1:部分 Python 运算符及对应的类方法，可以找到完整的列表【31】

新类的这些运算符的实现称为运算符重载。运算符重载的另一个例子是检查两个有理数是否相同的方法:

class RationalNumber:
...
    def __eq__(self, other):
        return self.denominator * other.numerator == 
            self.numerator * other.denominator

它是这样使用的:

p = RationalNumber(1, 2) # instantiation
q = RationalNumber(2, 4) # instantiation
p == q # True

属于不同类的对象之间的操作需要特别小心:

p = RationalNumber(1, 2) # instantiation
p + 5  # corresponds to p.__add__(5)  
5 + p  # returns an error

默认情况下，+运算符调用左操作数的方法__add__。我们对它进行了编程，使得它既允许int类型的对象，也允许RationalNumber类型的对象。在语句5 + p中，操作数被转换，并调用内置int类型的__add__方法。这个方法返回一个错误，因为它不知道如何处理有理数。这个案子可以用__radd__的方法处理，我们现在就用这个方法装备RationalNumber班。方法__radd__叫做反向加法。

反向操作

如果像+这样的操作应用于两个不同类型的操作数，则首先调用左操作数的对应方法(在本例中为__add__)。如果这引发异常，则调用右操作数的反向方法(此处为__radd__)。如果该方法不存在，则会引发TypeError异常。

考虑一个反向操作的例子。为了启用操作 5+ p ，其中 p 是RationalNumber的一个实例，我们定义如下:

class RationalNumber:
   ....
    def __radd__(self, other):
        return self + other

注意__radd__交换参数的顺序；self是RationalNumber类型的对象，而其他是必须转换的对象。

将类实例与方括号*、* ( 、)或[，]一起使用时，会调用一个特殊方法__call__或__getitem__中的一个方法，为该实例提供一个函数或一个可迭代对象的行为(请参考 T able 8.1 中的这些和其他特殊方法):

class Polynomial:
...
    def __call__(self, x):
        return self.eval(x)

现在可以如下使用:

p = Polynomial(...)
p(3.) # value of p at 3.

如果类提供了迭代器，__getitem__特殊方法就有意义了(在考虑下面的例子之前，建议参考第 9 章、迭代器一节迭代器)。

这个递归u_I+1= a₁T6】u_IT9】+aT12】0T14】u_IT18】-1 被称为一个三 - 术语递归。它在应用数学，特别是正交多项式的构造中起着重要的作用。我们可以通过以下方式将三项递归设置为一个类:

import itertools

class  Recursion3Term:
    def __init__(self, a0, a1, u0, u1):
        self.coeff = [a1, a0]
        self.initial = [u1, u0]
    def __iter__(self):
        u1, u0 = self.initial
        yield u0  # (see also Iterators section in Chapter 9) 
        yield u1
        a1, a0 = self.coeff
        while True :
            u1, u0 = a1 * u1 + a0 * u0, u1
            yield u1
    def __getitem__(self, k):
        return list(itertools.islice(self, k, k + 1))[0]

这里，__iter__方法定义了一个生成器对象，它允许我们使用类的一个实例作为迭代器:

r3 = Recursion3Term(-0.35, 1.2, 1, 1)
for i, r in enumerate(r3):
    if i == 7:
        print(r)  # returns 0.194167
        break

__getitem__方法使我们能够直接访问迭代，就像r3是一个列表一样:

r3[7] # returns 0.194167

请注意，我们在编码__getitem__方法时使用了itertools.islice(更多信息，请参考第 9 章、迭代的迭代器部分)。在第 5 章、高级数组概念的双变量函数一节中给出了一个使用__getitem__和切片以及函数ogrid的例子。

相互依赖的属性

只需为实例的属性赋值，就可以更改(或创建)它们。但是，如果其他属性依赖于刚刚更改的属性，则最好同时更改这些属性:

让我们考虑一个从三个给定点为平面三角形定义一个对象的类。建立这样一个类的第一次尝试如下:

class Triangle:
    def __init__(self,  A, B, C):
        self.A = array(A)
        self.B = array(B)
        self.C = array(C)
        self.a = self.C - self.B
        self.b = self.C - self.A
        self.c = self.B - self.A
    def area(self):
        return abs(cross(self.b, self.c)) / 2

这个三角形的一个实例是这样创建的:

tr = Triangle([0., 0.], [1., 0.], [0., 1.])

它的面积是这样计算的:

tr.area() # returns 0.5

如果我们改变一个属性，比如点 B ，对应的边 a 和 c 不会自动更新，计算的面积是错误的:

tr.B = [12., 0.]
tr.area() # still returns 0.5, should be 6 instead.

补救方法是定义一个在属性改变时执行的方法；这样的方法称为 setter 方法。相应地，人们可能会要求在请求属性值时执行的方法；这种方法称为 getter 方法。

属性函数

函数property将一个属性链接到这样的 getter、setter 和 deleter 方法。它也可以用于将文档字符串分配给属性:

attribute = property(fget = get_attr, fset = set_attr, 
                     fdel = del_attr, doc = string)

我们用 setter 方法继续前面的例子，并再次考虑Trinagle类。如果以下陈述包含在Triangle中

B = property(fget = get_B, fset = set_B, fdel = del_B, doc = ’The point B of a triangle’)

命令

tr.B = <something>

调用 setter 方法set_B。

让我们修改三角形类:

class Triangle:
    def __init__(self, A, B, C):
        self._A = array(A)
        self._B = array(B)
        self._C = array(C)
        self._a = self._C - self._B
        self._b = self._C - self._A
        self._c = self._B - self._A
    def area(self):
        return abs(cross(self._c, self._b)) / 2.
    def set_B(self, B):
        self._B = B
        self._a = self._C - self._B
        self._c = self._B - self._A
    def get_B(self):
        return self._B
    def del_Pt(self):
        raise Exception('A triangle point cannot be deleted')
    B = property(fget = get_B, fset = set_B, fdel = del_Pt)

如果属性B改变，则方法set_B将新值存储在内部属性_B中，并改变所有相关属性:

tr.B = [12., 0.]
tr.area() # returns 6.0

这里使用deleter方法的方式是为了防止删除属性:

del tr.B # raises an exception

使用下划线作为属性名称的前缀是一种约定，用于指示不是为直接访问而设计的属性。它们旨在保存由设置器和获取器处理的属性的数据。这些属性在其他编程语言的意义上不是私有的；它们只是不打算被直接访问。

结合和非结合方法

我们现在将仔细看看属于方法的属性。让我们考虑一个例子:

class A:
    def func(self,arg):
        pass

一个小小的检查向我们展示了func的性质在创建一个实例后是如何变化的:

A.func  # <unbound method A.func>
instA = A()  # we create an instance
instA.func  #  <bound method A.func of ... >

例如，调用A.func(3)会导致如下错误消息:

TypeError: func() missing 1 required positional argument: 'arg'

instA.func(3)按预期执行。创建实例时，func方法绑定到该实例。self参数将实例赋值。将方法绑定到实例使该方法可用作函数。在此之前，是没有用的。类方法，我们稍后会考虑，在这方面是不同的。

类属性

类声明中指定的属性称为类属性。考虑以下示例:

class Newton:
    tol = 1e-8 # this is a class attribute
    def __init__(self,f):
        self.f = f # this is not a class attribute
    ....

类属性对于模拟默认值很有用，并且可以在必须重置值时使用:

N1 = Newton(f)
N2 = Newton(g)

两个实例都有一个属性tol，其值在类定义中初始化:

N1.tol # 1e-8
N2.tol # 1e-8

更改类属性会自动影响所有实例的所有相应属性:

Newton.tol = 1e-10
N1.tol # 1e-10
N2.tol # 1e-10

更改一个实例的tol不会影响另一个实例:

N2.tol = 1.e-4
N1.tol  # still 1.e-10

但是现在N2.tol脱离了类属性。更改Newton.tol对N2.tol不再有任何影响；

Newton.tol = 1e-5 # now all instances of the Newton classes have 1e-5
N1.tol # 1.e-5
N2.tol # 1e-4 but not N2.

类方法

我们在前面关于绑定和未绑定方法的章节中看到了方法是如何绑定到类的实例或者作为未绑定方法保持状态的。类方法不同。它们总是绑定方法。他们与类本身息息相关。

我们将首先描述语法细节，然后给出一些例子来说明这些方法可以用来做什么。要表明一个方法是一个类方法，装饰行在方法定义之前:

@classmethod

虽然标准方法通过使用它们的第一个参数来引用实例，但是类方法的第一个参数引用类本身。按照惯例，标准方法的第一个参数叫做self，类方法的第一个参数叫做cls。

• 标准案例:

      class A:
          def func(self,*args):
               <...>

• 类方法案例:

      class B:
          @classmethod
          def func(cls,*args):
               <...>

实际上，类方法对于在创建实例之前执行命令可能是有用的，例如在预处理步骤中。请参见以下示例:

在本例中，我们展示了如何在创建实例之前使用类方法准备数据:

class Polynomial:
    def __init__(self, coeff):
        self.coeff = array(coeff)
    @classmethod
    def by_points(cls, x, y):
        degree = x.shape[0] - 1
        coeff = polyfit(x, y, degree)
        return cls(coeff) 
    def __eq__(self, other):
        return allclose(self.coeff, other.coeff)

该类被设计为通过指定其系数来创建多项式对象。或者，by_points类方法允许我们通过插值点定义多项式。即使没有多项式实例可用，我们也可以将插值数据转换为多项式系数:

p1 = Polynomial.by_points(array([0., 1.]), array([0., 1.]))
p2 = Polynomial([1., 0.])

print(p1 == p2)  # prints True

本章后面的示例中介绍了类方法的另一个示例。在该示例中，类方法用于访问与该类的几个(或所有)实例相关的信息。

子类化和继承

在本节中，我们将介绍面向对象编程的一些核心概念:抽象类、子类和继承。为了指导您理解这些概念，我们考虑另一个数学示例:求解微分方程的一步方法。普通初值问题的一般形式是

数据为右侧函数 f 、初始值 x ₀ 、关注区间【 t ₀ 、t _e 】。这个问题的解决是一个功能。一种数值算法将该解作为离散值 u _i 的向量 u 给出，近似为x(t_I)。这里，是自变量 t 的离散值，在物理模型中常表示时间。

一步法通过递归步骤构造解值 u _i :

这里，φ是表征单个方法的阶跃函数(参见【28】):

• 显式欧拉 :
• 中点法则 :
• 龙格-库塔 4 : 带

我们在这里所做的是描述数学算法的典型方式。我们首先根据方法的思想来描述它，以抽象的方式给出它的步骤。为了实际使用它，我们必须填写一个具体方法的参数，在这个例子中，函数φ。这也是面向对象编程中解释事物的方式。首先，我们用方法的抽象描述建立一个类:

class OneStepMethod:
    def __init__(self, f, x0, interval, N):
        self.f = f
        self.x0 = x0
        self.interval = [t0, te] = interval
        self.grid = linspace(t0, te, N)
        self.h = (te - t0) / N

    def generate(self):
        ti, ui = self.grid[0], self.x0
        yield ti, ui
        for t in self.grid[1:]:
            ui = ui + self.h * self.step(self.f, ui, ti)
            ti = t
            yield ti, ui

    def solve(self):
        self.solution = array(list(self.generate()))

    def plot(self):
        plot(self.solution[:, 0], self.solution[:, 1])

    def step(self, f, u, t):
        raise NotImplementedError()

这个抽象类及其方法被用作单个方法的模板:

class ExplicitEuler(OneStepMethod):
    def step(self, f, u, t):
        return f(u, t)

class MidPointRule(OneStepMethod):
    def step(self, f, u, t):
        return f(u + self.h / 2 * f(u, t), t + self.h / 2)

请注意，在类定义中，我们用作模板的抽象类的名称OneStepMethod是作为额外的参数给出的:

class ExplicitEuler(OneStepMethod)

该类称为父类。父类的所有方法和属性都被子类继承，只要它们没有被覆盖。如果它们在子类中被重新定义，它们将被覆盖。step方法在子类中被重新定义，而方法generate对于整个族是通用的，因此从父族继承而来。在考虑进一步的细节之前，我们将演示如何使用这三个类:

def f(x, t):
    return -0.5 * x

euler = ExplicitEuler(f, 15., [0., 10.], 20)
euler.solve()
euler.plot()
hold(True)
midpoint = MidPointRule(f, 15., [0., 10.], 20)

midpoint.solve()
midpoint.plot()

使用星形运算符可以避免常见参数列表的重复(详见第七章、函数中变参数一节):

...
argument_list = [f, 15., [0., 10.], 20]
euler = ExplicitEuler(*argument_list)
...
midpoint = MidPointRule(*argument_list)
...

请注意，抽象类从未用于创建实例。由于step方法没有完全定义，调用它会引发类型为NotImplementedError的异常。

有时必须访问父类的方法或属性。这是使用命令super完成的。当子类使用自己的__init__方法来扩展父类的__init__时，这很有用:

例如，让我们假设我们想要给每个求解器类一个带有求解器名称的字符串变量。为此，我们为求解器提供了一个__init__方法，因为它覆盖了父级的__init__方法。在两种方法都应该使用的情况下，我们必须通过命令super来引用父方法:

class ExplicitEuler(OneStepMethod):
    def __init__(self,*args, **kwargs):
        self.name='Explicit Euler Method'
        super(ExplicitEuler, self).__init__(*args,**kwargs)
    def step(self, f, u, t):
        return f(u, t)

请注意，可以显式使用父类的名称。而super的使用允许我们更改父类的名称，而不必更改对父类的所有引用。

封装

有时使用继承是不切实际的，甚至是不可能的。这激发了封装的使用。我们将通过考虑 Python 函数来解释封装的概念，即 Python 类型function的对象，我们将其封装在一个新的类Function中，并提供一些相关的方法:

class Function:
    def __init__(self, f):
        self.f = f
    def __call__(self, x):
        return self.f(x)
    def __add__(self, g):
        def sum(x):
            return self(x) + g(x)
        return type(self)(sum) 
    def __mul__(self, g): 
        def prod(x):
            return self.f(x) * g(x)
        return type(self)(prod)
    def __radd__(self, g):
        return self + g
    def __rmul__(self, g):
        return self * g

请注意，__add__和__mul__操作应该返回同一个类的实例。这是通过return type(self)(sum)声明实现的，在这种情况下，这是一种更一般的写作形式return Function(sum)。我们现在可以通过继承来派生子类:

作为一个例子，切比雪夫多项式可以在区间[1，-1]中通过以下方式计算:

。

我们构建一个切比雪夫多项式作为Function类的实例:

T5 = Function(lambda x: cos(5 * arccos(x)))
T6 = Function(lambda x: cos(6 * arccos(x)))

切比雪夫多项式在某种意义上是正交的:

使用这种结构可以很容易地进行检查:

import scipy.integrate as sci

weight = Function(lambda x: 1 / sqrt((1 - x ** 2)))
[integral, errorestimate] = 
        sci.quad(weight * T5 * T6, -1, 1) # (6.510878470473995e-17, 1.3237018925525037e-14)

没有封装，像写weight * T5 * T6一样简单的乘法函数是不可能的。

装饰类

在第 7 章、功能的功能作为装饰者一节中，我们看到了如何通过应用另一个功能作为装饰者来修改功能。在前面的例子中，我们看到了只要类被提供了__call__方法，它们是如何表现为函数的。我们将在这里用它来展示类如何被用作装饰器。

让我们假设我们想要以这样一种方式改变一些函数的行为，即在调用函数之前，打印所有的输入参数。这可能对调试有用。我们以这种情况为例来解释装饰器类的使用:

class echo:
    text = 'Input parameters of {name}n'+
        'Positional parameters {args}n'+
        'Keyword parameters {kwargs}n'
    def __init__(self, f):
        self.f = f
    def __call__(self, *args, **kwargs):
        print(self.text.format(name = self.f.__name__,
              args = args, kwargs = kwargs))
        return self.f(*args, **kwargs)

我们使用这个类来修饰函数定义，

@echo
def line(m, b, x):
    return m * x + b

像往常一样调用函数，

line(2., 5., 3.)
line(2., 5., x=3.)

在第二次调用中，我们获得了以下输出:

Input parameters of line
Positional parameters (2.0, 5.0)
Keyword parameters {'x': 3.0}

11.0

这个例子表明类和函数都可以用作装饰器。类允许更多的可能性，因为它们也可以用来收集数据。

事实上，我们注意到:

• 每个修饰的函数都会创建装饰器类的一个新实例。
• 一个实例收集的数据可以通过类属性保存并供另一个实例访问(参见第 8 章、类中的属性一节)。

最后一点强调与功能装饰者的区别。我们现在通过一个装饰器来展示这一点，该装饰器对函数调用进行计数，并将结果存储在以函数为键的字典中。

为了分析算法的性能，统计特定函数的调用可能是有用的。我们可以在不改变函数定义的情况下获取计数器信息。该代码是对【4】中给出的一个例子的轻微修改。

class CountCalls:
    """
    Decorator that keeps track of the number of times 
    a function is called.
    """
    instances = {} 
    def __init__(self, f):
        self.f = f
        self.numcalls = 0
        self.instances[f] = self
    def __call__(self, *args, **kwargs):
        self.numcalls += 1
        return self.f(*args, **kwargs)
    @classmethod
    def counts(cls):
        """
        Return a dict of {function: # of calls} for all 
        registered functions.
        """
        return dict([(f.__name__, cls.instances[f].numcalls) 
                                    for f in cls.instances])

这里，我们使用类属性CountCalls.instances来存储每个单独实例的计数器。让我们看看这个装饰器是如何工作的:

@CountCalls
def line(m, b, x):
    return m * x + b
@CountCalls 
def parabola(a, b, c, x):
    return a * x ** 2 + b * x + c
line(3., -1., 1.)
parabola(4., 5., -1., 2.)

CountCalls.counts() # returns {'line': 1, 'parabola': 1}
parabola.numcalls # returns 1

总结

现代计算机科学中最重要的编程概念之一是面向对象编程。在本章中，我们学习了如何将对象定义为类的实例，我们提供了方法和属性。方法的第一个参数，通常用self表示，起着重要而特殊的作用。您看到了可用于为自己的类定义基本操作的方法，如+和*。

虽然在其他编程语言中，属性和方法可以防止意外使用，但 Python 允许通过特殊的 getter 和 setter 方法隐藏属性和访问这些隐藏属性的技术。为此，你遇到了一个重要的功能，property。

练习

Ex。1 →给类RationalNumber写一个方法simplify。这个方法应该以元组的形式返回分数的简化版本。

Ex。2 →为了给结果提供置信区间，在数值数学中引入了一种特殊的微积分，即所谓的区间算法；(参考【3、14】)。定义一个名为Interval的类，为其提供加减乘除的方法(仅限正整数)。这些操作遵循以下规则:

。

为该类提供允许类型为 a + I、a I、I + a、I a 的操作的方法，其中 I 是一个区间， a 是一个整数或浮点数。首先将整数或浮点数转换为区间[a,a]。(提示:您可能希望为此使用函数装饰器；(参见第七章、功能中作为装饰者的功能一节)。此外，实现__contains__方法，该方法使您能够使用“区间”类型的对象I的语法x in I来检查给定的数字是否属于区间。通过对一个区间应用多项式f=lambda x: 25*x**2-4*x+1来测试你的类。

Ex。3 →考虑类作为装饰者一节下的例子。扩展这个例子，得到一个函数装饰器，计算某个函数被调用的频率。

Ex。4 →比较两种实现类 RationalNumber中反向加法__radd__的方法:特殊方法一节例子中给出的方法和这里给出的方法:

class RationalNumber:
    ....
    def __radd__(self, other):
        return other + self

你认为这个版本会有错误吗？错误是什么，你如何解释？通过执行以下命令来测试您的答案:

q = RationalNumber(10, 15)
5 + q

Ex。4 →将装饰者类CountCalls视为装饰者一节中的示例。为这个类提供一个方法reset，将字典中所有函数CountCalls.instances的计数器设置为零。如果字典换成空字典会怎么样？