Python 科学计算第二版(二)

原文:annas-archive.org/md5/58a42b6d23877a7910fec539d6659c22

译者:飞龙

协议:CC BY-NC-SA 4.0

高级数组概念

在本章中,我们将解释一些数组的高级概念。首先,我们将介绍数组视图的概念——这是每个 NumPy 程序员必须了解的,以避免难以调试的编程错误。然后将介绍布尔数组及比较数组的方法。此外,我们还将简要描述索引和向量化,解释一些特殊话题,比如广播和稀疏矩阵。

在本章中,我们将覆盖以下主题:

  • 数组视图与副本

  • 比较数组

  • 数组索引

  • 性能与向量化

  • 广播

  • 稀疏矩阵

第六章:5.1 数组视图与副本

为了精确控制内存的使用,NumPy 提供了数组视图的概念。视图是共享同一数据的大数组的较小数组。这就像是一个对单一对象的引用。

5.1.1 数组视图

最简单的视图示例是通过数组切片获得的:

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

前面的切片 vM 的视图。它与 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

如果一个数组拥有其数据,那么属性 base 的值为 None

M.base # None

5.1.2 切片作为视图

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

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

5.1.3 通过转置和重塑生成视图

其他一些重要操作也会返回视图。例如,转置一个数组返回的是一个视图:

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

5.1.4 数组副本

有时需要显式地请求复制数据。这可以通过 NumPy 的 array 函数轻松实现:

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

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

N.base is None # True

在本节中,你了解了数组视图的概念。NumPy 通过使用视图而不是数组副本来节省内存,尤其对于大型数组来说,这一点至关重要。另一方面,无意中使用视图可能会导致难以调试的编程错误。

5.2 比较数组

比较两个数组并不像看起来那么简单。考虑以下代码,它旨在检查两个矩阵是否相近:

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()

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

5.2.1 布尔数组

布尔数组对于高级数组索引非常有用(另见第 5.3.1 节:使用布尔数组索引)。布尔数组就是条目类型为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语句。解决方法是使用allany方法来创建一个简单的TrueFalse

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")

在这里,使用allany方法之一将导致一个“标量”布尔值,这现在允许在if语句中进行数组比较。

5.2.2 数组相等性检查

检查两个浮动数组的相等性并不简单,因为两个浮点数可能非常接近,但不完全相等。在 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

5.2.3 数组上的布尔操作

不能在布尔数组上使用andornot。实际上,这些运算符会强制将数组转换为布尔值,这是不允许的。相反,我们可以使用表 5.1 中给出的运算符来对布尔数组执行逐元素的逻辑操作:

逻辑运算符布尔数组的替代方法
A and BA & B
A or BA &#124; 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
data = data + deviation 
# do not forget the parentheses in the 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

在这个示例中,我们首先创建了一个数据向量,然后用一些从正态分布随机数中抽样的偏差对其进行扰动。我们展示了两种替代方法,用于找到具有大扰动绝对值的数据元素,最后,我们收集了仅对小扰动的小数据值。在这里,我们在处理数组data时使用了布尔数组,而不是索引。这一技术将在下一节中进行解释。

5.3 数组索引

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

5.3.1 使用布尔数组索引

根据数组的值,通常有必要仅访问和修改数组的部分内容。例如,你可能想要访问数组中所有的正元素。实际上,这是通过布尔数组实现的,布尔数组像掩码一样,仅选择数组中的部分元素。这样的索引结果总是一个向量。例如,考虑以下示例:

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

事实上,命令M[B]等同于M[B].flatten()。你可以用另一个向量替换结果向量。例如,你可以用零替换所有元素:

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被巧妙地转换为适当形状的向量(另见第 5.5 节:广播)。

5.3.2 使用where命令

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

基本结构是:

where(condition, a, b)

当条件为True时,这将返回来自a的值,而当条件为False时,返回来自b的值。

例如,考虑一个海维赛德函数:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/460877e6-45a6-42bb-b9a6-c002278d76c7.png

以及它与where命令的实现:

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, -x, x))
# [ 4., 2., 0, 2., 4.] ]
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]))

这个示例演示了如何找出布尔数组中值为True的元素的索引。命令where是一个非常实用的工具,用于在数组中查找满足给定条件的元素。

在本节中,你看到了布尔数组的各种使用场景。每当你的代码中包含基于条件和数组操作的for循环时,检查是否可以使用布尔数组的概念来帮助去除不必要的for循环,并至少提高代码的可读性。

5.4 性能与向量化

当谈到 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

测量函数的速度是科学计算中的一个重要方面。有关测量执行时间的详细信息,请参见第 15.3 节:测量执行时间

5.4.1 向量化

为了提高性能,通常需要将代码向量化。用 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循环和其他较慢的构造外,还有一个有用的函数叫做vectorized(参见第 4.8 节:作用于数组的函数)。它接收一个函数,并创建一个向量化版本,该版本使用函数对数组中的所有元素进行操作,尽可能使用函数。

考虑以下函数,我们将用它来演示如何将函数向量化:

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

以非向量化的方式在一个包含 100 个元素的向量上应用此函数!

[my_func(vk) for vk in v]

使用向量化方式时,它的速度几乎是传统方式的三倍:

vectorize(my_func)(v)

在本节中,我们展示了几个使用 NumPy 数组进行计算向量化的例子。强烈推荐积极使用这一概念,这不仅能加速代码的执行,还能提高代码的可读性。

5.5 广播

在 NumPy 中,广播指的是能够推测两个数组之间的共同兼容形状。例如,当向量(单维数组)和标量(零维数组)相加时,标量会被扩展为一个向量,以便进行加法操作。这个一般机制叫做广播。我们首先从数学角度回顾这一机制,然后给出 NumPy 中广播的精确规则。数学视角可能让受过数学训练的读者更容易理解广播,而其他读者可能想跳过数学细节,直接继续阅读 第 5.5.2 节: 广播数组

5.5.1 数学视角

广播在数学中经常进行,通常是隐式的。例如,表达式如 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/69891ec8-0b8d-4396-9e75-2192d35dd16f.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/92db4562-4d42-4820-a1e6-0d22d11a35e1.png 就是广播的应用。我们将在这一节中详细描述该技术。

我们牢记函数与 NumPy 数组之间的密切关系,如 第 4.2.1 节所描述: 数组作为函数

常数函数

广播的一个最常见的例子是函数与常数相加;如果 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/66a23cd3-5e96-4255-9e41-2491f350b5da.png 是一个标量,我们通常写作:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/ef17988a-4b0c-4195-8d65-d04ad4f6c07d.png

这是一种符号滥用,因为你不应该能将函数和常数相加。然而,常数会隐式地广播到函数。常数 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/371dcbce-b9e2-4331-8072-7038d84c16ef.png 的广播版本是函数 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/6c1b9471-8cc7-4600-bff9-b0720043b17e.png,由以下公式定义:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/f15a1f56-1c7a-429c-a38b-e1481f7436e1.png

现在,将两个函数相加是有意义的:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/dfba5757-8042-4e11-ad24-6e7b231cc428.png

我们并不是为了矫揉造作,而是因为数组也可能出现类似的情况,如下方的代码所示:

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

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

这个例子极为简单,因此我们将展示一些不那么明显的情况。

多变量函数

当构建多个变量的函数时,会出现一个更为复杂的广播示例。例如,假设我们有两个单变量函数,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/ac094d60-643a-4288-aaa0-a59873bf1149.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/fb48d608-1c7e-47ed-92f8-b8348723c3a9.png,我们希望根据以下公式构造一个新的函数, https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/568707f0-40b8-4f1a-9483-66297d555257.png

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/9ee8e56a-18e9-423a-b540-587e6a6997a9.png

这显然是一个有效的数学定义。我们希望将这个定义表示为两个变量的函数之和,定义为:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/a348796b-1a2e-4cfe-9fe5-b78373d683e4.png

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

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/fdb81504-466e-467b-abc2-3c8b78df504c.png

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

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

这在采样两个变量的函数时尤其有用,正如 5.5.3 节所示:典型示例

通用机制

我们已经看到了如何将函数与标量相加,以及如何从两个单变量函数构建一个双变量函数。现在,让我们聚焦于使这一切成为可能的通用机制。这个通用机制包括两个步骤:重塑扩展

首先,函数https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/b3f744ae-3eb2-4a0d-ad22-e109b90e4b30.png重塑为函数https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/744ec18e-d9c8-48ae-90bc-52778afdef85.png,它接受两个参数。其中一个参数是虚拟参数,我们约定将其视为零:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/ba7c8aec-ae56-4b15-877a-7936cbb1b6f6.png

从数学角度看,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/c1a4c9a1-5a05-4bb5-936f-fd31524e29dc.png的定义域现在是https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/ec7c650c-102d-4c3a-820b-99effb8c4af3.png。然后,函数https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/d442142f-b300-4b26-9f85-a94741d9e562.png以类似于下述方式重塑

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/af5cef49-26d4-4daf-b64c-119ed180760a.png

现在,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/e3bf73d2-a1a3-44fb-af5e-c0779496df1f.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/13afea47-b115-429d-b132-28d5905a513f.png都接受两个参数,尽管其中一个始终是零。接下来我们进行扩展步骤。这与将常数转化为常数函数的步骤相同。

函数https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/5d7e6e3f-a7df-408b-9161-38c50b311745.png扩展为:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/4bf12815-d1ac-4df8-8ec5-b93901699031.png

函数https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/0b274b2e-637a-4aaa-b3ce-7e34acb94f31.png被扩展为:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/de132c36-12e2-42e2-9f1d-8cc0137eb711.png

现在,两个变量的函数,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/a8620c7e-d435-4557-85d9-be15cbafb7d3.png,它由https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/eed2403f-0afd-4741-bfc4-11fd7ab13a03.png粗略定义,可以在不参考其参数的情况下定义:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/8de1d5f3-13e4-4d33-afff-95477ad21e9d.png

例如,让我们描述之前的机制在常数情况下的表现。常数是标量,也就是零参数的函数。因此,重塑步骤是定义一个(一空)变量的函数:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/15e63522-e49d-4dd7-8114-35a33fab1d8c.png

现在,扩展步骤简单地进行如下:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/2e4751c9-6b14-4d2e-a08b-95f9cb52e64e.png

约定

最后一个要素是关于如何将额外参数添加到函数中的约定,也就是如何自动执行重塑。根据约定,函数会通过在左边添加零来自动重塑。

例如,如果一个有两个参数的函数https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/dbd64242-e716-4805-b924-487009f23a49.png需要重塑为三个参数,则新函数将由以下方式定义:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/35e31143-4398-41a3-be36-f193b706925d.png

在看过广播的更数学化动机之后,我们现在展示它如何应用于 NumPy 数组。

5.5.2 广播数组

我们现在将重复观察,数组只是几个变量的函数(见第 4.2 节:数学基础)。因此,数组广播完全遵循上述为数学函数所解释的相同程序。广播在 NumPy 中是自动完成的。

在图 5.1 中,我们展示了将形状为(4, 3)的矩阵加到形状为(1, 3)的矩阵时发生了什么。结果矩阵的形状是(4, 3):

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/c2e7cf3e-17d3-4938-a31f-0456a76e19ec.png

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

广播问题

当 NumPy 给定两个形状不同的数组,并要求执行需要这两个形状相同的操作时,两个数组会广播到一个共同的形状。

假设两个数组的形状分别为!和!。广播由以下两步组成:

  1. 如果形状!比形状!短,即len(s1) < len(s2),那么在形状!的左侧会添加 1。这就是重塑。

  2. 当形状长度相同时,第一个数组会扩展以匹配形状s[2](如果可能的话)。

假设我们想将一个形状为!的向量加到一个形状为!的矩阵中。该向量需要广播。第一个操作是重塑;将向量的形状从(3,)转换为(1, 3)。第二个操作是扩展;将形状从(1, 3)转换为(4, 3)。

例如,假设一个大小为n的向量!需要广播到形状(m, n):

  1. https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/d012613e-decb-42fa-aba9-5d3fe7dc155b.png会自动重塑为(1, n)。

  2. https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/077eabe8-be98-4d56-b46c-6424aa64c9e1.png被扩展到(m, n)。

为了演示这一点,我们考虑由以下定义的矩阵:

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

和一个向量给定:

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

现在我们可以直接将Mv相加:

M + v # works directly

结果是这个矩阵:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/9f0e0965-3ef7-48ff-9dbf-25c7ab17f013.png

形状不匹配

不能自动将长度为n的向量v广播到形状(n, m)。这一点在图 5.2中有所说明。

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/6c77c88d-e3f2-4a51-a36b-d8ce883f6059.png

图 5.2:由于形状不匹配导致广播失败

广播失败,因为形状(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,请参见第 4.6.3 节:重塑。)结果是这个矩阵:

array([[111, 112, 113, 114],  
       [221, 222, 223, 224], 
       [331, 332, 333, 334]])

5.5.3 典型示例

让我们看看一些典型示例,在这些示例中,广播可能会派上用场。

重标定行

假设M是一个*https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/8c502202-49f7-4a8e-a850-142921b05540.png矩阵,我们希望将每行乘以一个系数。系数存储在长度为n*的向量coeff中。在这种情况下,自动重塑将不起作用,我们必须执行:

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

重标定列

这里的设置相同,但我们希望使用存储在长度为m的向量coeff中的系数来重新标定每一列。在这种情况下,自动重塑将起作用:

rescaled = M*coeff

显然,我们也可以手动进行重塑,并通过以下方式达到相同的结果:

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

两个变量的函数

假设!和!是向量,我们希望形成矩阵!,其中的元素是!。这对应于函数!。矩阵!仅由以下方式定义:

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

如果向量*https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/7c329343-7ce8-40d9-b58e-ee57831cc5b8.png*和!分别是!和!,那么结果是:

array([[2, 3, 4],
       [3, 4, 5]])

更一般地,假设我们希望采样函数!。假设向量!和!已经定义,采样值的矩阵!可以通过以下方式获得:

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

请注意,这通常与ogrid结合使用。从ogrid获得的向量已经便于广播。这允许以下优雅的函数采样!

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.1.5 节:特殊方法)。因此,它是用方括号而不是圆括号来使用的。

这两个命令是等效的:

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]

5.6 稀疏矩阵

非零元素较少的矩阵称为sparse matrices(稀疏矩阵)。稀疏矩阵通常出现在科学计算中,例如在数值求解偏微分方程时,用于描述离散的微分算子。

稀疏矩阵通常具有较大的维度,有时大到整个矩阵(包含零元素)甚至无法放入可用内存。这是为稀疏矩阵设计特殊数据类型的一个动机。另一个动机是能够避免零矩阵元素,从而提高操作的性能。

对于一般的非结构化稀疏矩阵,在线性代数中只有极少数的算法可用。大多数算法是迭代型的,基于稀疏矩阵的矩阵-向量乘法的高效实现。

稀疏矩阵的例子包括对角矩阵或带状矩阵。这些矩阵的简单模式使得存储策略变得直接;主对角线以及上下对角线存储在一维数组中。从稀疏表示转换为经典数组类型以及反向转换,可以通过 NumPy 命令diag完成。

一般来说,稀疏矩阵并没有如此简单的结构,因此描述稀疏矩阵需要特殊的技术和标准。

这种矩阵的例子如图 5.3 所示。在该图中,像素表示在 1250 × 1250 矩阵中的非零元素。

在本节中,我们介绍了两种稀疏矩阵的行和列导向类型,这两者都可以通过模块scipy.sparse使用。

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/f80bc23d-f23a-4f62-bd62-7fa64296e1ce.png

图 5.3:来自弹性板有限元模型的刚度矩阵。

5.6.1 稀疏矩阵格式

模块scipy.sparse提供了多种不同的稀疏矩阵存储格式。在这里,我们只描述最重要的几种:CSR、CSC 和 LIL。LIL 格式应当用于生成和修改稀疏矩阵;CSR 和 CSC 是矩阵-矩阵和矩阵-向量运算的高效格式。

压缩稀疏行格式(CSR)

压缩稀疏行CSR)格式使用三个数组:dataindptrindices

  • 一维数组data按顺序存储所有非零值。它的元素数量等于非零元素的数量,通常用变量nnz表示。

  • 一维数组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

让我们来看一个例子:

矩阵的 CSR 格式:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3637523e-fddf-44da-8c67-6286fe700b5a.png

由以下三个数组给出:

data = array([1., 2., 3., 1., 4.])
indptr = array([0, 2, 2, 3, 5])
indices = array([0, 2, 0, 0, 3]

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

  • 使用二维数组作为参数

  • 使用 scipy.sparse 中的其他稀疏格式之一的矩阵

  • 使用形状参数 (m,n) 来生成 CSR 格式的零矩阵

  • 通过一个 1D 数组 data 和一个形状为 (2,len(data)) 的整数数组 ij,其中 ij[0,k] 是行索引,ij[1,k] 是矩阵中 data[k] 的列索引

  • 三个参数 dataindptrindices 可以直接传递给构造函数

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

考虑上面的示例;在 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

压缩稀疏列格式(CSC)

CSR 格式有一个列导向的对应物——压缩稀疏列CSC)格式。与 CSR 格式相比,它唯一的不同是 indptrindices 数组的定义,现在它们与列相关。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

基于行的链表格式(LIL)

链表稀疏格式按行存储非零矩阵项,存储在一个列表 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

这些操作在其他稀疏矩阵格式中不建议使用,因为它们非常低效。

5.6.2 生成稀疏矩阵

NumPy 命令 eyeidentitydiagrand 都有它们的稀疏版本。它们需要一个额外的参数,指定结果矩阵的稀疏矩阵格式。

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

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

在研究了不同的稀疏矩阵格式之后,我们现在转向稀疏矩阵的特殊方法,主要是转换方法。

5.6.3 稀疏矩阵方法

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

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

可以通过方法issparseisspmatrix_lilisspmatrix_csrisspmatrix_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或一个 1D 的 NumPyarray

考虑以下函数,我们将在此展示如何向量化一个函数:

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.sparsedot命令。

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

5.7 小结

视图的概念是本章你应该学习的重要主题之一。缺少这一主题会使你在调试代码时遇到困难。布尔数组在本书的多个地方出现。它们是避免冗长的if语句和循环处理数组时的简洁工具。在几乎所有大型计算项目中,稀疏矩阵都会成为一个问题。你已经看到如何处理它们以及可用的相关方法。

绘图

在 Python 中,绘图可以通过 matplotlib 模块的 pyplot 部分来完成。使用 matplotlib,你可以创建高质量的图形和图表,并可视化你的结果。Matplotlib 是开源的、免费的软件。Matplotlib 网站还包含了优秀的文档和示例,详见 35。在本节中,我们将展示如何使用最常见的功能。接下来的示例假设你已经导入了该模块:

from matplotlib.pyplot import *

如果你希望在 IPython 中使用绘图命令,建议在启动 IPython shell 后立即运行 magic command %matplotlib。这将为 IPython 准备好交互式绘图功能。

第七章:6.1 使用基本绘图命令绘制图形

本节中,我们将通过基本命令创建图形。这是学习如何使用 Python 绘制数学对象和数据图形的入门点。

6.1.1 使用 plot 命令及其一些变体

标准的绘图函数是 plot。调用 plot(x,y) 会创建一个图形窗口,并绘制出 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/0bebbcfd-cc30-44ee-a0a3-5e07f610d63b.png 作为 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/17e2c94a-5fd0-415b-a4e3-217aebb2fa23.png 的函数图像。输入参数是等长的数组(或列表)。也可以使用 plot(y),在这种情况下,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/77346714-da93-4af9-97f6-a5ccf7e135d1.png 中的值将会根据其索引绘制,也就是说,plot(y)plot(range(len(y)),y) 的简写。

下面是一个示例,展示如何使用 200 个样本点和每隔四个点设置标记来绘制 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/d15e859d-6238-4437-987e-6f989435ab76.png

# 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):

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/f8bca15b-e7cd-42d6-bc6b-1dc72d497cd7.png

图 6.1:绘制函数 sin(x) 的图像,并显示网格线

如你所见,标准的绘图是一条实心的蓝色曲线。每个坐标轴都会自动缩放以适应数值,但也可以手动设置。颜色和绘图选项可以在前两个输入参数后提供。在这里,r* 表示红色星形标记。格式设置将在下一节中详细讨论。命令 title 在绘图区域上方添加标题文本字符串。

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

可以通过使用 legend 函数并为每个绘图调用添加标签来解释多个图形。以下示例使用 polyfitpolyval 命令拟合多项式并绘制结果,同时添加图例:

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

# plot the two 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):

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/a14516df-6320-4707-9dde-d9e1f12eb6ce.png

图 6.2:两个多项式拟合同一点集

作为基础绘图的最终示例,我们演示了如何绘制散点图和二维对数图。

这是一个二维点散点图的示例:

# 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')

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/1abff8af-22ad-4ced-b8a8-6d55a08d1841.png

图 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')

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/79c6a098-ec45-455b-8a16-454ae29cd3b2.png

图 6.3(b):一个带有对数 xy 坐标轴的图形示例

前面图中展示的示例(图 6.3(a)图 6.3(b)) 使用了plotloglog的一些参数,这些参数允许特殊的格式化。在下一节中,我们将更详细地解释这些参数。

6.1.2 格式设置

图形和绘图的外观可以通过样式和定制设置为你想要的效果。一些重要的变量包括linewidth(控制绘图线条的粗细),xlabelylabel(设置坐标轴标签),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='-')

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/eb6fd063-8472-4188-8a3e-a95226d58bad.png

表 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,density=True)
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(fr'Histogram with $\mu={mu}, \sigma={sigma}')

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/9d2b730a-47c1-40ff-b0cd-b5ff99327b78.png

图 6.4:具有 50 个区间的正态分布图,绿色曲线表示真实分布

结果图形看起来与图 6.4类似。标题以及其他任何文本都可以使用 LaTeX 格式化来显示数学公式。LaTeX 格式化被包含在一对$符号之间。另请注意,使用format方法进行的字符串格式化;参见第 2.4.3 节,字符串格式化

有时,字符串格式化的括号会与 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*numpy.random.standard_normal(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)

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/4ece0cd7-35f1-4ea1-80a4-606d8f8b4633.png

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

函数avg使用 NumPy 的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')

6.1.3 使用 meshgrid 和等高线

一个常见的任务是对矩形区域内的标量函数进行图形化表示:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/bf7a60e2-808a-4314-8834-89fbc8f13bed.png

为此,我们首先需要在矩形上生成一个网格!。这是通过命令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))

XY是形状为(n,m)的数组,其中X[i,j]Y[i,j]包含网格点的坐标![],如图 6.6所示:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/f5b1481f-1499-4c11-8758-b8049dd174cb.png

图 6.6:一个由 meshgrid 离散化的矩形。

一个由meshgrid离散化的矩形将在下一节中用于可视化迭代过程,而我们将在这里用它绘制函数的等高线。这是通过命令contour来完成的。

作为示例,我们选择了罗森布鲁克香蕉函数:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3ea8cf3c-24a2-416a-a130-2cefe007da25.png

它用于挑战优化方法,见[27]。函数值向着一个香蕉形的山谷下降,该山谷本身慢慢下降,最终到达函数的全局最小值(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

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/297117a6-be35-49fb-bd9a-e882ec0a0d08.png

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

在前面的示例中,使用了lambda关键字表示的匿名函数来保持代码简洁。匿名函数的解释见第 7.7 节,匿名函数。如果未将级别作为参数传递给contour,该函数会自动选择合适的级别。

contourf函数执行与contour相同的任务,但根据不同的级别填充颜色。等高线图非常适合可视化数值方法的行为。我们通过展示优化方法的迭代过程来说明这一点。

我们继续前面的示例,并描绘了通过 Powell 方法生成的罗森布鲁克函数最小值的步骤,我们将应用该方法来找到罗森布鲁克函数的最小值:

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应用 Powell 方法找到最小值。通过给定的起始值*https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/899aeac0-4ed8-4a07-aa68-aefa0abcfbf9.png。迭代过程在等高线图中以子弹点表示;见图 6.8

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/b8264e8e-e652-40ce-b542-cd5933ccb2b3.png

图 6.8:罗森布鲁克函数的等高线图,展示了优化方法的搜索路径

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

6.1.4 生成图像和轮廓

让我们看一些将数组可视化为图像的示例。以下函数将为曼德尔布罗特分形创建一个颜色值矩阵,另见[20]。这里,我们考虑一个依赖于复数参数的固定点迭代,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/7da969b2-9b13-46c2-8302-ac87dcd6ea89.png

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/94a571d8-b85c-4c81-b346-4715ff922487.png

根据此参数的选择,它可能会或可能不会创建一个有界的复数值序列,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/79494086-b942-40d5-8ddd-9d5909102f58.png

对于每个https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/8f79f3af-5180-4652-be25-3ac505e9c0e1.png的值,我们检查![]是否超过了预定的界限。如果在maxit次迭代内仍然低于该界限,则认为序列是有界的。

请注意,在以下代码片段中,meshgrid用于生成一个复数参数值矩阵,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/23c90700-fa16-49ff-a92c-4d06ad22376d.png

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') 来关闭坐标轴,因为这对于图像来说可能不太有用。

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/509af838-00f1-4ec8-acce-bf44d789faa4.png

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

默认情况下,imshow 使用插值来使图像看起来更漂亮。当矩阵较小时,这一点尤为明显。下图显示了使用以下方法的区别:

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

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

在第二个示例中,像素值只是被复制了,见 [30]

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/ae936a8e-16d1-4fd1-8a10-bc1a5c07efba.png

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

有关使用 Python 处理和绘制图像的更多详细信息。

在了解了如何以“命令方式”制作图表之后,我们将在接下来的部分中考虑一种更面向对象的方法。虽然稍微复杂一些,但它打开了广泛的应用范围。

6.2 直接使用 Matplotlib 对象

到目前为止,我们一直在使用 matplotlib 的 pyplot 模块。这个模块使我们可以直接使用最重要的绘图命令。通常,我们感兴趣的是创建一个图形并立即显示它。然而,有时我们希望生成一个图形,稍后可以通过更改某些属性来修改它。这要求我们以面向对象的方式与图形对象进行交互。在这一节中,我们将介绍修改图形的一些基本步骤。要了解 Python 中更复杂的面向对象绘图方法,您需要离开 pyplot,直接进入 matplotlib,并参考其广泛的文档。

6.2.1 创建坐标轴对象

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

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.7 节 中解释这个构造,匿名函数。实际上,这两个绘图命令将列表 ax.lines 填充了两个 Lines2D 对象:

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(左图) 所示。

6.2.2 修改线条属性

我们刚刚通过标签标识了一个特定的线条对象。它是列表 ax.lines 中的一个元素,索引为 il。它的所有属性都被收集在一个字典中:

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,(右)*所示:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/fb85293c-d3b5-49dc-b8b7-c982db0cf062.png

图 6.11:幅度调制正弦函数(左)和数据点被破坏的曲线(右)

6.2.3 制作注释

一个有用的坐标轴方法是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让箭头呈现为样条曲线或其他形状。

文本属性,甚至是围绕文本的边框框,可以通过额外的关键字参数传递给annotate方法,见图 6.12,(左)

尝试注释有时需要多次尝试,我们需要丢弃其中的一些。因此,我们将注释对象赋值给一个变量,这样就可以通过其remove方法移除注释:

annot1.remove()

6.2.4 填充曲线之间的区域

填充是一个理想的工具,用于突出曲线之间的差异,例如预期数据上的噪声和近似函数与精确函数之间的差异。

填充是通过坐标轴方法fill_between完成的:

ax.fill_between(x,y1,y2)

对于下一个图,我们使用了以下命令:

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

在上一章中,我们已经了解了 NumPy 方法where。在这里的上下文中,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给出。

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

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/67407da9-48ee-4d76-9731-3596707d1634.png

图 6.12:带有注释和填充区域的幅度调制正弦函数(左),以及仅通过使用where参数部分填充区域的修改图(右)

如果你自己测试这些命令,记得在尝试部分填充之前移除完整填充,否则你将看不到任何变化:

axf.remove()

相关的填充命令是fill用于填充多边形,fill_betweenx用于填充水平方向的区域。

6.2.5 定义刻度和刻度标签

在演讲、海报和出版物中的图形,如果没有过多不必要的信息,看起来会更美观。你希望引导观众关注那些包含信息的部分。在我们的例子中,我们通过去除x轴和y轴的刻度,并引入与问题相关的刻度标签来清理图像:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/cc6a84f4-c69e-4d7d-ad83-e4bf7c748600.png

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

图 6.13中的刻度线是通过以下命令设置的。注意使用 LaTeX 方式设置带有希腊字母的标签:

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所示。

6.2.6 设置脊柱使你的图更具启发性——一个综合示例

脊柱是显示坐标的带有刻度和标签的线条。如果不采取特别措施,Matplotlib 会将它们放置为四条线——底部、右侧、顶部和左侧,形成由坐标轴参数定义的框架。

通常,图像在没有框架时看起来更好,而且脊柱有时可以放置在更具教学意义的位置。在这一节中,我们展示了改变脊柱位置的不同方法。

让我们从一个指导示例开始,见 6.14。

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/edcc56e5-af01-4b23-918a-f5ea1413c186.png

图 6.14:一个 Matplotlib 图,具有非自动的脊柱位置设置

在这个例子中,我们选择只显示四个脊柱中的两个。

我们通过使用set_visible方法取消选择了顶部和右侧的脊柱,并通过使用set_position方法将左侧和底部的脊柱放置在数据坐标中:


fig = figure(1)
ax = fig.add_axes((0.,0.,1,1))
ax.spines["left"].set_position(('data',0.))
ax.spines["bottom"].set_position(("data",0.))
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
x=linspace(-2*pi,2*pi,200)
ax.plot(x,arctan(10*x), label=r'$\arctan(10 x)$')
ax.legend()

脊柱携带刻度和刻度标签。通常,它们是自动设置的,但手动设置它们往往更有优势。

在以下示例中,我们甚至利用了两组刻度线的可能性,并且设置了不同的放置参数。Matplotlib 将这两组刻度线分别称为’minor’和’major’。其中一组用于水平对齐y轴左侧的刻度标签:

ax.set_xticks([-2*pi,-pi,pi,2*pi])
ax.set_xticklabels([r"$-2\pi$",r"$-\pi$",r"$\pi$", r"$2\pi$"])
ax.set_yticks([pi/4,pi/2], minor=True)
ax.set_yticklabels([r"$\pi/4$", r"$\pi/2$"], minor=True)
ax.set_yticks([-pi/4,-pi/2], minor=False,)
ax.set_yticklabels([r"$-\pi/4$", r"$-\pi/2$"], minor=False) # major label set
ax.tick_params(axis='both', which='major', labelsize=12)
ax.tick_params(axis='y', which='major',pad=-35) # move labels to the right
ax.tick_params(axis='both', which='minor', labelsize=12)

结果如图 6.15 所示。

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/dd735ee6-c02a-436f-b715-02b266371708.png

图 6.15:改变刻度和标签的位置

这个例子可以通过添加更多的轴和注释来进一步展开。我们参考了练习 7图 6.20

到目前为止,我们考虑的是二维图。接下来,我们将在下一节讨论三维数学对象的可视化。

6.3 制作三维图

有一些有用的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.16 中看到。

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/9827dbe5-eab1-4a5d-b0ba-1a75281e3a42.png

图 6.16:使用 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)

alpha 值设置透明度。表面图如图 6.17所示。

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/acf60e0f-5670-40bc-863f-5156dc8d1cb3.png

图 6.17:绘制表面网格的示例

你还可以在任意坐标投影中绘制等高线,如以下示例所示:

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])。该命令适用于二维图。然而,axis([-40,40,-40,40,-40,40]) 则无效。对于三维图,你需要使用面向对象的命令版本 ax.set_xlim3d(-40,40)。同样,设置坐标轴标签时也有类似的命令。对于二维图,你可以使用 xlabel('X axis')ylabel('Y axis'),但没有 zlabel 命令。对于三维图,你需要使用 ax.set_xlabel('X axis') 和类似的命令设置其他标签,如前面的示例所示。

这段代码生成的图形如下:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/90bbde60-4d49-4366-9efe-c6304645e46c.png

图 6.18:带有额外等高线图的三维图,在三个坐标投影中展示

有许多选项可供设置图形的外观,包括颜色和表面透明度。mplot3d 文档网站 [23] 中有详细信息。

数学对象有时通过一系列图片甚至电影来动态可视化效果更佳。这是下一节的主题。

6.4 从图形生成电影

如果你有演变的数据,可能希望将其保存为电影,并在图形窗口中显示,类似于命令 savefig。一种方法是使用模块 visvis,请参阅 [37]。

这里是一个使用隐式表示演变圆形的简单示例。令圆形由一个函数的零水平表示,

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/53d5e400-f103-4d73-8d95-1d9f771a7e14.png 一个函数 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3cad0783-6c67-4da2-8bed-6f3de9a998f8.png 的零水平。

另外,考虑到盘面 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/4f11df5a-0c9f-413a-9548-2dc4796ac29b.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/6f82cbed-2f3f-4614-9e23-eca35223238a.png 的零集内。如果 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/2b7a8922-2cce-4811-8ef9-6fce4d09ab58.png 的值以速率 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/66d4131e-1c59-423c-8242-933277ff6eba.png 递减,则圆圈将以速率 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/cc1d309a-3955-4e88-980a-1c4eff3b32a9.png 向外移动。

这可以实现为:

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)

结果是一个黑色圆圈逐渐扩大的动画(*.swf 文件),如 图 6.19 所示。

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/0683ec0e-d3db-4236-aac8-c127bc6ebfd6.png https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/af022cbe-d334-4cd1-8b0f-a3179199221e.png https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/cb4d0486-dedf-456f-a4ad-60e895a754cd.png https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/03a96259-8454-4516-85ef-f3a9a87f52ff.png

图 6.19:演变圆圈的示例

在这个示例中,使用了一组数组来创建动画。模块 visvis 也可以保存 GIF 动画,并且在某些平台上,可以生成 AVI 动画(*.gif*.avi 文件)。此外,还可以直接从图形窗口捕捉电影帧。然而,这些选项要求系统安装更多的包(例如,PyOpenGLPIL,即 Python 图像库)。有关更多细节,请参阅 visvis 官方网页上的文档。

另一种选择是使用 savefig 创建图像,为每一帧生成一张图像。以下代码示例创建了一系列 200 张图片文件,这些文件可以合并成一个视频:

# 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, aspect='auto')
    gray()
    axis('off')
    savefig('circle_evolution_{:d}.png'.format(iteration))
    f -= 1

这些图像可以使用标准视频编辑软件进行合并,例如 Mencoder 或 ImageMagick。该方法的优点是你可以通过保存高分辨率图像来制作高分辨率视频。

6.5 小结

图形表示是展示数学结果或算法行为最紧凑的形式。本章为你提供了绘图的基本工具,并介绍了一种更精细的面向对象的图形对象工作方式,例如图形、坐标轴和线条。

在本章中,你学习了如何绘制图形,不仅是经典的 x/y 图,还有 3D 图和直方图。我们还为你提供了制作影片的前菜。你还看到了如何修改图形,视其为图形对象,并使用相关的方法和属性进行设置、删除或修改。

6.6 练习

示例 1: 编写一个函数,给定椭圆的中心坐标 (x,y),半轴 ab 以及旋转角度,绘制椭圆 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/41ba9462-862e-42f9-9719-4723a30ee33d.png

示例 2: 编写一个简短的程序,接受一个二维数组,例如前面的曼德尔布罗特轮廓图,并迭代地将每个值替换为其邻居的平均值。在图形窗口中更新数组的轮廓图,以动画形式展示轮廓的演变。解释其行为。

示例 3: 考虑一个 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/ab1805a6-ee8e-44b4-a1e3-043c630fc03e.png 矩阵或整数值图像。映射

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/6fb1559c-a1fb-45f5-b477-e3ba88bd5d7e.png

是一个点阵网格映射到自身的示例。这个方法有一个有趣的特性,它通过剪切然后使用模函数mod将超出图像的部分移回图像内,进而使图像随机化,最终恢复到原始图像。依照以下顺序实施,

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/58f5266a-a958-4494-a86e-0fec0a7dca8d.png,

并将前几个步骤!保存为文件或绘制到图形窗口中。

作为示例图像,你可以使用经典的 512×512 Lena 测试图像,该图像来自scipy.misc

from scipy.misc import lena
I = lena()

结果应该如下所示:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/1c2ca42a-82c1-49a3-a149-d04714f1d85e.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/862a99ab-c967-4f3d-9632-cffe8776bc6c.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/e4cea046-8a57-4eba-a979-c88ea1971f44.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/122156fe-df35-41e5-a1e2-6b26cb575482.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3173fa79-21df-4333-8883-f524e70bc76a.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/f7dc2216-0045-488d-a339-020a20ed5fa6.png
01128256511512

计算xy映射,并使用数组索引(参见第 5.3 节:数组索引)来复制像素值。

Ex. 4: 读取并绘制图像。SciPy 提供了imread函数(位于scipy.misc模块中)来读取图像(参见第 14.6 节: 读取和写入图像)。编写一个简短的程序,从文件中读取图像,并在原始图像上叠加给定灰度值的图像轮廓。你可以通过像这样平均颜色通道来获得图像的灰度版本:mean(im,axis=2)

Ex. 5: 图像边缘。二维拉普拉斯算子的零交叉是图像边缘的一个很好指示。修改前一个练习中的程序,使用scipy.ndimage模块中的gaussian_laplacelaplace函数来计算二维拉普拉斯算子,并将边缘叠加到图像上。

Ex. 6: 通过使用orgid代替meshgrid,重新编写曼德博集合分形示例第 6.1.4 节:生成图像和轮廓。参见第 5.5.3 节对ogrid的解释,典型示例orgidmgridmeshgrid之间有什么区别?

Ex. 7:6.20 中,研究了使用反正切函数来近似跳跃函数(Heaviside 函数)。该曲线的一部分被放大以可视化近似的质量。通过你自己的代码重现这张图。

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/b9ee6bc2-3edd-4d2b-b244-2d86f5b611aa.png

图 6.20:使用反正切函数近似 Heaviside 函数(跳跃函数)

函数

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

本章将涉及以下主题:

  • 数学中的函数与 Python 中的函数

  • 参数和参数值

  • 返回值

  • 递归函数

  • 函数文档

  • 函数是对象

  • 匿名函数 – 关键字lambda

  • 函数作为装饰器

第八章:7.1 数学中的函数与 Python 中的函数

在数学中,函数表示为一个映射,它唯一地将域*!中的每个元素与范围中的对应元素*相联系。

这通过*!*来表示。

或者,在考虑特定元素*!时,可以写成*。

这里,*!被称为函数的名称,而是其应用于时的值。这里,有时被称为*的参数。在考虑 Python 中的函数之前,让我们先看一个示例。

例如,。这个函数将两个实数映射到它们的差值。

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

在这种情况下,返回的是一个实数。在处理函数时,我们需要区分两个不同的步骤:

  • 函数的定义

  • 函数的求值,也就是计算给定值的*!对于*

第一步只需要执行一次,而第二步可以针对不同的参数执行多次。

编程语言中的函数遵循类似的概念,并将其应用于各种类型的输入参数,例如字符串、列表、浮动数或任何对象。我们通过再次考虑给定的示例来演示函数的定义:

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

关键字def表示我们将定义一个函数。subtract是函数的名称,x1x2是它的参数。冒号表示我们正在使用一个代码块。函数返回的值跟在关键字return后面。

现在,我们可以评估这个函数。该函数在其参数被输入参数替代后被调用:

r = subtract(5.0, 4.3)

结果0.7被计算并赋值给变量r

7.2 参数和参数值

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

7.2.1 传递参数——通过位置和关键字

我们将再次考虑之前的例子,其中函数有两个参数,分别是x1x2

它们的名字用于区分这两个数,这两个数在此情况下不能互换,否则会改变结果。第一个参数定义了从中减去第二个参数的数字。当subtract函数被调用时,每个参数都被替换为一个参数值。参数的顺序很重要;参数可以是任何对象。例如,我们可以调用如下:

z = 3 
e = subtract(5,z)

除了这种标准的调用函数方式,即通过位置传递参数,有时使用关键字传递参数可能会更方便。参数的名称就是关键字;考虑以下示例:

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

在这里,参数是通过名称分配给参数的,而不是通过调用中的位置。两种调用函数的方式可以结合使用,使得位置参数排在前面,关键字参数排在后面。我们通过使用函数plot来演示这种方式,plot函数在第 6.1 节:基本绘图中有描述:

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

7.2.2 改变参数值

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

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]

这样的函数错误地使用其参数来返回结果。我们强烈劝阻使用这种构造,并建议你在函数内部不要更改输入参数(有关更多信息,请参见第 7.2.4 节:默认参数)。

7.2.3 访问定义在局部命名空间外的变量

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

全局变量在处理闭包时非常有用;请参见第 7.7 节中的相关示例:匿名函数—— 关键字 lambda*。

7.2.4 默认参数

有些函数可能有很多参数,其中一些参数可能只在非标准情况下才有意义。如果参数可以自动设置为标准(默认)值,那将是非常实用的。

我们通过查看模块scipy.linalg中的命令norm来演示默认参数的使用。它计算矩阵和向量的各种范数。更多关于矩阵范数的信息,请参见[10 ,§2.3]

以下用于计算 Frobenius 范数的调用是等效的:

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]

回想一下,列表是可变对象。

7.2.5 可变数量的参数

列表和字典可以用来定义或调用具有可变数量参数的函数。我们可以定义一个列表和一个字典,如下所示:

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

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

plot(*data,**style)

*开头的变量名,例如前面示例中的*data,意味着将一个列表解包以向函数提供其参数。通过这种方式,列表生成位置参数。类似地,带有**前缀的变量名,例如示例中的**style,将解包一个字典为关键字参数;见 图 7.1

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/af27c774-d30e-4912-82c4-83562b036b9d.png

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

你也可能想要使用反向过程,在这种情况下,所有给定的位置参数会被打包成一个列表,所有的关键字参数会被打包成一个字典并传递给函数。在函数定义中,这通过分别以***作为前缀的参数来表示。你经常会在代码文档中看到*args**kwargs这两个参数;参见图 7.2

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/351df170-0fda-4879-a49e-9bf80d2526e8.png

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

7.3 返回值

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 formedcite

(另请参见欧拉公式,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/de874ce0-e2fa-4001-9103-6dd2ec877c4a.png)。

在这里,我们使用了 NumPy 函数sqrt(x)来计算数字x的平方根,并使用arctan2(x, y)来表示https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/b3c0157c-8a54-45ae-99c2-05f5c7ef11d5.png

让我们试一下我们的函数:

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)

我们可以通过调用在练习部分的练习 1中定义的polar_to_comp函数来测试我们的函数。

如果一个函数没有return语句,它将返回值None。有很多情况,函数不需要返回任何值。这可能是因为传递给函数的变量可能会被修改。例如,考虑下面的函数:

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

前面的函数不返回任何内容,因为它修改了作为可变参数传递的一个对象。有很多方法的行为也是如此。仅列举列表方法,appendextendreversesort这些方法都不返回任何内容(即它们返回None)。当一个对象通过这种方式被方法修改时,称为就地修改。很难知道一个方法是否会改变一个对象,除非查看代码或文档。

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

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

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

7.4 递归函数

在数学中,许多函数是递归定义的。在本节中,我们将展示如何在编写函数时使用这个概念。这使得程序与其数学对应物之间的关系变得非常清晰,这可能有助于提高程序的可读性。

然而,我们建议谨慎使用这种编程技巧,尤其是在科学计算中。在大多数应用中,更直接的迭代方法通常更高效。通过以下示例,这一点将立刻变得清晰。

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

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/599c7304-96a1-4b1c-8d58-e91f7d52874b.png

这种递归需要初始化,即 ![]

在 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)

要计算 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/30b7eadf-af02-468d-84cb-081ad13a56cb.png,函数可以像这样调用:

chebyshev(5, 0.52) # returns 0.39616645119999994

这个示例还展示了浪费计算时间的巨大风险。随着递归层级的增加,函数评估的数量呈指数增长,而且这些评估中的大多数只是先前计算的重复结果。虽然使用递归程序来展示代码与数学定义之间的紧密关系可能很诱人,但生产代码通常会避免使用这种编程技巧(参见练习6 节中的练习部分)。我们还提到了一种叫做记忆化(memoization)的技巧,它将递归编程与缓存技术相结合,以保存重复的函数评估,详见[22]

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

  • 基本情况;这里是前两个 if 分支

  • 递归主体,在这个主体中,函数本身会使用较小级别的参数一次或多次被调用

执行递归函数时经过的层数称为递归深度。这个值不应过大,否则计算可能变得低效,并且在极端情况下,会抛出以下错误:

RuntimeError: maximum recursion depth exceeded

最大的递归深度取决于你使用的计算机的内存。这个错误也会在函数定义缺少初始化步骤时发生。我们鼓励仅在非常小的递归深度下使用递归程序(更多信息,请参见第 9.7.2 节:递归)。

7.5 函数文档

你应该在函数开始时使用一个字符串来记录文档,这个字符串叫做 文档字符串

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__。你应该在文档字符串中提供的最基本信息是函数的目的以及输入和输出对象的描述。有一些工具可以通过收集程序中的所有文档字符串来自动生成完整的代码文档(更多信息,请参见 Sphinx 的文档,[32])。

7.6 函数是对象

函数是对象,就像 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的信息:

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

7.6.1 部分应用

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

函数 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/9502da7b-86bf-48e3-8564-7821236cc621.png 可以看作是一个双变量的函数。通常,你会把 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/338903f2-95ad-4f1e-a84c-eafbfeb69218.png 视为一个固定的参数,而不是一个自由变量,属于一族函数 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/eeaff405-ef7d-405b-8f02-e3bff86d22fb.png

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/61ae8234-2d80-4209-af59-a3d132be42cb.png

这个解释将一个双变量函数简化为一个单变量函数 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3be52ca7-be0d-413c-990e-475b16061773.png,其中固定了一个参数值 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/49a1eac7-3ffe-48ef-9a3d-3f4b0cca152d.png。通过固定(冻结)函数的一个或多个参数来定义一个新函数的过程称为部分应用。

使用 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)

fomega=make_sine(0.25)
fomega(3) # returns -1.0

在最后一行,新的函数在 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/7dab2cdb-e183-49f4-a87d-9369f0a1240f.png 被求值。

7.6.2 使用闭包

从函数是对象的角度出发,可以通过编写一个函数来实现部分应用,这个函数本身返回一个新的函数,并且输入参数的数量减少。例如,函数make_sine可以定义如下:

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

在这个例子中,内嵌函数mysine可以访问变量freq;它既不是该函数的局部变量,也没有通过参数列表传递给它。Python 允许这样的构造,参见章节 13.1,命名空间

7.7 匿名函数 —— 关键字 lambda

关键字lambda在 Python 中用于定义匿名函数,也就是没有名字、由单一表达式描述的函数。你可能只想对一个可以通过简单表达式表示的函数执行某个操作,而不需要给这个函数命名,也不需要通过冗长的def块来定义它。

名字lambda源自微积分和数学逻辑的一个特殊分支,即 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/8acb12a5-c98d-4653-ac41-9e4fba93f10a.png-微积分。

我们通过数值评估以下积分来演示lambda函数的使用:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/f7b1ca59-f9c8-4ae1-965b-1fe67f6f9331.png

我们使用 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

7.7.1 lambda构造总是可以替换的

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

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

使用这种构造的主要原因是对于非常简单的函数来说,完整的函数定义会显得过于繁琐。

lambda函数提供了创建闭包的第三种方式,正如我们通过继续前面的例子![]所演示的那样。

我们使用第 7.6.1 节中的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.) )

7.8 函数作为装饰器

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

假设我们有一个函数用来确定矩阵的稀疏度:

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

如果未以数组对象作为输入调用此函数,则会返回一个错误。更准确地说,它将无法与没有实现reshape方法的对象一起工作。例如,how_sparse函数无法与列表一起工作,因为列表没有reshape方法。以下辅助函数修改任何具有一个输入参数的函数,以便尝试将其类型转换为数组:

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

因此,修改后的函数how_sparse = cast2array(how_sparse)可以应用于任何可以转换为数组的对象。如果how_sparse的定义用这个类型转换函数进行装饰,也可以实现相同的功能:

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

要定义一个装饰器,你需要一个可调用的对象,例如一个修改被装饰函数定义的函数。其主要目的包括:

  • 通过将不直接服务于函数功能的部分分离来增加代码可读性(例如,记忆化)

  • 将一组类似函数的公共前言和尾部部分放在一个共同的地方(例如,类型检查)

  • 为了能够方便地开关函数的附加功能(例如,测试打印或追踪)

还建议考虑使用functools.wraps,详情请见[8]

7.9 小结

函数不仅是使程序模块化的理想工具,还能反映出数学思维。你已经学习了函数定义的语法,并了解如何区分定义函数和调用函数。

我们将函数视为可以被其他函数修改的对象。在处理函数时,了解变量的作用域以及如何通过参数将信息传递到函数中是非常重要的。

有时,定义所谓的匿名函数非常方便。为此,我们引入了关键字lambda

7.10 练习

例 1:编写一个函数 polar_to_comp,该函数接收两个参数 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/e834ef4c-2926-4ad9-9504-de4b2234bc13.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/4f3160df-59d2-427a-b19a-51eb90212f65.png,并返回复数 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/6e1ade86-768c-4fe4-9963-7bde1348f589.png。使用 NumPy 函数 exp 来计算指数函数。

例 2:在 Python 模块 functools 的描述中[8],你会找到以下 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

解释并测试此函数。

例 3:为函数 how_sparse 编写一个装饰器,该装饰器通过将小于 1.e-16 的元素设置为零来清理输入矩阵 A(参考 第 7.8 节:作为装饰器的函数)。

例 4:一个连续函数 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/edafb426-3175-4e5d-af76-2d21dc0acd01.png,其 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3c9e124d-b89a-4263-ab14-04b3293dfaea.png 在区间 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/bff068a9-97e1-48a0-a3f5-3f16aa97fd84.png 内改变符号,并且在该区间内至少有一个根(零)。可以通过二分法找到此根。该方法从给定区间开始,接着检查子区间中的符号变化。

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/40f2d3ed-3d13-44bd-b5a0-7fe826182297.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/a4aff68c-105a-4e2e-93e1-a93487e3f53e.png

如果第一个子区间内符号发生变化,则 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/aeb1a071-16be-4a0e-8dc7-06b3fb165cb3.png 将重新定义为:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/55964a4d-c0f2-4cac-9220-584b8c3f556f.png

否则,它将以相同的方式重新定义为:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/da78bd63-e2ee-44d6-a9e2-7ce4c207a3b2.png

该过程会重复进行,直到区间的长度 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/ac340bd1-abae-44ba-b90f-af0835bce7e6.png 小于给定的容差。

例 5:可以使用欧几里得算法计算两个整数的最大公约数,该算法通过以下递归实现:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/4a74e62b-aed2-4ba1-8520-8a46208dac88.png

编写一个函数,用于计算两个整数的最大公约数。再编写一个函数,利用以下关系计算这两个数的最小公倍数:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/31b854bc-725e-4380-8491-1321cefce918.png

例 6:研究切比雪夫多项式的递归实现。参考第 7.4 节中的示例:递归函数。将程序改写为非递归形式,并研究计算时间与多项式次数的关系(另见timeit模块)。

在数学中,当我们写 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/b1869c68-1acd-4903-8bde-c96d445836db.png 时,我们指的是一个数学对象,对于这个对象,我们知道很多初等微积分的方法。例如:

这些方法不仅适用于 sin,还适用于其他足够光滑的函数。然而,也有其他数学对象,例如数字 5,对于这些对象,这些方法是没有意义的。

具有相同方法的对象被归为抽象类,例如函数。可以应用于函数的每个语句和方法,特别适用于 sincos

此类的其他例子可能是有理数,存在分母和分子方法;区间,具有左边界和右边界方法;无限序列,我们可以询问它是否有极限,等等。

在这种情况下,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/796b437a-5221-4fa9-af25-8c1efb98be48.png 被称为实例。数学表达式 设 g 为一个函数… 在这种情况下被称为实例化。这里,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/b441010f-f760-498d-83b1-19bbbacf79ca.png 是函数的名称,这是可以分配给它的众多属性之一。另一个属性可能是它的定义域。

数学对象 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/a10ad167-28bd-4fc6-9c5f-8a658b818ccb.png,但我们也可以为 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/5b332529-ae82-48ca-86f5-8858c45f7ae7.png 定义特殊的方法。例如,我们可能会要求 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/4d4879f1-5b0b-4547-8323-dc47f88332e5.pngs 的系数。这些方法可以用来定义多项式类。由于多项式是函数,它们还继承了函数类的所有方法。

在数学中,我们经常使用相同的运算符符号来表示完全不同的运算。例如,在 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/b3530656-c09b-4bda-87e2-34e9dac3fec7.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/23af4c8e-d3ba-4c2a-b5be-9d1e88f67112.png 中,运算符符号 + 的含义是不同的。通过使用相同的符号,强调了与对应数学运算的相似性。我们通过将这些术语从面向对象编程引入数学示例中,来应用它们,如类、实例和实例化、继承、方法、属性以及运算符重载。

具体来说,在本章中,我们将涵盖以下主题:

  • 类的简介

  • 绑定方法和非绑定方法

  • 类属性和类方法

  • 子类和继承

  • 封装

  • 类作为装饰器

本章将展示这些概念如何在 Python 中使用,并从一些基础知识开始。

第九章:8.1 类简介

本节介绍了类的最常见术语及其在 Python 中的实现。

首先,我们设置一个指导性示例。

8.1.1 一个指导性示例:有理数

我们将通过有理数的示例来说明类的概念,即形如 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/66f8f5f2-953a-4695-ba5d-a6ed3021626e.png 的数字,其中 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/830cf1ec-1cf4-43f7-9771-6b50c9b2ed63.pnghttps://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/99434b20-95a8-41b5-b0fa-79fbd50f80f7.png 是整数。

以下图示给出了一个类声明的示例:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/cecb5307-8e17-4e70-92d6-ac4927583d45.png

图 8.1:类声明示例

在这里,我们仅使用有理数作为类概念的示例。若要在 Python 中进行有理数运算,请使用 Python 模块 fractions

8.1.2 定义类并创建实例

类的定义是通过一个块命令来完成的,该命令使用 class 关键字、类名以及块中的一些语句(见 图 8.1):

class RationalNumber: 
    pass

通过 r = RationalNumber() 创建该类的实例(或换句话说,创建一个 RationalNumber 类型的对象),并且查询 type(r) 返回的结果是 <class'__main__.RationalNumber'>。如果我们想检查一个对象是否是该类的实例,可以使用以下代码:

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

到目前为止,我们已经生成了一个类型为 RationalNumber 的对象,但该对象尚未包含任何数据。此外,尚未定义任何方法来对这些对象进行操作。接下来的章节将讨论这一主题。

8.1.3 __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 类型的对象。这条语句做了两件事:

  1. 它首先创建一个空对象 q

  2. 然后,它应用 __init__ 函数;也就是说,执行 q.__init__(10, 20)

__init__ 的第一个参数指的是新对象本身。在函数调用时,这个第一个参数被对象的实例替代。这适用于类的所有方法,而不仅仅是特殊方法 __init__。这个第一个参数的特殊作用体现在约定中,通常将其命名为 self

在前面的示例中,__init__ 函数定义了新对象的两个属性,numeratordenominator

8.1.4 属性和方法

使用类的一个主要原因是可以将对象归为一组并绑定到一个共同的对象上。当我们查看有理数时就看到了这一点;denominatornumerator 是两个我们绑定到 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*step*z.imag

一旦实例定义好,我们就可以设置、更改或删除该实例的属性。语法与普通变量相同:

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

更改或删除属性可能会产生不希望出现的副作用,甚至可能导致对象失效。我们将在第 8.2 节中学习更多内容:相互依赖的属性。由于函数也是对象,我们也可以将函数作为属性使用;它们被称为实例的方法:

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

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

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

同样,这个方法将 self 作为第一个(也是唯一的)参数,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 只接受一个参数。

8.1.5 特殊方法

特殊方法 __repr__ 使我们能够定义对象在 Python 解释器中的表现方式。对于有理数,这个方法的可能定义如下:

class RationalNumber:
    ...
    def __repr__(self): 
        return f'{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__!=qD__ne__
<=__le__<__lt__
>=__ge__>__gt__
()__call__[]__getitem__

表 8.1:一些 Python 运算符及其对应的类方法

为新类实现这些运算符被称为运算符重载。运算符重载的另一个示例是一个方法,用于检查两个有理数是否相等:

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异常。

为了实现操作![],其中![]是RationalNumber的一个实例,我们将__radd__定义为:

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

注意,__radd__会交换参数的顺序;selfRationalNumber类型的对象,而other是必须转换的对象。

模拟函数调用和可迭代对象的方法

使用类实例和括号或中括号()[],调用的是特殊方法__call____getitem__,使得该实例表现得像函数或可迭代对象;另见表 8.1

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

现在可以按如下方式使用:

p = Polynomial(...)    # Creating a polynomial object
p(3.) # value of p at 3.

特殊方法__getitem__是有意义的,如果类提供了迭代器(建议你在考虑以下示例之前,回顾第 9.2.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.3.2 节:“迭代器工具”)。

8.2 彼此依赖的属性

实例的属性可以通过简单地为其赋值来修改(或创建)。然而,如果其他属性依赖于刚刚改变的属性,最好同时修改它们。

为了证明这一点,我们考虑一个例子:假设我们定义一个类,该类通过三个给定点来定义平面三角形对象。第一次尝试建立这样一个类可能是如下所示:

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,对应的边ac不会自动更新,计算出的面积也会出错:

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

一种解决方法是定义一个方法,该方法在属性改变时执行;这样的一个方法称为设置器方法。相应地,你可能需要一个在请求属性值时执行的方法;这样的一个方法称为获取器方法。我们现在将解释这两种方法是如何定义的。

8.2.1 函数属性

特殊函数property将属性与获取器、设置器和删除器方法关联。它也可以用来为属性分配文档字符串:

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

我们继续之前的例子,使用设置器方法,再次考虑Triangle类。如果在该类的定义中包括以下语句,那么命令tr.B = <某些值>会调用设置器方法set_B

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

让我们相应地修改Triangle类:

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

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

del tr.B # raises an exception

使用下划线作为属性名的前缀是一种约定,用来表示这些属性不打算被直接访问。它们用于存储由设置器和获取器处理的属性的数据。这些属性并不像其他编程语言中的私有属性那样真正私密,它们只是没有直接访问的设计意图。

8.3 绑定与未绑定的方法

现在我们将更详细地查看作为方法的属性。我们来看一个例子:

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被分配为实例的值。将方法绑定到实例上,使得该方法可以作为函数使用。在此之前,它没有任何用途。类方法(我们将在第 8.4.2 节:“类方法”中讨论)在这一点上与之不同。

8.4 类属性和类方法

到目前为止,我们已经看到了绑定到类实例上的属性和方法。在这一节中,我们介绍类属性和类方法。它们允许在实例创建之前访问方法和数据。

8.4.1 类属性

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

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 tol=1e-5
N1.tol # 1.e-5
N2.tol # 1.e-4 
# N2.tol is now detached and therefore not altered

8.4.2 类方法

我们在第 8.3 节中看到过:绑定和未绑定的方法,方法要么绑定到类的实例,要么保持未绑定的方法状态。类方法则不同,它们始终是绑定方法。它们绑定到类本身。

我们将首先描述语法细节,然后给出一些示例,展示这些方法的用途。

为了表明一个方法是类方法,装饰器行应当出现在方法定义之前:

@classmethod

标准方法通过其第一个参数引用实例,而类方法的第一个参数引用的是类本身。按照约定,标准方法的第一个参数称为 self,类方法的第一个参数称为 cls

以下是标准情况的示例:

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

这与 classmethod 的示例对比:

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 允许我们通过插值点定义一个多项式。

即使没有 Polynomial 的实例,我们也可以将插值数据转化为多项式系数:

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

print(p1 == p2)  # prints True

第 8.7 节中展示了类方法的另一个示例:类作为装饰器。在那里,类方法用于访问与该类的多个(或所有)实例相关的信息。

8.5 子类和继承

在这一节中,我们介绍一些面向对象编程中的核心概念:抽象类子类继承。为了帮助你理解这些概念,我们考虑了另一个数学示例:求解微分方程的单步方法。

普通初值问题的通用形式如下:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/984b6f95-b6c1-45b0-b4b9-a33b4a888ee6.png

数据包括右侧函数 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/ae32d475-3b7f-42fc-ba19-35f6c70f1e90.png,初始值 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/110bf954-0df2-465d-a569-e653d9b7bfec.png,以及感兴趣的区间 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/84ce4328-754a-4cae-a02c-6149358e7cf5.png

该问题的解是一个函数https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/6564bfaf-0da4-44ae-ae55-6e5e89f4e9ca.png,这些离散值 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/c623b321-35c6-4b22-9667-2847e637a603.png 是对https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/b47fccf1-7514-4e00-bfb2-adfc999d7379.png 的离散化值,在物理模型中,通常表示时间。

一步法通过递归步骤构造解的值 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/39df427c-bf50-4e98-9e2d-be6e53c2b29f.png

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/b8b77557-6e14-4250-951f-c876172b3cad.png

在这里,https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/c29823d2-ff3c-45aa-a1c1-b0ca44432554.png是一个步进函数,用于描述各个方法,详情请参见[28]

我们在这里所做的,是描述数学算法的典型方式。我们首先通过其思想描述了一种方法,以抽象的方式给出了步骤。要实际使用它,我们必须填写具体方法的参数,在这个例子中,就是函数![]。这也是面向对象编程中常见的解释方式。首先,我们设定一个类,提供方法的抽象描述:

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()

你可以通过使用星号操作符避免重复编写常见的参数列表(有关更多细节,请参见第 7.2.5 节: 可变参数数量):

...
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则允许我们在不更改所有父类引用的情况下更改父类的名称。

8.6 封装

有时使用继承是不实际的,甚至是不可能的。这促使我们使用封装。

我们将通过考虑 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)更为通用。我们现在可以通过继承来派生子类。

以切比雪夫多项式为例。它们可以在区间![]内计算:

![]

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

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

切比雪夫多项式是正交的,意思是:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/5828d012-2424-4a0c-9c78-6055184edd76.png

这可以通过这种构造轻松检查:

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这样的乘法函数是不可能实现的。

8.7 类作为装饰器

在第 7.8 节:函数作为装饰器中,我们看到了如何通过将另一个函数作为装饰器来修改函数。在第 8.1.5 节:特殊方法中,我们看到只要类提供了__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

8.8 总结

现代计算机科学中最重要的编程概念之一是面向对象编程。在本章中,我们学习了如何将对象定义为类的实例,并为其提供方法和属性。方法的第一个参数,通常表示为self,在其中扮演着重要且特殊的角色。你会看到一些方法,它们可以用来为自定义类定义基本运算,比如+*

虽然在其他编程语言中,属性和方法可以防止被意外使用,但 Python 允许一种技巧来隐藏属性,并通过特殊的 getter 和 setter 方法访问这些隐藏的属性。为此,你会遇到一个重要的函数property

8.9 练习

  1. RationalNumber类编写一个 simplify方法。该方法应返回该分数的简化版本,形式为一个元组。

  2. 为了提供带有置信区间的结果,数值数学中引入了一种特殊的计算方法,称为区间算术。定义一个名为Interval的类,并为其提供加法、减法、除法、乘法和幂运算(仅限正整数)的相关方法。这些运算遵循以下规则:

https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/f6b0c186-7c26-4084-8b79-d623f4f52a86.png

为这个类提供方法,使得可以进行a + I, a I, I + a, I a类型的运算,其中I是一个区间,a是整数或浮点数。首先,将整数或浮点数转换为区间[a,a]。(提示:你可能想使用函数装饰器来实现这一点;请参见章节 7.8:函数作为装饰器。)此外,实现__contains__方法,它使你能够使用x in I语法检查某个数字是否属于区间I,其中IInterval类型的对象。通过将多项式f=lambda x: 25*x**2-4*x+1应用于一个区间来测试你的类。

  1. 考虑第 8.7 节中的示例:类作为装饰器。扩展这个示例,创建一个函数装饰器,用于统计某个函数被调用的次数;另见第 7.8 节:函数作为装饰器

  2. 比较两种在RationalNumber类中实现反向加法__radd__的方法:一种是第 8.1.5 节中的示例:特殊方法,另一种是此处给出的实现:

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

你预期这个版本会出错吗?错误是什么,你如何解释它?通过执行以下代码测试你的答案:

q = RationalNumber(10, 15)
5 + q
  1. 考虑装饰器类CountCalls,如第 8.7 节中的示例:类作为装饰器。为这个类提供一个方法reset,该方法将字典CountCalls.instances中所有函数的计数器重置为0。如果将字典替换为空字典,会发生什么?
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值