原文：
annas-archive.org/md5/58a42b6d23877a7910fec539d6659c22译者：飞龙
测试
在本章中，我们将重点讨论科学编程中的两个方面的测试。第一个方面是科学计算中经常遇到的测试什么的问题。第二个方面涉及如何测试的问题。我们将区分手动测试和自动化测试。手动测试是每个程序员用来快速检查程序是否按预期工作的方式。自动化测试则是这一概念的精细化和自动化版本。我们将介绍一些适用于自动化测试的工具，并特别关注科学计算中的应用。
第十六章：15.1 手动测试
在代码开发过程中，你会做很多小的测试，以测试其功能。这可以称为手动测试。通常，你会通过在交互式环境中手动测试函数，来验证给定的函数是否做了它应该做的事。例如，假设你实现了二分法算法。它是一个找到标量非线性函数零点（根）的方法。为了启动该算法，必须给定一个区间，且该区间的边界上的函数值具有不同符号（详情见第 7.10 节：习题，了解更多信息）。
然后，你将测试该算法的实现，通常是通过检查：
-
当函数在区间边界处具有不同符号时，问题就得到了解决。
-
当函数在区间边界处具有相同符号时，是否抛出异常。
手动测试，尽管它看起来是必要的，但并不令人满意。一旦你确信代码做了它应该做的事，你会编写相对较少的示范示例来说服他人代码的质量。在那个阶段，你通常会对开发过程中进行的测试失去兴趣，它们可能会被遗忘或甚至删除。每当你更改了某个细节，导致事情不再正常工作时，你可能会后悔之前的测试已经不再可用。
15.2 自动化测试
开发任何代码的正确方法是使用自动化测试。其优势在于：
-
每次代码重构后，以及在发布任何新版本之前，自动化地重复大量测试。
-
对代码使用的静默文档记录。
-
记录代码的测试覆盖率：在更改之前，事情是否正常工作？某个特定方面是否从未经过测试？
程序中的变化，特别是其结构上的变化，但不影响其功能，称为代码重构。
我们建议在编码的同时开发测试。良好的测试设计本身就是一门艺术，而且很少有投资能像投资于良好测试那样，保证在开发时间节省方面获得如此好的回报。
现在，我们将通过考虑自动化测试方法来实现一个简单的算法。
15.2.1 测试二分法算法
让我们来研究一下二分法算法的自动化测试。通过该算法，可以找到一个实值函数的零点。它在第 7.10 节的练习 4中有所描述：练习。该算法的实现可以具有以下形式：
def bisect(f, a, b, tol=1.e-8):
"""
Implementation of the bisection algorithm
f real valued function
a,b interval boundaries (float) with the property
f(a) * f(b) <= 0
tol tolerance (float)
"""
if f(a) * f(b)> 0:
raise ValueError("Incorrect initial interval [a, b]")
for i in range(100):
c = (a + b) / 2.
if f(a) * f(c) <= 0:
b = c
else:
a = c
if abs(a - b) < tol:
return (a + b) / 2
raise Exception('No root found within the given tolerance {tol}')
我们假设该内容存储在名为bisection.py的文件中。作为第一个测试用例，我们测试该函数的零点是否能够找到！。
def test_identity():
result = bisect(lambda x: x, -1., 1.)
expected = 0\.
assert allclose(result, expected),'expected zero not found'
test_identity()
在这段代码中，您第一次接触到 Python 关键字assert。如果它的第一个参数返回False，它将引发AssertionError异常。它的可选第二个参数是一个包含附加信息的字符串。我们使用allclose函数来测试浮点数的相等性。
让我们对测试函数的某些特性进行评论。我们使用断言来确保如果代码没有按预期行为执行，异常将被抛出。我们必须手动运行测试，在test_identity()这一行中。
有许多工具可以自动化这种调用；我们将在第 15.2.2 节中看到其中的一种：使用 unittest 模块。
现在我们设置第二个测试，检查当函数在区间两端具有相同符号时，bisect是否会抛出异常。现在，我们假设抛出的异常是ValueError。在以下示例中，我们将检查初始区间！。对于二分法算法，它应该满足符号条件：
def test_badinput():
try:
bisect(lambda x: x,0.5,1)
except ValueError:
pass
else:
raise AssertionError()
test_badinput()
在这种情况下，如果抛出的异常不是ValueError类型，将引发AssertionError。有一些工具可以简化之前的构造，以检查是否抛出了异常。
另一个有用的测试是边界情况测试。在这里，我们测试可能会产生数学上未定义的情况或程序员未预见的程序状态的参数或用户输入。例如，如果两个边界相等，会发生什么？如果！，会发生什么？
以下代码是此类边界测试的示例：
def test_equal_boundaries():
result = bisect(lambda x: x, 0., 0.)
expected = 0.
assert allclose(result, expected), \
'test equal interval bounds failed'
def test_reverse_boundaries():
result = bisect(lambda x: x, 1., -1.)
expected = 0.
assert allclose(result, expected),\
'test reverse int_erval bounds failed'
test_equal_boundaries()
test_reverse_boundaries()
测试检查程序单元是否符合其规范的要求。在前面的例子中，我们假设规范要求在情况下！，这两个值应该默默地交换。并且这就是我们测试的内容。另一种方式是指定这种情况被视为错误输入，用户必须进行修正。在这种情况下，我们将测试是否抛出了适当的异常，例如ValueError。
15.2.2 使用 unittest 模块
Python 模块unittest大大简化了自动化测试。该模块要求我们重写之前的测试以保持兼容性。
第一个测试需要重写成一个class，如下所示：
from bisection import bisect
import unittest
class TestIdentity(unittest.TestCase):
def test(self):
result = bisect(lambda x: x, -1.2, 1.,tol=1.e-8)
expected = 0.
self.assertAlmostEqual(result, expected)
if __name__=='__main__':
unittest.main()
让我们看看与之前实现的区别。首先，测试现在是一个方法，并且是类的一部分。类必须继承自 unittest.TestCase。测试方法的名称必须以 test 开头。请注意，我们现在可以使用 unittest 包中的一个断言工具，即 assertAlmostEqual。最后，使用 unittest.main 运行测试。我们建议将测试写在与要测试的代码分开的文件中。因此，它以 import 开头。测试通过并返回如下：
Ran 1 test in 0.002s
OK
如果我们使用一个宽松的容差参数，例如 1.e-3，测试失败时将会报告：
F
======================================================================
FAIL: test (__main__.TestIdentity)
----------------------------------------------------------------------
Traceback (most recent call last):
File "<ipython-input-11-e44778304d6f>", line 5, in test
self.assertAlmostEqual(result, expected)
AssertionError: 0.00017089843750002018 != 0.0 within 7 places
----------------------------------------------------------------------
Ran 1 test in 0.004s
FAILED (failures=1)
测试可以并且应该作为测试类的方法进行分组，如下例所示：
import unittest
from bisection import bisect
class TestIdentity(unittest.TestCase):
def identity_fcn(self,x):
return x
def test_functionality(self):
result = bisect(self.identity_fcn, -1.2, 1.,tol=1.e-8)
expected = 0.
self.assertAlmostEqual(result, expected)
def test_reverse_boundaries(self):
result = bisect(self.identity_fcn, 1., -1.)
expected = 0.
self.assertAlmostEqual(result, expected)
def test_exceeded_tolerance(self):
tol=1.e-80
self.assertRaises(Exception, bisect, self.identity_fcn,
-1.2, 1.,tol)
if __name__=='__main__':
unittest.main()
在这里，在最后一个测试中我们使用了 unittest.TestCase.assertRaises 方法。它测试是否正确引发了异常。它的第一个参数是异常类型，例如 ValueError、Exception，第二个参数是预期引发异常的函数的名称。其余的参数是该函数的参数。
命令 unittest.main() 创建了一个 TestIdentity 类的实例，并执行那些以 test 开头的方法。
15.2.3 测试的 setUp 和 tearDown 方法
unittest.TestCase 类提供了两个特殊方法：setUp 和 tearDown，它们在每次调用测试方法之前和之后执行。这在测试生成器时很有用，因为生成器在每次测试后会被消耗完。我们通过测试一个程序来演示这一点，该程序检查文件中给定字符串首次出现的行：
class StringNotFoundException(Exception):
pass
def find_string(file, string):
for i,lines in enumerate(file.readlines()):
if string in lines:
return i
raise StringNotFoundException(
f'String {string} not found in File {file.name}.')
我们假设这段代码保存在名为 find_in_file.py 的文件中。
一个测试必须准备一个文件，打开并在测试后删除，如下例所示：
import unittest
import os # used for, for example, deleting files
from find_in_file import find_string, StringNotFoundException
class TestFindInFile(unittest.TestCase):
def setUp(self):
file = open('test_file.txt', 'w')
file.write('bird')
file.close()
self.file = open('test_file.txt', 'r')
def tearDown(self):
self.file.close()
os.remove(self.file.name)
def test_exists(self):
line_no=find_string(self.file, 'bird')
self.assertEqual(line_no, 0)
def test_not_exists(self):
self.assertRaises(StringNotFoundException, find_string,
self.file, 'tiger')
if __name__=='__main__':
unittest.main()
在每个测试之前执行 setUp，在每个测试之后执行 tearDown。
在创建测试用例时设置测试数据
方法 setUp 和 tearDown 在每个测试方法的前后执行。这在测试方法会改变数据时是必需的。它们保证在执行下一个测试之前，测试数据能够恢复到原始状态。
然而，也经常会有一种情况，测试不改变测试数据，你希望通过仅一次设置数据来节省时间。这可以通过类方法 setUpClass 来完成。
以下代码块简要说明了 setUpClass 方法的使用。你也许想再次查看 第 8.4 节：类属性和类方法。
import unittest
class TestExample(unittest.Testcase):
@classmethod
def setUpClass(cls):
cls.A=....
def Test1(self):
A=self.A
# assert something
....
def Test2(self):
A=self.A
# assert something else
15.2.4 测试参数化
我们经常希望用不同的数据集重复相同的测试。当使用 unittest 功能时，这要求我们自动生成带有相应方法注入的测试用例。
为此，我们首先构建一个测试用例，其中包含一个或多个将在后续设置测试方法时使用的方法。我们将再次考虑二分法，并检查它返回的值是否真的是给定函数的零点。
我们首先构建测试用例和将用于测试的方法，如下所示：
class Tests(unittest.TestCase):
def checkifzero(self,fcn_with_zero,interval):
result = bisect(fcn_with_zero,*interval,tol=1.e-8)
function_value=fcn_with_zero(result)
expected=0.
self.assertAlmostEqual(function_value, expected)
然后，我们动态创建测试函数作为该类的属性：
test_data=[
{'name':'identity', 'function':lambda x: x,
'interval' : [-1.2, 1.]},
{'name':'parabola', 'function':lambda x: x**2-1,
'interval' :[0, 10.]},
{'name':'cubic', 'function':lambda x: x**3-2*x**2,
'interval':[0.1, 5.]},
]
def make_test_function(dic):
return lambda self :\
self.checkifzero(dic['function'],dic['interval'])
for data in test_data:
setattr(Tests, f"test_{data['name']}", make_test_function(data))
if __name__=='__main__':
unittest.main()
在此示例中，数据以字典列表的形式提供。函数make_test_function动态生成一个测试函数，使用特定的数据字典与之前定义的checkifzero方法进行测试。最后，使用命令setattr将这些测试函数作为Tests类的方法。
15.2.5 断言工具
在本节中，我们收集了用于引发AssertionError的最重要工具。我们看到了命令assert和unittest中的三个工具，分别是assertAlmostEqual、assertEqual和assertRaises。以下表格（表 15.1）总结了最重要的断言工具及相关模块：
| 断言工具和应用示例 | 模块 |
|---|---|
assert 5 == 5 | – |
assertEqual(5.27, 5.27) | unittest.TestCase |
assertAlmostEqual(5.24, 5.2, places = 1) | unittest.TestCase |
assertTrue(5 > 2) | unittest.TestCase |
assertFalse(2 < 5) | unittest.TestCase |
assertRaises(ZeroDivisionError, lambda x: 1/x, 0.) | unittest.TestCase |
assertIn(3, {3, 4}) | unittest.TestCase |
assert_array_equal(A, B) | numpy.testing |
assert_array_almost_equal(A, B, decimal=5) | numpy.testing |
assert_allclose(A, B, rtol=1.e-3, atol=1.e-5) | numpy.testing |
表 15.1：Python、unittest 和 NumPy 中的断言工具
15.2.6 浮点数比较
两个浮点数不应使用==运算符进行比较，因为计算结果通常由于舍入误差略有偏差。为测试目的，存在许多用于测试浮点数相等性的工具。
首先，allclose检查两个数组是否几乎相等。它可以在测试函数中使用，如下所示：
self.assertTrue(allclose(computed, expected))
这里，self指的是一个unittest.Testcase实例。numpy包中的testing模块也有测试工具。可以通过以下方式导入：
import numpy.testing
测试两个标量或两个数组是否相等，可以使用numpy.testing.assert_array_allmost_equal或numpy.testing.assert_allclose。这两种方法在描述所需精度的方式上有所不同，如上表表 15.1所示。
https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/561e22e7-9410-4f8e-98c0-5f5472a83c2e.png因式分解将给定矩阵分解为一个正交矩阵https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3ec22240-f164-4c4c-91e2-08c4baac3ee3.png和一个上三角矩阵https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/79db8e65-2824-4d59-8a3a-731f50166e27.png，如下所示的例子所示：
import scipy.linalg as sl
A=rand(10,10)
[Q,R]=sl.qr(A)
方法应用是否正确？我们可以通过验证https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/daa3b63f-d531-4c68-af13-88bd2340ad86.png确实是一个正交矩阵来检查：
import numpy.testing as npt
npt.assert_allclose(
Q.T @ self.Q,identity(Q.shape[0]),atol=1.e-12)
此外，我们可能会通过检查是否进行了一项基本检查来执行一个合理性测试，方法是检查 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/93b48b98-bfcb-48d8-863c-e9a6af4c7585.png：
import numpy.testing as npt
npt.assert_allclose(Q @ R,A))
所有这些可以通过如下方式汇总到测试用例unittest中：
import unittest
import numpy.testing as npt
from scipy.linalg import qr
from scipy import *
class TestQR(unittest.TestCase):
def setUp(self):
self.A=rand(10,10)
[self.Q,self.R]=qr(self.A)
def test_orthogonal(self):
npt.assert_allclose(
self.Q.T @ self.Q,identity(self.Q.shape[0]),
atol=1.e-12)
def test_sanity(self):
npt.assert_allclose(self.Q @ self.R,self.A)
if __name__=='__main__':
unittest.main()
请注意，assert_allclose中的参数atol默认为零，这在处理具有小元素的矩阵时常常会导致问题。
15.2.7 单元测试与功能性测试
到目前为止，我们只使用了功能性测试。功能性测试检查功能是否正确。对于二分法算法，当存在零时，该算法实际上会找到它。在这个简单的例子中，什么是单元测试并不完全清楚。虽然看起来有点牵强，但仍然可以为二分法算法编写单元测试。它将展示单元测试如何通常导致更具模块化的实现。
因此，在二分法方法中，我们想要检查，例如，在每一步是否正确选择了区间。如何做到这一点呢？请注意，由于当前实现将算法隐藏在函数内部，这实际上是不可能的。一个可能的解决方法是仅运行一次二分法算法的步骤。由于所有步骤相似，我们可以认为我们已经测试了所有可能的步骤。我们还需要能够检查算法当前步骤中的当前边界a和b。因此，我们必须将要运行的步骤数作为参数添加，并改变函数的返回接口。我们将按如下方式进行：
def bisect(f,a,b,n=100):
...
for iteration in range(n):
...
return a,b
注意，我们必须更改现有的单元测试以适应这一变化。现在我们可以添加一个单元测试，如下所示：
def test_midpoint(self):
a,b = bisect(identity,-2.,1.,1)
self.assertAlmostEqual(a,-0.5)
self.assertAlmostEqual(b,1.)
15.2.8 调试
在测试时，有时需要调试，特别是当不能立即明确为何某个测试未通过时。在这种情况下，能够在交互式会话中调试特定测试是非常有用的。然而，由于unittest.TestCase类的设计，这使得测试用例对象的实例化变得不容易。解决方案是只为调试目的创建一个特殊实例。
假设在之前的TestIdentity类的示例中，我们想要测试test_functionality方法。可以按如下方式实现：
test_case = TestIdentity(methodName='test_functionality')
现在这个测试可以单独运行，命令是：
test_case.debug()
这将运行这个单独的测试，并允许进行调试。
15.2.9 测试发现
如果你编写一个 Python 包，多个测试可能分散在包的各个部分。模块discover会找到、导入并运行这些测试用例。命令行中的基本调用方式是：
python -m unittest discover
它开始在当前目录查找测试用例，并递归目录树向下查找名称中包含'test'字符串的 Python 对象。该命令接受可选参数。最重要的参数是-s来修改起始目录，-p来定义识别测试的模式：
python -m unittest discover -s '.' -p 'Test*.py'
15.3 测量执行时间
为了做出代码优化决策，通常需要比较几种代码替代方案，并根据执行时间决定优先使用哪种代码。此外，在比较不同算法时，讨论执行时间是一个重要问题。在本节中，我们展示了一种简单易用的计时方法。
15.3.1 使用魔法函数进行计时
测量单个语句的执行时间最简单的方法是使用 IPython 的魔法函数 %timeit。
IPython 外壳为标准 Python 添加了额外的功能。这些额外的功能被称为魔法函数。
由于单个语句的执行时间可能非常短，因此将语句放入循环中并执行多次。通过取最小的测量时间，确保计算机上其他任务不会对测量结果产生过多影响。
让我们考虑提取数组中非零元素的四种替代方法，如下所示：
A=zeros((1000,1000))
A[53,67]=10
def find_elements_1(A):
b = []
n, m = A.shape
for i in range(n):
for j in range(m):
if abs(A[i, j]) > 1.e-10:
b.append(A[i, j])
return b
def find_elements_2(A):
return [a for a in A.reshape((-1, )) if abs(a) > 1.e-10]
def find_elements_3(A):
return [a for a in A.flatten() if abs(a) > 1.e-10]
def find_elements_4(A):
return A[where(0.0 != A)]
使用 IPython 的魔法函数 %timeit 测量时间的结果如下：
In [50]: %timeit -n 50 -r 3 find_elements_1(A)
50 loops, best of 3: 585 ms per loop
In [51]: %timeit -n 50 -r 3 find_elements_2(A)
50 loops, best of 3: 514 ms per loop
In [52]: %timeit -n 50 -r 3 find_elements_3(A)
50 loops, best of 3: 519 ms per loop
In [53]: %timeit -n 50 -r 3 find_elements_4(A)
50 loops, best of 3: 7.29 ms per loop
参数 -n 控制在计时前语句执行的次数，-r 参数控制重复次数。
15.3.2 使用 Python 模块 timeit 进行计时
Python 提供了 timeit 模块，可用于测量执行时间。它要求首先构造一个时间对象。该对象是通过两个字符串构造的：一个包含设置命令的字符串和一个包含要执行命令的字符串。
我们采用与前面示例中相同的四种替代方案。现在，数组和函数的定义写入一个名为 setup_statements 的字符串，并按照如下方式构造四个计时对象：
import timeit
setup_statements="""
from scipy import zeros
from numpy import where
A=zeros((1000,1000))
A[57,63]=10.
def find_elements_1(A):
b = []
n, m = A.shape
for i in range(n):
for j in range(m):
if abs(A[i, j]) > 1.e-10:
b.append(A[i, j])
return b
def find_elements_2(A):
return [a for a in A.reshape((-1,)) if abs(a) > 1.e-10]
def find_elements_3(A):
return [a for a in A.flatten() if abs(a) > 1.e-10]
def find_elements_4(A):
return A[where( 0.0 != A)]
"""
experiment_1 = timeit.Timer(stmt = 'find_elements_1(A)',
setup = setup_statements)
experiment_2 = timeit.Timer(stmt = 'find_elements_2(A)',
setup = setup_statements)
experiment_3 = timeit.Timer(stmt = 'find_elements_3(A)',
setup = setup_statements)
experiment_4 = timeit.Timer(stmt = 'find_elements_4(A)',
setup = setup_statements)
定时器对象有一个 repeat 方法。它接受两个参数 repeat 和 number。该方法在一个循环中执行定时器对象的语句，测量时间，并根据 repeat 参数重复该实验。
我们继续前面的示例，并按照如下方式测量执行时间：
t1 = experiment_1.repeat(3,5)
t2 = experiment_2.repeat(3,5)
t3 = experiment_3.repeat(3,5)
t4 = experiment_4.repeat(3,5)
# Results per loop in ms
min(t1)*1000/5 # 615 ms
min(t2)*1000/5 # 543 ms
min(t3)*1000/5 # 546 ms
min(t4)*1000/5 # 7.26 ms
与使用 timeit 方法的示例不同，我们获取了所有测量结果的列表。由于计算时间可能根据计算机的整体负载而有所不同，因此该列表中的最小值可以视为执行语句所需计算时间的良好近似。
15.3.3 使用上下文管理器进行计时
最后，我们介绍第三种方法。它展示了上下文管理器的另一种应用。我们首先构造一个用于测量经过时间的上下文管理器对象，如下所示：
import time
class Timer:
def __enter__(self):
self.start = time.time()
# return self
def __exit__(self, ty, val, tb):
end = time.time()
self.elapsed=end-self.start
print(f'Time elapsed {self.elapsed} seconds')
return False
回顾一下，__enter__ 和 __exit__ 方法使得这个类成为一个上下文管理器。__exit__ 方法的参数 ty、val 和 tb 在正常情况下为 None。如果在执行过程中抛出异常，它们将分别包含异常类型、异常值和追踪信息。返回值 False 表示异常尚未被捕获。
我们将展示如何使用上下文管理器来测量前一个例子中四种方法的执行时间：
with Timer(): find_elements_1(A)
这样将显示类似 Time elapsed 15.0129795074 ms 的消息。
如果需要将计时结果存储在变量中，enter 方法必须返回 Timer 实例（取消注释 return 语句），并且必须使用 with ... as ... 结构：
with Timer() as t1:
find_elements_1(A)
t1.elapsed # contains the result
15.4 总结
没有测试就没有程序开发！我们展示了组织良好且有文档的测试的重要性。一些专业人士甚至通过首先指定测试来开始开发。一个有用的自动测试工具是 unittest 模块，我们已详细解释。虽然测试可以提高代码的可靠性，但性能提升需要通过分析代码的瓶颈来实现。不同的编码方法可能导致性能差异很大。我们展示了如何测量计算时间以及如何定位代码中的瓶颈。
15.5 习题
例 1： 如果存在一个矩阵 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/2a9104de-2907-4372-acd8-224c1bc3e52c.png，使得 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/9eb00513-26a8-4b51-bc70-69e23ce0e356.png，则两个矩阵 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/2d8e6f34-43b2-4e70-82aa-0403dfc9d4e9.png 被称为相似。矩阵 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/bda68266-6133-47e1-af9f-7ceb31af3bb3.png 和 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/b1ee0969-d1f6-4ca0-b59e-952e73467c51.png 具有相同的特征值。编写一个测试，检查两个矩阵是否相似，通过比较它们的特征值。这个是功能测试还是单元测试？
例 2： 创建两个大维度的向量。比较计算它们点积的不同方法的执行时间：
-
SciPy 函数：
v @ w -
生成器和求和：
sum((x*y for x,y in zip(v,w))) -
综合列表和求和：
sum([x*y for x,y in zip(v,w)])
例 3： 设！为一个向量。该向量 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/bc63c3f4-50dc-4b1c-be52-85d26cfad5bb.png的分量为
被称为 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/8ac291cf-4704-4bef-8bc3-988982c3b547.png的移动平均。确定计算 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/76df3e80-0678-45d1-9303-c3644ca27e70.png的两个方法中哪个更快：
v = (u[:-2] + u[1:-1] + u[2:]) / 3
或者
v = array([(u[i] + u[i + 1] + u[i + 2]) / 3
for i in range(len(u)-3)])
符号计算 - SymPy
本章将简要介绍如何使用 Python 进行符号计算。市场上有许多强大的软件可以执行符号计算，例如 Maple™或 Mathematica™。但有时，可能更倾向于在自己熟悉的语言或框架中进行符号计算。在本书的这一阶段，我们假设这种语言是 Python，因此我们寻找一个 Python 工具——模块 SymPy。
对 SymPy 的完整描述，如果有可能的话，将填满整本书，而这并非本章的目的。相反，我们将通过一些指导性示例来为您指明通向该工具的道路，展示其作为 NumPy 和 SciPy 的补充工具的潜力。
第十七章：16.1 什么是符号计算？
迄今为止，本书中所有的计算都属于所谓的数值计算。这些是主要基于浮动点数的一系列操作。数值计算的特点是，结果是精确解的近似值。
符号计算通过对公式或符号进行变换，像代数或微积分中所教的那样，将其转化为其他公式。这些变换的最后一步可能需要插入数字并进行数值计算。
我们通过计算这个定积分来说明两者之间的区别：
从符号上看，通过考虑被积函数的原始函数，可以对这个表达式进行变换：
我们现在通过插入积分界限来获得定积分的公式：
这被称为积分的封闭形式表达式。极少数数学问题有解可以用封闭形式表达式表示。这是积分的精确值，没有任何近似。此外，通过将实数表示为浮动点数，也不会引入误差，否则会产生舍入误差。
近似和舍入误差在最后一刻发挥作用，当需要对这个表达式进行求值时。平方根和arctan只能通过数值方法近似求值。这样的求值给出的是最终结果，精度达到某个特定的（通常是未知的）程度：
另一方面，数值计算会直接通过某种近似方法（例如辛普森法则）来近似定积分，并给出数值结果，通常还会给出误差估计。在 Python 中，这可以通过以下命令完成：
from scipy.integrate import quad
quad(lambda x : 1/(x**2+x+1),a=0, b=4)
它们返回值![]和误差边界的估计![]*。
以下图（图 16.1）展示了数值和符号近似的比较：
图 16.1：符号和数值求积
16.1.1 在 SymPy 中展开一个示例
首先，我们将展开之前的示例，并解释步骤。
首先，我们必须导入模块：
from sympy import *
init_printing()
第二个命令确保公式如果可能的话以图形方式呈现。接着，我们生成一个符号并定义被积函数：
x = symbols('x')
f = Lambda(x, 1/(x**2 + x + 1))
x 现在是一个类型为 Symbol 的 Python 对象，f 是一个 SymPy Lambda 函数（注意命令以大写字母开头）。
现在我们开始进行积分的符号计算：
integrate(f(x),x)
根据你的工作环境，结果的呈现方式不同；见下图 (图 16.2)，展示了在不同环境下 SymPy 公式的两种不同结果：
图 16.2：在两种不同环境下，SymPy 公式的两张截图
我们可以通过微分来检查结果是否正确。为此，我们为原始函数分配一个名称，并对其关于 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/1618222e-1bbe-4fb4-8342-4f4d2261bf0a.png 进行微分：
pf = Lambda(x, integrate(f(x),x))
diff(pf(x),x)
获得的结果如下：
这可以通过以下命令进行简化：
simplify(diff(pf(x),x))
到
这是我们预期的结果。
定积分通过以下命令获得：
pf(4) - pf(0)
简化后通过 simplify 命令得到的输出是：
为了获得数值结果，我们最终将这个表达式求值为一个浮动点数：
(pf(4)-pf(0)).evalf() # returns 0.9896614396123
16.2 SymPy 的基本元素
在这里，我们介绍了 SymPy 的基本元素。你会发现，如果已经熟悉 Python 中的类和数据类型，会更有利于理解。
16.2.1 符号 – 所有公式的基础
在 SymPy 中构建公式的基本构建元素是符号。正如我们在引言示例中看到的，符号是通过命令 symbols 创建的。这个 SymPy 命令通过给定的字符串生成符号对象：
x, y, mass, torque = symbols('x y mass torque')
这实际上是以下命令的简写：
symbol_list=[symbols(l) for l in 'x y mass torque'.split()]
然后进行解包操作以获取变量：
x, y, mass, torque = symbol_list
命令的参数定义了符号的字符串表示形式。所选符号的变量名通常与其字符串表示相同，但这并不是语言要求的：
row_index=symbols('i',integer=True)
print(row_index**2) # returns i**2
在这里，我们也定义了符号假定为整数。
可以通过非常紧凑的方式定义一整套符号：
integervariables = symbols('i:l', integer=True)
dimensions = symbols('m:n', integer=True)
realvariables = symbols('x:z', real=True)
类似地，可以通过以下方式定义索引变量的符号：
A = symbols('A1:3(1:4)')
这给出了一个符号元组：
索引的范围规则是我们在本书中使用切片时看到的规则（参见第 3.1.1 节：切片*，了解更多详细信息）。
16.2.2 数字
Python 直接对数字进行运算，并且引入了不可避免的舍入误差。这些误差会妨碍所有的符号计算。通过sympify数字可以避免这一问题：
1/3 # returns 0.3333333333333333
sympify(1)/sympify(3) # returns '1/3'
sympify命令将整数转换为类型为sympy.core.numbers.Integer的对象。
代替将1/3写作两个整数的运算，它也可以通过Rational(1,3)直接表示为一个有理数。
16.2.3 函数
SymPy 区分已定义函数和未定义函数。未定义函数（这一术语可能有些误导）指的是那些没有特殊性质的通用函数，虽然它们是已定义的 Python 对象。
一个具有特殊性质的函数例子是atan或本章入门示例中使用的Lambda函数。
请注意，同一数学函数的不同实现有不同的名称：sympy.atan和scipy.arctan。
未定义函数
通过给symbols命令一个额外的类参数，可以创建未定义函数的符号：
f, g = symbols('f g', cls=Function)
同样的效果可以通过使用构造函数Function来实现：
f = Function('f')
g = Function('g')
对于未定义的函数，我们可以评估微积分的通用规则。
例如，让我们评估以下表达式：
这是通过以下命令在 Python 中符号计算得到的：
x = symbols('x')
f, g = symbols('f g', cls=Function)
diff(f(x*g(x)),x)
执行时，前面的代码返回以下输出：
这个例子展示了如何应用乘积法则和链式法则。
我们甚至可以使用一个未定义的函数作为多个变量的函数，例如：
x = symbols('x:3')
f(*x)
这将返回以下输出：
请注意使用星号操作符来解包元组以形成带有参数的f；请参见第 7.2.5 节：可变数量的参数。
通过列表推导，我们可以构造一个包含所有偏导数的列表![]：
[diff(f(*x),xx) for xx in x]
这将返回一个列表，包含https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/6b258487-c295-41cf-89ad-f94b0fef08e5.png的梯度）：
该命令也可以通过使用Function对象的diff方法来重写：
[f(*x).diff(xx) for xx in x]
另一种方法是泰勒级数展开：
x = symbols('x')
f(x).series(x,0,n=4)
这将返回泰勒公式，并且包含通过兰道符号表示的余项：
16.2.4 基本函数
SymPy 中的基本函数例子包括三角函数及其反函数。以下例子展示了simplify如何作用于包含基本函数的表达式：
x = symbols('x')
simplify(cos(x)**2 + sin(x)**2) # returns 1
这是另一个使用基本函数的例子：
atan(x).diff(x) - 1./(x**2+1) # returns 0
如果你同时使用 SciPy 和 SymPy，我们强烈建议你将它们放在不同的命名空间中：
import numpy as np
import sympy as sym
# working with numbers
x=3
y=np.sin(x)
# working with symbols
x=sym.symbols('x')
y=sym.sin(x)
16.2.5 Lambda 函数
在第 7.7 节：匿名函数中，我们看到如何在 Python 中定义所谓的匿名函数。SymPy 的对应命令是Lambda。请注意两者的区别；lambda是一个关键字，而Lambda是一个构造函数。
命令Lambda接受两个参数，一个是函数的自变量符号，另一个是用于求值的 SymPy 表达式。
这里有一个示例，定义了空气阻力（也叫做拖力）作为速度的函数：
C,rho,A,v=symbols('C rho A v')
# C drag coefficient, A coss-sectional area, rho density
# v speed
f_drag = Lambda(v,-Rational(1,2)*C*rho*A*v**2)
f_drag以图形表达式的形式显示：
这个函数可以通过提供一个参数来以常规方式进行求值：
x = symbols('x')
f_drag(2)
f_drag(x/3)
这将得到以下表达式：
还可以通过仅提供Lambda的第一个参数为一个元组，来创建多个变量的函数，例如如下所示：
x,y=symbols('x y')
t=Lambda((x,y),sin(x) + cos(2*y))
调用这个函数有两种方式，可以通过直接提供多个参数来完成：
t(pi,pi/2) # returns -1
或者通过解包元组或列表：
p=(pi,pi/2)
t(*p) # returns -1
SymPy 中的矩阵对象甚至使我们能够定义向量值函数：
F=Lambda((x,y),Matrix([sin(x) + cos(2*y), sin(x)*cos(y)]))
这使得我们能够计算雅可比矩阵：
F(x,y).jacobian((x,y))
这将输出以下表达式：
如果有更多变量，使用更紧凑的形式来定义函数会更加方便：
x=symbols('x:2')
F=Lambda(x,Matrix([sin(x[0]) + cos(2*x[1]),sin(x[0])*cos(x[1])]))
F(*x).jacobian(x)
16.3 符号线性代数
符号线性代数由 SymPy 的matrix数据类型支持，我们将首先介绍它。然后我们将展示一些线性代数方法，作为符号计算在这一领域广泛应用的示例。
16.3.1 符号矩阵
我们在讨论向量值函数时简要介绍了matrix数据类型。在那里，我们看到了它最简单的形式，能够将一个列表的列表转换为矩阵。为了举个例子，让我们构造一个旋转矩阵：
phi=symbols('phi')
rotation=Matrix([[cos(phi), -sin(phi)],
[sin(phi), cos(phi)]])
在使用 SymPy 矩阵时，我们需要注意，操作符*执行的是矩阵乘法，而不是像 NumPy 数组那样的逐元素乘法。
之前定义的旋转矩阵可以通过使用矩阵乘法和矩阵的转置来检查其正交性：
simplify(rotation.T*rotation -eye(2)) # returns a 2 x 2 zero matrix
前面的示例展示了如何转置一个矩阵以及如何创建单位矩阵。或者，我们本可以检查其逆矩阵是否为其转置，方法如下：
simplify(rotation.T - rotation.inv())
设置矩阵的另一种方法是通过提供符号列表和形状：
M = Matrix(3,3, symbols('M:3(:3)'))
这将创建以下矩阵：
创建矩阵的第三种方法是通过给定函数生成其条目。语法如下：
Matrix(number of rows,number of colums, function)
我们通过考虑一个 Toeplitz 矩阵来举例说明前面的矩阵。它是一个具有常数对角线的矩阵。给定一个https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/296cf8da-ce9c-4a65-85c5-6acd0ae169c0.png，其元素定义为：
在 SymPy 中，可以直接使用此定义来定义矩阵：
def toeplitz(n):
a = symbols('a:'+str(2*n))
f = lambda i,j: a[i-j+n-1]
return Matrix(n,n,f)
执行前面的代码会得到toeplitz(5)：
我们可以清楚地看到所需的结构；所有沿副对角线和超对角线的元素都相同。我们可以根据 Python 语法（见第 3.1.1 节：切片）通过索引和切片访问矩阵元素：
M[0,2]=0 # changes one element
M[1,:]=Matrix(1,3,[1,2,3]) # changes an entire row
16.3.2 SymPy 中线性代数方法的示例
线性代数中的基本任务是求解线性方程组：
让我们用符号方法处理一个 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/980b1643-67ce-4e96-aa4e-5cfcf5d6595b.png 矩阵：
A = Matrix(3,3,symbols('A1:4(1:4)'))
b = Matrix(3,1,symbols('b1:4'))
x = A.LUsolve(b)
这个相对较小问题的输出已经只是可读的，我们可以在以下的图形表达式中看到：
同样，使用simplify命令有助于我们检测抵消项并收集公共因子：
simplify(x)
这将导致以下输出，结果看起来更好：
当矩阵维度增加时，符号计算变得非常慢。对于大于 15 的维度，甚至可能出现内存问题。
下一幅图（图 16.3）展示了符号求解和数值求解线性系统之间 CPU 时间的差异：
图 16.3：数值和符号求解线性系统的 CPU 时间
16.4 替代
我们首先考虑一个简单的符号表达式：
x, a = symbols('x a')
b = x + a
如果我们设置 x = 0 会发生什么？我们观察到 b 并没有改变。我们所做的是改变了 Python 变量 x，它现在不再引用符号对象，而是引用了整数对象 0。由字符串 'x' 表示的符号保持不变，b 也没有变化。
另一方面，通过将符号替换为数字、其他符号或表达式，来改变表达式是通过一种特殊的替代方法完成的，以下代码展示了这一点：
x, a = symbols('x a')
b = x + a
c = b.subs(x,0)
d = c.subs(a,2*a)
print(c, d) # returns (a, 2a)
此方法接受一个或两个参数。以下两个语句是等价的：
b.subs(x,0)
b.subs({x:0}) # a dictionary as argument
将字典作为参数允许我们一步完成多个替代：
b.subs({x:0, a:2*a}) # several substitutions in one
由于字典中的项目没有固定顺序——我们永远无法知道哪个是第一个——因此需要确保对项目的排列不会影响替代结果。因此，在 SymPy 中，替代操作首先在字典中进行，然后再在表达式中进行。以下示例演示了这一点：
x, a, y = symbols('x a y')
b = x + a
b.subs({a:a*y, x:2*x, y:a/y})
b.subs({y:a/y, a:a*y, x:2*x})
两种替代方法返回相同的结果：
![]
定义多个替代的第三种方法是使用旧值/新值对的列表：
b.subs([(y,a/y), (a,a*y), (x,2*x)])
也可以将整个表达式替换为其他表达式：
n, alpha = symbols('n alpha')
b = cos(n*alpha)
b.subs(cos(n*alpha), 2*cos(alpha)*cos((n-1)*alpha)-cos((n-2)*alpha))
为了说明矩阵元素的替代，我们再次取 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3a72dcd0-f73c-4b62-8d90-2b07c3477916.png Toeplitz 矩阵：
考虑替换T.subs(T[0,2],0)。它改变了位置[0, 2]处的符号对象，即符号https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/ac3eb574-96b2-42d8-9e91-627d68e3c16b.png。它还出现在其他两个地方，这些地方会被这个替换自动影响。
给定的表达式是结果矩阵：
或者，我们可以为该符号创建一个变量并在替换中使用它：
a2 = symbols('a2')
T.subs(a2,0)
作为一个更复杂的替换示例，让我们考虑如何将 Toeplitz 矩阵转换为三对角 Toeplitz 矩阵*.* 这可以通过以下方式完成：
首先，我们生成一个符号列表，选择我们要替换的符号；然后使用zip命令生成一对对的列表。最后，我们通过给出旧值/新值对的列表进行替换，如前所述：
symbs = [symbols('a'+str(i)) for i in range(19) if i < 3 or i > 5]
substitutions=list(zip(symbs,len(symbs)*[0]))
T.subs(substitutions)
这会得到以下矩阵作为结果：
16.5 评估符号表达式
在科学计算中，通常需要先进行符号运算，然后将符号结果转换为浮点数。
评估符号表达式的核心工具是evalf。它通过使用以下方法将符号表达式转换为浮点数：
pi.evalf() # returns 3.14159265358979
结果对象的数据类型是Float（注意大写），这是一个 SymPy 数据类型，允许使用任意位数（任意精度）的浮点数。
默认精度对应于 15 位数字，但可以通过给evalf一个额外的正整数参数来改变它。
通过指定所需精度的数字位数：
pi.evalf(30) # returns 3.14159265358979323846264338328
使用任意精度的一个结果是，数字可以非常小，也就是说，打破了经典浮点数表示的限制；请参见第 2.2.2 节：浮点数。
有趣的是，用Float类型的输入来评估 SymPy 函数会返回一个与输入精度相同的Float。我们将在一个来自数值分析的更复杂示例中演示这一事实的使用。
16.5.1 示例：牛顿法收敛阶的研究
一个迭代方法，如果迭代*https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/e551ea10-3819-42ef-b8b1-3884edbf01bc.png收敛，并且存在一个正的常数https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/d1e134fb-86be-4b63-ab18-e936219620fe.png，使得：
牛顿法，当从一个好的初值开始时，其收敛阶为https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/d5d4c2b0-7c21-4ae7-abe1-13d87cd15d2e.png，对于某些问题，甚至可以达到https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3778c99a-3238-4758-81e5-ffa04d05c1e4.png。应用牛顿法解决问题https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3e25a412-60af-451a-9075-b63f22edc726.png时，给出以下迭代方案：
该过程的收敛速度是立方收敛；也就是说，q = 3。
这意味着正确数字的数量会随着每次迭代从上一轮迭代中三倍增加。为了演示立方收敛并数值求解常数，![]使用标准的 16 位数字float数据类型几乎无法实现。
以下代码使用 SymPy 并结合高精度求值，将立方收敛研究推向极致：
import sympy as sym
x = sym.Rational(1,2)
xns=[x]
for i in range(1,9):
x = (x - sym.atan(x)*(1+x**2)).evalf(3000)
xns.append(x)
结果如下图所示（图 16.4），显示了每次迭代正确数字的数量是如何从上一轮迭代中三倍增加的：
图 16.4：对应用于https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/188cd499-6f66-4039-bcfd-e039c4de31da.png的牛顿法收敛性的研究
这种极高精度要求（3,000 位数字！）使我们能够以如下方式评估前面序列的七项，从而演示立方收敛：
import numpy as np
# Test for cubic convergence
print(np.array(np.abs(np.diff(xns[1:]))/np.abs(np.diff(xns[:-1]))**3,
dtype=np.float64))
[ 0.41041618, 0.65747717, 0.6666665, 0.66666667, 0.66666667, 0.66666667, 0.66666667]}
16.5.2 将符号表达式转换为数值函数
正如我们所见，符号表达式的数值求解分为三个步骤：首先，我们进行一些符号计算，然后通过数字替换变量，最后使用evalf进行浮动点数的求值。
进行符号计算的原因通常是我们希望进行参数研究。这要求在给定的参数范围内修改参数。这要求符号表达式最终被转换为数值函数。
对多项式系数的参数依赖性研究
我们通过一个插值示例展示了符号/数值参数研究，以介绍 SymPy 命令lambdify。
让我们考虑任务，即对数据https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/1186c128-483c-4d51-8df7-207d483a2da8.png 和 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/58696cc5-47bd-4b9d-a3f5-0c0344b1abda.png进行插值。在这里，https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/99eedf3c-a08b-4a54-85e3-061b29472f88.png是一个自由参数，我们将在区间https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/a32da2e1-a0f2-4cd7-aa5b-badedd42495b.png上变化。
二次插值多项式具有依赖于该参数的系数：
使用 SymPy 和练习 3中描述的单项式方法，如第 4.11 节中的练习给出了这些系数的封闭公式：
t=symbols('t')
x=[0,t,1]
# The Vandermonde Matrix
V = Matrix([[0, 0, 1], [t**2, t, 1], [1, 1,1]])
y = Matrix([0,1,-1]) # the data vector
a = simplify(V.LUsolve(y)) # the coefficients
# the leading coefficient as a function of the parameter
a2 = Lambda(t,a[0])
我们为插值多项式的主系数![]获得了一个符号函数：
现在是将表达式转换为数值函数的时候了，例如，为了生成一个图形。这是通过lamdify函数完成的。该函数接受两个参数，一个是自变量，另一个是 SymPy 函数。
在我们的 Python 示例中，我们可以编写：
leading_coefficient = lambdify(t,a2(t))
现在可以通过以下命令绘制该函数，例如：
import numpy as np
import matplotlib.pyplot as mp
t_list= np.linspace(-0.4,1.4,200)
ax=mp.subplot(111)
lc_list = [leading_coefficient(t) for t in t_list]
ax.plot(t_list, lc_list)
ax.axis([-.4,1.4,-15,10])
ax.set_xlabel('Free parameter $t$')
ax.set_ylabel('$a_2(t)$')
图 16.5 是此参数研究的结果，我们可以清楚地看到由于多个插值点（这里在 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/97299a92-07c9-427b-943c-7dcf074dcf70.png 或 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/266b98d5-623e-4ca0-9e10-2a00f2f33642.png）而产生的奇点：
图 16.5：多项式系数依赖于插值点位置的关系
16.6 小结
在本章中，你初步了解了符号计算的世界，并领略了 SymPy 的强大功能。通过学习这些例子，你掌握了如何设置符号表达式、如何操作符号矩阵，并学习了如何进行简化处理。通过处理符号函数并将其转化为数值计算，你最终建立了与科学计算和浮点结果之间的联系。你体验了 SymPy 的强大，它通过与 Python 的完美集成，提供了强大的构造和易读的语法。
将最后一章视为开胃菜，而非完整菜单。我们希望你对未来在科学计算和数学中的精彩编程挑战产生兴趣。
与操作系统交互
本节是 Python 简介的附加内容。它将 Python 放入计算机操作系统的上下文中，并展示了如何在命令窗口（也称为控制台）中使用 Python。我们演示系统命令和 Python 命令如何互动以及如何创建新应用程序。
在许多其他操作系统和各种 方言 中，桌面计算机和笔记本电脑上使用的三个主要操作系统是：Windows 10、macOS 和 Linux。在本章中，我们讨论如何在 Linux 环境中使用 Python。所有示例均在广泛分发的 Ubuntu 环境中进行测试。这里介绍的原则同样适用于其他大型操作系统。但本章内容仅限于 Linux，作为一个完全自由获取的环境。
首先，假设我们已经编写并测试了一个 Python 程序，并且现在希望直接从控制台窗口而不是从诸如 IPython 的 Python shell 运行它。
本章涵盖以下主题：
-
在 Linux shell 中运行 Python 程序
-
模块
sys -
如何从 Python 执行 Linux 命令
第十八章：17.1 在 Linux shell 中运行 Python 程序
打开终端应用程序时，会获得一个带有命令提示符的窗口；参见 图 17.1：
图 17.1：Ubuntu 20.04 中的终端窗口
终端窗口带有命令提示符，通常由用户名和计算机名称前缀，后跟目录名称。这取决于个人设置。
要执行名为 myprogram.py 的 Python 命令文件，您有两个选择：
-
执行命令
python myprogram.py -
直接执行命令
myprogram.py
第二种变体需要一些准备工作。首先，您必须允许执行该文件，其次，您必须告诉系统需要执行该文件的命令。通过以下命令来允许执行该文件：
chmod myprogram.py o+x
chmod 代表更改文件模式。该命令后跟文件名，最后是所需的新模式，在这里是 o+x。
在该示例中给出的模式代表“授予（+）文件所有者（o）执行（x）该文件的权限”。我们假设您是文件的所有者，并且它位于当前目录中。
然后，我们必须找到计算机上命令 python 的位置。通过运行以下命令来完成：
which python
如果您像在 第 1.1 节 中描述的那样通过 Anaconda 安装 Python，则会获得有关命令 python 在您的系统上位置的信息，例如 /home/claus/anaconda3/bin/python。
这些信息必须写在 Python 脚本的第一行，以告知系统通过哪个程序可以将文本文件转换为可执行文件。以下是一个示例：
#! /home/claus/anaconda3/bin/python
a=3
b=4
c=a+b
print(f'Summing up {a} and {b} gives {c}')
Python 视角中的第一行只是一个注释，会被忽略。但是 Linux 会将第一行中shebang组合符号#!之后的部分视为将文件转化为可执行文件所需的命令。这里，它学习使用位于目录/home/claus/anaconda3/bin/中的python命令来实现这一点。
我们可以通过逻辑路径来定义 Python 解释器的位置，而不是通过绝对路径，这样可以使代码更具可移植性。你可以将其提供给使用 Linux 系统的同事，而无需修改：
#! /usr/bin/env python # a logical path to Python (mind the blank)
现在你可以直接在控制台执行示例代码；见 图 17.2：
图 17.2：在 Linux 终端中执行示例文件 example.py
我们需要在命令前加上./，这告诉操作系统在当前目录中查找该命令。如果没有这个前缀，Linux 会在搜索路径中列出的某个目录中查找文件example.py。
在接下来的部分中，我们将展示如何直接从命令行传递参数给 Python 程序。
17.2 sys 模块
模块sys提供了与系统命令进行通信的工具。我们可以直接从命令行传递带有参数的 Python 脚本，并将结果输出到控制台。
17.2.1 命令行参数
为了说明命令行参数的使用，我们考虑以下代码段，我们将其保存到名为demo_cli.py的文件中：
#! /usr/bin/env python
import sys
text=f"""
You called the program{sys.argv[0]}
with the {len(sys.argv)-1} arguments, namely {sys.argv[1:]}"""
print(text)
在通过chmod o+x demo_cli.py给予文件执行权限后，我们可以在 Shell 中带参数执行它；见 图 17.3：
图 17.3：在终端命令行上带三个参数执行 Python 脚本
控制台中给出的三个参数可以通过列表sys.argv在 Python 脚本中访问。这个列表中的第一个元素——索引为0的元素——是脚本的名称。其他元素是给定的参数，作为字符串形式。
参数是传递给 Python 脚本的调用。它们不应与脚本执行过程中用户输入的内容混淆。
17.2.2 输入与输出流
在前面的示例中，我们使用了命令print来显示终端中生成的消息（甚至是在 Python Shell 中）。脚本的先前输入通过参数和变量sys.argv获取。与print相对的命令是input，它提示从终端（或 Python Shell）获取数据。
在第 14.1 节：文件处理中，我们看到了如何通过文件对象和相关方法向脚本提供数据以及如何输出数据。模块sys使得可以将键盘当作输入的文件对象（例如，readline、readlines），而将控制台当作输出的文件对象（例如，write、writelines）。
在 UNIX 中，信息流通过三个流组织：
-
标准输入流：
STDIN -
标准输出流：
STDOUT -
标准错误流：
STDERR
这些流对应于文件对象，可以通过sys.stdin、sys.stdout和sys.stderr在 Python 中访问。
为了举个例子，我们考虑一个小脚本pyinput.py，它计算一些数字的总和：
#!/usr/bin/env python3
from numpy import *
import sys
# Terminate input by CTRL+D
a=array(sys.stdin.readlines()).astype(float)
print(f'Input: {a}, Sum {sum(a)}'
语句sys.stdin.readlines()创建了一个生成器。命令array会迭代这个生成器，直到用户输入结束符号，该符号在 Linux 系统上是CTRL-D，在 Windows 系统上是CTRL-Z。请参见图 17.4中的截图：
图 17.4：使用sys.stdin执行脚本的终端截图。
请注意，给定的结束符号 CTRL-D 是不可见的。
重定向流
标准输入正在等待来自键盘的数据流。但输入可以通过文件进行重定向。通过在 Linux 中使用重定向符号<可以实现这一点。我们通过使用与之前相同的脚本，但这次由数据文件intest.txt提供数据来展示这一点；请参见图 17.5：
图 17.5：展示输入重定向（sys.stdin）的截图
脚本本身无需修改。无论哪种方式都可以使用。
输出也是如此。默认情况下，输出会显示在终端中，但这里也有将输出重定向到文件的选项。在这种情况下，重定向符号是>；请参见图 17.6：
图 17.6：重定向输入和重定向输出的截图
现在，创建了一个名为result.txt的文件，并将脚本的输出写入该文件。如果该文件已经存在，其内容将被覆盖。如果输出应该附加到已有文件的末尾，则必须使用重定向符号>>。
最后，有时可能希望将输出与错误或警告信息分开。这就是 Linux 提供两个输出流通道的原因，分别是sys.stdout和sys.stderr。默认情况下，两者都指向同一个位置，即终端。但通过使用重定向符号，错误信息可以例如写入文件，而主输出则显示在屏幕上，或者反之亦然。
为了演示这一点，我们修改示例pyinput.py，以便在没有提供输入时生成错误信息：
#!/usr/bin/env python3
from numpy import *
import sys
# Terminate input by CTRL+D
a=array(sys.stdin.readlines()).astype(float)
print(f'Input: {a}, Sum {sum(a)}')
if a.size == 0:
sys.stderr.write('No input given\n')
在终端窗口中，重定向输入和错误输出的脚本典型调用呈现在图 17.7中：
图 17.7：带有标准输入和标准错误输出重定向的终端窗口截图
如果输入文件为空，错误信息会写入文件 error.txt，而输出则显示在终端窗口中。
错误消息是来自未捕获异常的消息，包括语法错误和显式写入 sys.stderr 的文本。
在表 17.1中，汇总了不同的重定向情况：
| 任务 | Python 对象 | 重定向符号 | 替代符号 |
|---|---|---|---|
| 数据输入 | sys.stdin | < | |
| 数据输出 | sys.stdout | > | 1> |
| 数据输出附加到文件 | sys.stdout | >> | 1>> |
| 错误输出 | sys.stderr | 2> | |
| 错误输出附加到文件 | sys.stderr | 2>> | |
| 所有输出 | sys.stdout, sys.stderr | &> | |
| 所有输出附加到文件 | sys.stdout, sys.stderr | &>> |
表 17.1：不同重定向场景的汇总
在 Linux 命令与 Python 脚本之间建立管道
在上一节中，我们展示了如何将 Python 程序的输入和输出重定向到文件。当不同 Python 程序之间，或者 Python 程序与 Linux 命令之间的数据流时，数据通常通过文件传递。如果数据不在其他地方使用或应当保留供以后使用，这个过程就显得冗长：仅仅为了直接将信息从一段代码传递到另一段代码，需要创建、命名和删除文件。替代方案是使用 Linux 管道，让数据从一个命令直接流向另一个命令。
让我们从一个纯 Linux 示例开始，然后将管道构造应用于 Python。
Linux 命令 ifconfig 显示有关 Linux 计算机当前网络配置的大量信息。在这些信息中，你可以找到 IP 地址（即当前使用的网络地址）。例如，要自动判断计算机（如笔记本电脑）是否通过某个网络单元连接到家庭网络，而不是连接到外部网络，你可能需要扫描 ifconfig 的输出，查找包含网络适配器标识符（例如 wlp0s20f3）的行，并在其下几行中查找包含网络前缀的字符串（如 192.168）。如果找到这个字符串，输出应仅显示包含该字符串的行；即行数计数应返回 1。如果计算机未连接到家庭网络，则行数计数返回 0。
我们使用三个命令：
-
ifconfig用于显示完整的网络配置信息。 -
grep用于查找包含特定模式的行。该行及根据参数-A的要求，显示一些后续行。 -
wc用于对文件执行各种计数操作。参数-l指定计数行数。
这些命令的输出会直接传递给下一个命令。这是通过使用管道符号 | 来实现的：
ifconfig|grep -A1 wlp0s20f3 | grep 192.168|wc -l
此命令行在家庭网络中执行时仅在屏幕上显示 1。所有中间输出直接传递到下一个命令，而无需显示任何内容，并且不使用任何文件来临时存放信息，直到下一个命令读取它。一个命令的标准输出，stdout，成为下一个命令的标准输入，stdin。
这也适用于 Python 脚本直接调用。
我们通过在管道中演示 Python 脚本来继续上一个示例。在这里，我们简单地使用 Python 生成一个友好的消息：
#!/usr/bin/env python3
import sys
count = sys.stdin.readline()[0]
status = '' if count == '1' else 'not'
print(f"I am {status} at home")
现在，我们可以通过添加此脚本扩展管道；请参见 图 17.8 中的截图：
图 17.8：带有五个管道和一个 Python 脚本的命令链
在此示例中，我们在终端窗口中执行了一系列 Linux 程序和一个 Python 脚本。或者，可以让 Python 脚本调用 UNIX 命令。这将在下一节中演示。
17.3 如何从 Python 执行 Linux 命令
在上一节中，我们看到如何从 Linux 终端执行 Python 命令。在本节中，我们考虑如何在 Python 程序中执行 Linux 命令的重要情况。
17.3.1 模块 subprocess 和 shlex
要在 Python 中执行系统命令，首先需要导入模块 subprocess。此模块提供的高级工具是 run。使用此工具，您可以快速访问 Python 中的 Linux 命令，并处理它们的输出。
更复杂的工具是 Popen，我们将在解释如何在 Python 中模拟 Linux 管道时进行简要介绍。
完整的过程：subprocess.run
我们将通过最标准和简单的 UNIX 命令 ls 来演示此工具——列出目录内容的命令。它带有各种可选参数；例如，ls -l 以扩展信息显示列表。
要在 Python 脚本中执行此命令，我们使用 subprocess.run。最简单的用法是仅使用一个参数，将 Linux 命令拆分为几个文本字符串的列表：
import subprocess as sp
res = sp.run(['ls','-l'])
模块 shlex 提供了一个特殊工具来执行此分割：
_import shlex
command_list = shlex.split('ls -l') # returns ['ls', '-l']
它还尊重文件名中的空格，并且不将其用作分隔符。
命令 run 显示 Linux 命令的结果和 subprocess.CompletedProcess 对象 res。
以这种方式执行 UNIX 命令是相当无用的。大多数情况下，您希望处理输出。因此，必须将输出提供给 Python 脚本。为此，必须将可选参数 capture_output 设置为 True。通过这种方式，UNIX 命令的 stdout 和 stderr 流将作为返回对象的属性可用：
import subprocess as sp
import shlex
command_list = shlex.split('ls -l') # returns ['ls', '-l']
res = sp.run(command_list, capture_output=True)
print(res.stdout.decode())
注意，此处使用方法 decode 将字节字符串解码为标准字符串格式。
如果 Linux 命令返回错误，属性 res.returncode 将获得非零值，并且 res.stderr 包含错误消息。
Python 的做法是抛出一个错误。可以通过将可选参数check设置为True来实现：
import subprocess as sp
import shlex
command ='ls -y'
command_list = shlex.split(command) # returns ['ls', '-y']
try:
res = sp.run(command_list, capture_output=True, check=True)
print(res.stdout.decode())
except sp.CalledProcessError:
print(f"{command} is not a valid command")
创建进程：subprocess.Popen
当你对一个需要用户输入才能终止的进程应用subprocess.run时会发生什么？
这样一个程序的简单示例是xclock。它打开一个新窗口，显示一个时钟，直到用户关闭窗口。
由于命令subprocess.run创建了一个CompletedProcess对象，下面的 Python 脚本：
import subprocess as sp
res=sp.run(['xclock'])
启动一个进程并等待它结束，也就是说，直到有人关闭带有时钟的窗口；参见图 17.9：
图 17.9：xclock 窗口
这对subprocess.Popen有影响。它创建了一个 _Popen对象。进程本身变成了一个 Python 对象。它不需要完成才能成为一个可访问的 Python 对象：
import subprocess as sp
p=sp.Popen(['xclock'])
进程通过用户在时钟窗口上的操作或通过显式终止进程来完成，命令如下：
p.terminate()
使用Popen，我们可以在 Python 中构建 Linux 管道。下面的脚本：
import subprocess as sp_
p1=sp.Popen(['ls', '-l'],stdout=sp.PIPE)
cp2=sp.run(['grep','apr'], stdin=p1.stdout, capture_output=True)
print(cp2.stdout.decode())
对应于 UNIX 管道：
ls -l |grep 'apr
它显示了在四月最后访问的目录中的所有文件。
模块subprocess有一个对象PIPE，它接受第一个进程p1的输出。然后，它被作为stdin传递给命令run。该命令然后返回一个带有stdout属性（类型为bytes）的CompletedProcess对象cp2。最后，调用方法decode可以很好地打印结果。
另外，也可以通过使用两个Popen进程来实现相同的效果。第二个进程也使用管道，可以通过方法communicate进行显示：
import subprocess as sp
p1=sp.Popen(['ls', '-l'],stdout=sp.PIPE)
p2=sp.Popen(['grep','apr'], stdin=p1.stdout, stdout=sp.PIPE)
print(p2.communicate()[0].decode)
方法communicate返回一个元组，其中包含stdout和stderr上的输出。
17.4 小结
在本章中，我们演示了 Python 脚本与系统命令的交互。Python 脚本可以像系统命令一样被调用，或者 Python 脚本本身可以创建系统进程。本章基于 Linux 系统，如 Ubuntu，仅作为概念和可能性的演示。它允许将科学计算任务置于应用场景中，在这些场景中，通常需要将不同的软件结合起来。甚至硬件组件可能也会涉及其中。
用于并行计算的 Python
本章涉及并行计算和模块mpi4py。复杂且耗时的计算任务通常可以拆分为子任务，如果有足够的计算能力，这些子任务可以同时执行。当这些子任务彼此独立时，进行并行计算尤其高效。需要等待另一个子任务完成的情况则不太适合并行计算。
考虑通过求积法则计算一个函数的积分任务：
使用![]。如果评估![]非常耗时，且![]很大，那么将问题拆分为两个或多个较小的子任务会更有利：
我们可以使用几台计算机，并将必要的信息提供给每台计算机，使其能够执行各自的子任务，或者我们可以使用一台带有所谓多核架构的计算机。
一旦子任务完成，结果会传送给控制整个过程并执行最终加法的计算机或处理器。
我们将在本章中以此为指导示例，涵盖以下主题：
-
多核计算机和计算机集群
-
消息传递接口（MPI）
第十九章：18.1 多核计算机和计算机集群
大多数现代计算机都是多核计算机。例如，本书写作时使用的笔记本电脑配备了 Intel® i7-8565U 处理器，具有四个核心，每个核心有两个线程。
这意味着什么？处理器上的四个核心允许并行执行四个计算任务。四个核心，每个核心有两个线程，通常被系统监视器计为八个 CPU。在本章中，只有核心数量才是重要的。
这些核心共享一个公共内存——你的笔记本的 RAM——并且每个核心有独立的缓存内存：
图 18.1：具有共享和本地缓存内存的多核架构
缓存内存由核心最优使用，并以高速访问，而共享内存可以被一个 CPU 的所有核心访问。在其上方是计算机的 RAM 内存，最后是硬盘，也是共享内存。
在下一部分，我们将看到如何将计算任务分配到各个核心，以及如何接收结果并进一步处理，例如，存储到文件中。
另一种并行计算的设置是使用计算机集群。在这里，一个任务被划分为可以并行化的子任务，这些子任务被发送到不同的计算机，有时甚至跨越长距离。在这种情况下，通信时间可能会非常重要。只有当处理子任务的时间相对于通信时间较长时，使用计算机集群才有意义。
18.2 消息传递接口（MPI）
在多核计算机或分布式内存的计算机集群上编程需要特殊的技术。我们在这里描述了消息传递以及 MPI 标准化的相关工具。这些工具在不同的编程语言中相似，例如 C、C++和 FORTRAN，并通过mpi4py模块在 Python 中实现。
18.2.1 前提条件
你需要先通过在终端窗口执行以下命令来安装此模块：
conda install mpi4py
可以通过在你的 Python 脚本中添加以下一行来导入该模块：
import mpi4py as mpi
并行化代码的执行是通过终端使用命令mpiexec完成的。假设你的代码存储在文件script.py中，在一台具有四核 CPU 的计算机上执行此代码，可以在终端窗口通过运行以下命令来实现：
mpiexec -n 4 python script.py
或者，为了在一个包含两台计算机的集群上执行相同的脚本，可以在终端窗口运行以下命令：
mpiexec --hostfile=hosts.txt python script.py
你需要提供一个文件hosts.txt，其中包含你想绑定到集群的计算机的名称或 IP 地址，以及它们的核心数：
# Content of hosts.txt
192.168.1.25 :4 # master computer with 4 cores
192.168.1.101:2 # worker computer with 2 cores
Python 脚本（这里是script.py）必须被复制到集群中的所有计算机上。
18.3 将任务分配到不同的核心
当在多核计算机上执行时，我们可以认为mpiexec将给定的 Python 脚本复制到相应数量的核心上，并运行每个副本。例如，考虑一下包含命令print("Hello it's me")的单行脚本print_me.py，当通过mpiexec -n 4 print_me.py执行时，它会在屏幕上显示相同的消息四次，每次来自不同的核心。
为了能够在不同的核心上执行不同的任务，我们必须能够在脚本中区分这些核心。
为此，我们创建了一个所谓的通信实例，它组织了世界之间的通信，即输入输出单元，如屏幕、键盘或文件，与各个核心之间的通信。此外，每个核心都会被分配一个标识编号，称为 rank：
from mpi4py import MPI
comm=MPI.COMM_WORLD # making a communicator instance
rank=comm.Get_rank() # querrying for the numeric identifyer of the core
size=comm.Get_size() # the total number of cores assigned
通信器属性 size 指的是在mpiexec语句中指定的进程总数。
现在我们可以为每个核心分配一个独立的计算任务，就像在下一个脚本中那样，我们可以称之为basicoperations.py：
from mpi4py import MPI
comm=MPI.COMM_WORLD # making a communicator instance
rank=comm.Get_rank() # querrying for the numeric identifyer of the core
size=comm.Get_size() # the total number of cores assigned
a=15
b=2
if rank==0:
print(f'Core {rank} computes {a}+{b}={a+b}')
if rank==1:
print(f'Core {rank} computes {a}*{b}={a*b}')
if rank==2:
print(f'Core {rank} computes {a}**{b}={a**b}')
这个脚本可以通过在终端输入以下命令来执行：
mpiexec -n 3 python basicoperations.py
我们得到三个消息：
Core 0 computes 15+2=17
Core 2 computes 15**2=225
Core 1 computes 15*2=3
所有三个进程都有各自的任务，并且是并行执行的。显然，将结果打印到屏幕上是一个瓶颈，因为屏幕是所有三个进程共享的。
在下一节中，我们将看到进程间是如何进行通信的。
18.3.1 进程间信息交换
进程间有不同的发送和接收信息的方法：
-
点对点通信
-
单对多和多对单
-
多对多
在本节中，我们将介绍点对点、单对多和多对单通信。
向邻居讲话并让信息沿街道传递，这是前面列出的第一种通信类型的日常生活示例，而第二种可以通过新闻广播来说明，一人讲话并广播给一大群听众。单对多和多对单通信
https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/13f24ef2-b180-4b43-b9b3-8498ba57e23a.png https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/1aea42da-4134-4564-beab-a2904d20fdf9.png
图 18.2：点对点通信和单对多通信
在接下来的子节中，我们将研究这些不同的通信类型在计算上下文中的应用。
18.3.2 点对点通信
点对点通信将信息流从一个进程引导到指定的接收进程。我们首先通过考虑乒乓情况和电话链情况来描述方法和特点，并解释阻塞的概念。
点对点通信应用于科学计算，例如在分布域上的随机游走或粒子追踪应用，这些域被划分为多个子域，每个子域对应一个可以并行执行的进程数。
在这个乒乓示例中，我们假设有两个处理器相互发送一个整数，并将其值增加一。
我们从创建一个通信对象并检查是否有两个可用的进程开始：
from mpi4py import MPI
comm=MPI.COMM_WORLD # making a communicator instance
rank=comm.Get_rank() # querying for the numeric identifier of the core
size=comm.Get_size() # the total number of cores assigned
if not (size==2):
raise Exception(f"This examples requires two processes. \
{size} processes given.")
然后我们在两个进程之间来回发送信息：
count = 0
text=['Ping','Pong']
print(f"Rank {rank} activities:\n==================")
while count < 5:
if rank == count%2:
print(f"In round {count}: Rank {rank} says {text[count%2]}""
"and sends the ball to rank {(rank+1)%2}")
count += 1
comm.send(count, dest=(rank+1)%2)
elif rank == (count+1)%2:
count = comm.recv(source=(rank+1)%2)
信息通过通信器的send方法发送。在这里，我们提供了要发送的信息以及目的地。通信器确保将目的地信息转换为硬件地址；可以是你计算机的一个 CPU 核心，或主机的一个 CPU 核心。
另一台机器通过通信方法comm.recv接收信息。它需要知道信息来自哪里。在后台，它通过释放数据通道上的信息缓冲区，告诉发送方信息已经被接收。发送方在继续操作之前，需要等待此信号。
两个语句if rank == count%2和elif rank == (count+1)%2确保处理器交替进行发送和接收任务。
这是我们保存为pingpong.py文件并使用以下命令执行的短脚本输出：
mpiexec -n 2 python pingpong.py
在终端中，这会生成以下输出：
Rank 0 activities:
==================
In round 0: Rank 0 says Ping and sends the ball to rank 1
In round 2: Rank 0 says Ping and sends the ball to rank 1
In round 4: Rank 0 says Ping and sends the ball to rank 1
Rank 1 activities:
==================
In round 1: Rank 1 says Pong and sends the ball to rank 0
In round 3: Rank 1 says Pong and sends the ball to rank 0
可以发送或接收什么类型的数据？由于命令 send 和 recv 以二进制形式传递数据，因此它们首先会将数据进行 pickle（参见 第 14.3 节：Pickling）。大多数 Python 对象都可以被 pickled，但例如 lambda 函数不能。也可以 pickle 缓冲数据，例如 NumPy 数组，但直接发送缓冲数据更高效，正如我们将在下一小节中看到的那样。
请注意，在进程之间发送和接收函数可能有其原因。由于方法 send 和 recv 仅传递对函数的引用，因此这些引用必须在发送和接收的处理器上存在。因此，以下 Python 脚本会返回一个错误：
from mpi4py import MPI
comm=MPI.COMM_WORLD # making a communicator instance
rank=comm.Get_rank() # querying for the numeric identifier of the core
size=comm.Get_size() # the total number of cores assigned
if rank==0:
def func():
return 'Function called'
comm.send(func, dest=1)
if rank==1:
f=comm.recv(source=0) # <<<<<< This line reports an error
print(f())One-to-all and all-to-one communication
由语句 recv 抛出的错误信息是 AttributeError: Can't get attribute 'func'。这是由于 f 引用了 func 函数，而该函数在秩为 1 的处理器上没有定义。正确的做法是为两个处理器都定义该函数：
from mpi4py import MPI
comm=MPI.COMM_WORLD # making a communicator instance
rank=comm.Get_rank() # querying for the numeric identifier of the core
size=comm.Get_size() # the total number of cores assigned
def func():
return 'Function called'
if rank==0:
comm.send(func, dest=1)
if rank==1:
f=comm.recv(source=0)
print(f())
18.3.3 发送 NumPy 数组
命令 send 和 recv 是高级命令。这意味着它们在幕后执行工作，节省了程序员的时间并避免了可能的错误。它们会在内部推导出数据类型和所需通信缓冲区数据量后分配内存。这是在较低层次上基于 C 结构完成的。
NumPy 数组是对象，它们本身利用了类似 C 缓冲区的对象，因此在发送和接收 NumPy 数组时，可以通过在底层通信对等方 Send 和 Recv 中使用它们来提高效率（注意大小写！）。
在以下示例中，我们从一个处理器发送数组到另一个处理器：
from mpi4py import MPI
comm=MPI.COMM_WORLD # making a communicator instance
rank=comm.Get_rank() # querying for the numeric identifier of the core
size=comm.Get_size() # the total number of cores assigned
import numpy as np
if rank==0:
A = np.arange(700)
comm.Send(A, dest=1)
if rank==1:
A = np.empty(700, dtype=int) # This is needed for memory allocation
# of the buffer on Processor 1
comm.Recv(A, source=0) # Note, the difference to recv in
# providing the data.
print(f'An array received with last element {A[-1]}')
需要注意的是，在两个处理器上，必须分配缓冲区的内存。在这里，通过在处理器 0 上创建一个包含数据的数组以及在处理器 1 上创建一个具有相同大小和数据类型但包含任意数据的数组来完成这项工作。
此外，我们可以看到命令 recv 在输出中的区别。命令 Recv 通过第一个参数返回缓冲区。这是可能的，因为 NumPy 数组是可变的。
18.3.4 阻塞和非阻塞通信
命令 send 和 recv 及其缓冲区对应的 Send 和 Recv 是所谓的阻塞命令。这意味着，当相应的发送缓冲区被释放时，命令 send 才算完成。释放的时机取决于多个因素，例如系统的特定通信架构和要传输的数据量。最终，命令 send 在相应的命令 recv 接收到所有信息后才被认为是已释放的。如果没有这样的命令 recv，它将永远等待。这就形成了死锁情况。
以下脚本演示了可能发生死锁的情况。两个进程同时发送。如果要传输的数据量太大，无法存储，命令 send 就会等待相应的 recv 来清空管道，但由于等待状态，recv 永远不会被调用。这就是死锁。
from mpi4py import MPI
comm=MPI.COMM_WORLD # making a communicator instance
rank=comm.Get_rank() # querrying for the numeric identifier of the core
size=comm.Get_size() # the total number of cores assigned
if rank==0:
msg=['Message from rank 0',list(range(101000))]
comm.send(msg, dest=1)
print(f'Process {rank} sent its message')
s=comm.recv(source=1)
print(f'I am rank {rank} and got a {s[0]} with a list of \
length {len(s[1])}')
if rank==1:
msg=['Message from rank 1',list(range(-101000,1))]
comm.send(msg,dest=0)
print(f'Process {rank} sent its message')
s=comm.recv(source=0)
print(f'I am rank {rank} and got a {s[0]} with a list of \
length {len(s[1])}')
注意，执行这段代码可能不会导致你的计算机死锁，因为要通信的数据量非常小。
在这种情况下，避免死锁的直接解决办法是交换命令 recv 和 send 在 一个 处理器上的执行顺序：
from mpi4py import MPI
comm=MPI.COMM_WORLD # making a communicator instance
rank=comm.Get_rank() # querrying for the numeric identifier of the core
size=comm.Get_size() # the total number of cores assigned
if rank==0:
msg=['Message from rank 0',list(range(101000))]
comm.send(msg, dest=1)
print(f'Process {rank} sent its message')
s=comm.recv(source=1)
print(f'I am rank {rank} and got a {s[0]} with a list of \
length {len(s[1])}')
if rank==1:
s=comm.recv(source=0)
print(f'I am rank {rank} and got a {s[0]} with a list of \
length {len(s[1])}')
msg=['Message from rank 1',list(range(-101000,1))]
comm.send(msg,dest=0)
print(f'Process {rank} sent its message')
print(f'I am rank {rank} and got a {s[0]} with a list of \
length {len(s[1])}')
18.3.5 一对多与多对一通信
当一个依赖于大量数据的复杂任务被分解为子任务时，数据也必须分成与相关子任务相关的部分，并且结果必须汇总并处理成最终结果。
以两个向量的标量积为例，考虑将其划分为！的子任务：https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/4093f1e5-0763-4c6d-9c86-3d791dbe381b.png
使用 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/83232989-a3d1-4bce-8ab1-8c511977a52d.png 所有子任务在初始数据的各个部分上执行相同的操作，结果必须汇总，并可能执行剩余的操作。
我们需要执行以下步骤：
-
创建向量
u和v -
将它们分成 m 个子向量，且每个子向量的元素个数平衡，即当
N能被m整除时，每个子向量包含 ![] 元素，否则一些子向量会包含更多元素。 -
将每个子向量传递给 “它的” 处理器
-
在每个处理器上执行子向量的标量积
-
收集所有结果
-
汇总结果
步骤 1、2 和 6 在一个处理器上运行，即所谓的 根 处理器。在以下示例代码中，我们选择排名为 0 的处理器来执行这些任务。步骤 3、4 和 5 在所有处理器上执行，包括根处理器。对于步骤 3 中的通信，mpi4py 提供了命令 scatter，而用于收集结果的命令是 gather。
准备通信数据
首先，我们来看看步骤 2。编写一个脚本，将一个向量分成 m 个平衡元素的部分，是一个不错的练习。这里有一个建议的脚本实现，当然还有很多其他的实现方式：
def split_array(vector, n_processors):
# splits an array into a number of subarrays
# vector one dimensional ndarray or a list
# n_processors integer, the number of subarrays to be formed
n=len(vector)
n_portions, rest = divmod(n,n_processors) # division with remainder
# get the amount of data per processor and distribute the res on
# the first processors so that the load is more or less equally
# distributed
# Construction of the indexes needed for the splitting
counts = [0]+ [n_portions + 1 \
if p < rest else n_portions for p in range(n_processors)]
counts=numpy.cumsum(counts)
start_end=zip(counts[:-1],counts[1:]) # a generator
slice_list=(slice(*sl) for sl in start_end) # a generator comprehension
return [vector[sl] for sl in slice_list] # a list of subarrays
由于本章是本书中的最后几章之一，我们已经看到很多可以用于这段代码的工具。我们使用了 NumPy 的累积和 cumsum。我们使用了生成器 zip，通过运算符 * 解包参数，以及生成器推导式。我们还默默地引入了数据类型 slice，它允许我们在最后一行以非常简洁的方式执行分割步骤。
命令——scatter 和 gather
现在我们已经准备好查看整个脚本，来解决我们的演示问题——标量积：
from mpi4py import MPI
import numpy as np
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
nprocessors = comm.Get_size()
import splitarray as spa
if rank == 0:
# Here we generate data for the example
n = 150
u = 0.1*np.arange(n)
v = - u
u_split = spa.split_array(u, nprocessors)
v_split = spa.split_array(v, nprocessors)
else:
# On all processor we need variables with these names,
# otherwise we would get an Exception "Variable not defined" in
# the scatter command below
u_split = None
v_split = None
# These commands run now on all processors
u_split = comm.scatter(u_split, root=0) # the data is portion wise
# distributed from root
v_split = comm.scatter(v_split, root=0)
# Each processor computes its part of the scalar product
partial_dot = u_split@v_split
# Each processor reports its result back to the root
partial_dot = comm.gather(partial_dot,root=0)
if rank==0:
# partial_dot is a list of all collected results
total_dot=np.sum(partial_dot)
print(f'The parallel scalar product of u and v'
f'on {nprocessors} processors is {total_dot}.\n'
f'The difference to the serial computation is \
{abs(total_dot-u@v)}')
如果此脚本存储在文件parallel_dot.py中，使用五个处理器执行的命令如下：
mexec -n 5 python parallel_dot.py
在这种情况下，结果如下：
The parallel scalar product of u and v on 5 processors is -11137.75.
The difference to the serial computation is 0.0
本示例演示了使用scatter将特定信息发送到每个处理器。要使用此命令，根处理器必须提供一个包含与可用处理器数量相同元素的列表。每个元素包含要传送到某个处理器的数据，包括根处理器本身。
反向过程是gather。当所有处理器完成此命令时，根处理器将得到一个包含与可用处理器数量相同元素的列表，每个元素包含其对应处理器的结果数据。
在最后一步，根处理器再次独自工作，后处理此结果列表。上面的示例将所有列表元素求和并显示结果。
并行编程的艺术在于避免瓶颈。理想情况下，所有处理器都应保持忙碌，并且应该同时开始和结束。这就是为什么我们前面描述的脚本splitarray将工作负载大致平均分配给处理器的原因。此外，代码应以这样的方式组织，即根处理器独自工作的开始和结束阶段，相对于所有处理器同时执行的计算密集型部分来说是很短的。
最终数据归约操作——命令 reduce
并行标量积示例是许多其他任务的典型示例，展示了结果处理方式：来自所有处理器的数据量在最后一步被归约为一个单一的数字。在这里，根处理器将所有处理器的部分结果相加。命令reduce可以有效地用于此任务。我们通过让reduce在一步中完成聚集和求和来修改之前的代码。以下是修改后的前几行代码：
......... modification of the script above .....
# Each processor reports its result back to the root
# and these results are summed up
total_dot = comm.reduce(partial_dot, op=MPI.SUM, root=0)
if rank==0:
print(f'The parallel scalar product of u and v'
f' on {nprocessors} processors is {total_dot}.\n'
f'The difference to the serial computation \
is {abs(total_dot-u@v)}')
其他常用的归约操作有：
-
MPI.MAX或MPI.MIN：部分结果的最大值或最小值 -
MPI.MAXLOC或MPI.MINLOC：部分结果的最大值位置或最小值位置 -
MPI.PROD：部分结果的乘积 -
MPI.LAND或MPI.LOR：部分结果的逻辑与/逻辑或
向所有处理器发送相同的消息
另一个集体命令是广播命令bcast。与scatter不同，它用于将相同的数据发送到所有处理器。它的调用方式与scatter类似：
data = comm.bcast(data, root=0)
但是，发送的是总数据，而不是分割数据的列表。同样，根处理器可以是任何处理器。它是准备广播数据的处理器。
缓冲数据
类似地，mpi4py为类似 NumPy 数组的缓冲数据提供了相应的集体命令，通过大写命令来实现：scatter/Scatter、gather/Gather、reduce/Reduce、bcast/Bcast。
18.4 总结
在本章中，我们了解了如何在不同的处理器上并行执行相同脚本的副本。消息传递允许这些不同进程之间进行通信。我们看到了点对点通信，以及两种不同的分布式集体通信类型：一对多和多对一。本章中展示的命令是由 Python 模块mpi4py提供的，这是一个 Python 封装器，用于实现 C 语言中的 MPI 标准。
经过本章的学习后，你现在能够编写自己的并行编程脚本，并且你会发现我们这里只描述了最基本的命令和概念。进程分组和信息标记只是我们遗漏的两个概念。许多这些概念对于特殊和具有挑战性的应用非常重要，但它们对于本介绍来说过于具体。
综合示例
在本章中，我们将提供一些综合性和较长的示例，并简要介绍理论背景以及示例的完整实现。在这里，我们希望向您展示本书中定义的概念如何在实践中应用。
本章将涵盖以下主题：
-
多项式
-
多项式类
-
谱聚类
-
求解初值问题
第二十章：19.1 多项式
首先，我们将通过设计一个多项式类来展示迄今为止所介绍的 Python 构造的强大功能。
注意，这个类在概念上与类 numpy.poly1d 不同。
我们将提供一些理论背景，这将引导我们列出需求，然后给出代码并附带一些注释。
19.1.1 理论背景
或者，可以将多项式写成：
使用基函数：
有无穷多种表示方法，但我们这里只限制于这三种典型表示。
多项式可以通过插值条件确定：
给定不同的值 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/df8177e3-c79a-47ce-bc35-4799dd8db630.png 和任意值 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/ac1fbb98-6657-4d5e-9e9f-8a7a253adf9d.png 作为输入。在拉格朗日公式中，插值多项式是直接可用的，因为其系数即为插值数据。牛顿表示法中的插值多项式系数可以通过递推公式获得，称为分差公式：
和
单项式表示法中插值多项式的系数通过解线性系统获得：
一个具有给定多项式https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/2baead6e-97d1-4d9e-8b27-d18720ef5356.png（或其倍数）作为特征多项式的矩阵被称为伴随矩阵。伴随矩阵的特征值即为多项式的零点（根）。通过先建立其伴随矩阵，再使用scipy.linalg.eig计算特征值，可以构建一个计算https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/5b023d20-84ee-4e91-b5e8-f6620cfb7579.png零点的算法。牛顿表示下的多项式的伴随矩阵如下所示：
19.1.2 任务
我们现在可以制定一些编程任务：
- 编写一个名为
PolyNomial的类，具有points、degree、coeff和basis属性，其中：
-
为类提供一个方法，用于在给定点上评估多项式。
-
为类提供一个名为
plot的方法，用于在给定区间内绘制多项式。 -
编写一个名为
__add__的方法，返回两个多项式的和。需要注意的是，只有在单项式情况下，才能通过简单地将系数相加来计算和。 -
编写一个方法，计算表示为单项式形式的多项式的系数。
-
编写一个方法，计算多项式的伴随矩阵。
-
编写一个方法，通过计算伴随矩阵的特征值来计算多项式的零点。
-
编写一个方法，计算给定多项式的https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/9ca472ad-554e-4ea9-b316-765cf23dbeac.png^(th) 导数。
-
编写一个方法，检查两个多项式是否相等。可以通过比较所有系数来检查相等性（零的首项系数不应影响结果）。
19.1.3 多项式类
现在我们基于单项式形式设计一个多项式基类。多项式可以通过给定其相对于单项式基的系数，或通过给定插值点列表来初始化，如下所示：
import scipy.linalg as sl
import matplotlib.pyplot as mp
class PolyNomial:
base='monomial'
def __init__(self,**args):
if 'points' in args:
self.points = array(args['points'])
self.xi = self.points[:,0]
self.coeff = self.point_2_coeff()
self.degree = len(self.coeff)-1
elif 'coeff' in args:
self.coeff = array(args['coeff'])
self.degree = len(self.coeff)-1
self.points = self.coeff_2_point()
else:
self.points = array([[0,0]])
self.xi = array([1.])
self.coeff = self.point_2_coeff()
self.degree = 0
新类的__init__方法使用了构造**args，如第 7.2.5 节中讨论的：可变参数个数。如果没有给定参数，则默认假定为零多项式。如果多项式由插值点给出，计算系数的方法是通过求解范德蒙矩阵系统，方法如下：
def point_2_coeff(self):
return sl.solve(vander(self.x),self.y)
从![] 给定的系数，![] 插值点由以下方式构造：
def coeff_2_point(self):
points = [[x,self(x)] for x in linspace(0,1,self.degree+1)]
return array(points)
命令self(x)进行多项式评估，这通过提供__call__方法来实现：
def __call__(self,x):
return polyval(self.coeff,x)
（参见第 8.1.5 节中的__call__方法示例：特殊方法。）这里，该方法使用了 NumPy 命令polyval。接下来的步骤，我们只需为了方便添加两个方法，并用property装饰器装饰它们（参见第 7.8 节：作为装饰器的函数*)*：
@property
def x(self):
return self.points[:,0]
@property
def y(self):
return self.points[:,1]
让我们解释一下这里发生了什么。我们定义了一个方法来提取定义多项式时使用的数据的 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/51485369-65e9-4da0-811d-1947fbb24f27.png值。类似地，也定义了一个提取数据 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3a9672e4-ff7e-43b3-8049-7fc79880852c.png值的方法。通过property装饰器，调用该方法的结果呈现得就像它是多项式的一个属性一样。这里有两种编码选择：
- 我们使用方法调用：
def x(self):
return self.points[:,0]
这通过调用：p.x()提供对 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/4b12e4f2-67b2-4d6f-b995-0911ae5ca323.png值的访问。
- 我们使用
property装饰器。它允许我们通过以下语句轻松访问 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/35b8b3fb-62ea-4a59-9044-bcb3cffd718c.png值：p.x。我们在这里选择第二种变体。
定义__repr__方法始终是一个好习惯（参见第 8.1.5 节：特殊方法）。至少对于快速检查结果，这个方法非常有用：
def __repr__(self):
txt = f'Polynomial of degree {self.degree} \n'
txt += f'with coefficients {self.coeff} \n in {self.base} basis.'
return txt
现在，我们提供一个用于绘制多项式的方法，如下所示：
margin = .05
plotres = 500
def plot(self,ab=None,plotinterp=True):
if ab is None: # guess a and b
x = self.x
a, b = x.min(), x.max()
h = b-a
a -= self.margin*h
b += self.margin*h
else:
a,b = ab
x = linspace(a,b,self.plotres)
y = vectorize(self.__call__)(x)
mp.plot(x,y)
mp.xlabel('$x$')
mp.ylabel('$p(x)$')
if plotinterp:
mp.plot(self.x, self.y, 'ro')
请注意使用vectorize命令（参见第 4.8 节：作用于数组的函数）。__call__方法是特定于单项式表示的，如果多项式以其他方式表示，则必须进行更改。这对于计算多项式的伴随矩阵也是如此：
def companion(self):
companion = eye(self.degree, k=-1)
companion[0,:] -= self.coeff[1:]/self.coeff[0]
return companion
一旦伴随矩阵可用，给定多项式的零点就是它的特征值：
def zeros(self):
companion = self.companion()
return sl.eigvals(companion)
为此，必须首先导入模块scipy.linalg，并命名为sl。
19.1.4 多项式类的使用示例
让我们给出一些使用示例。
首先，我们从给定的插值点创建一个多项式实例：
p = PolyNomial(points=[(1,0),(2,3),(3,8)])
相对于单项式基，多项式的系数作为p的一个属性提供：
p.coeff # returns array([ 1., 0., -1.]) (rounded)
这对应于多项式 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/7d283009-8ea6-4ff2-b454-6ef8dba0c847.png 。通过p.plot((-3.5,3.5))获得的多项式默认图形如下所示（图 19.1）：
图 19.1：多项式绘图方法的结果
最后，我们计算多项式的零点，在此案例中是两个实数：
pz = p.zeros() # returns array([-1.+0.j, 1.+0.j])
结果可以通过在这些点处评估多项式来验证：
p(pz) # returns array([0.+0.j, 0.+0.j])
19.1.5 牛顿多项式
类 NewtonPolyNomial 定义了一个基于牛顿基的多项式。我们通过使用命令 super 让它继承自多项式基类的一些常用方法，例如 polynomial.plot、polynomial.zeros，甚至 __init__ 方法的部分内容（参见 第 8.5 节：子类与继承）：
class NewtonPolynomial(PolyNomial):
base = 'Newton'
def __init__(self,**args):
if 'coeff' in args:
try:
self.xi = array(args['xi'])
except KeyError:
raise ValueError('Coefficients need to be given'
'together with abscissae values xi')
super(NewtonPolynomial, self).__init__(**args)
一旦给定插值点，系数的计算就通过以下方式进行：
def point_2_coeff(self):
return array(list(self.divdiff()))
我们使用了分割差分来计算多项式的牛顿表示，这里编程为生成器（参见 第 9.3.1 节：生成器 和 第 9.4 节：列表填充模式）：
def divdiff(self):
xi = self.xi
row = self.y
yield row[0]
for level in range(1,len(xi)):
row = (row[1:] - row[:-1])/(xi[level:] - xi[:-level])
if allclose(row,0): # check: elements of row nearly zero
self.degree = level-1
break
yield row[0]
让我们简单检查一下这个是如何工作的：
# here we define the interpolation data: (x,y) pairs
pts = array([[0.,0],[.5,1],[1.,0],[2,0.]])
pN = NewtonPolynomial(points=pts) # this creates an instance of the
# polynomial class
pN.coeff # returns the coefficients array([0\. , 2\. , -4\. , 2.66666667])
print(pN)
函数 print 执行基类的 __repr__ 方法并返回以下文本：
Polynomial of degree 3
with coefficients [ 0\. 2\. -4\. 2.66666667]
in Newton basis.
多项式评估与基类的相应方法不同。方法 NewtonPolyNomial.__call__ 需要重写 Polynomial.__call__：
def __call__(self,x):
# first compute the sequence 1, (x-x_1), (x-x_1)(x-x_2),...
nps = hstack([1., cumprod(x-self.xi[:self.degree])])
return self.coeff@nps
最后，我们给出了伴随矩阵的代码，该代码重写了父类的相应方法，如下所示：
def companion(self):
degree = self.degree
companion = eye(degree, k=-1)
diagonal = identity(degree,dtype=bool)
companion[diagonal] = self.x[:degree]
companion[:,-1] -= self.coeff[:degree]/self.coeff[degree]
return companion
请注意布尔数组的使用。练习将进一步在此基础上构建。
19.2 谱聚类
特征向量的一个有趣应用是用于数据聚类。通过使用由距离矩阵导出的矩阵的特征向量，可以将未标记的数据分成不同的组。谱聚类方法因使用该矩阵的谱而得名。对于 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/e27f1d49-c801-4539-8009-9cae5853b3cf.png 元素（例如，数据点之间的成对距离）的距离矩阵是一个 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/a36e8b17-7e0b-488d-96b3-29c79c5257f2.png 对称矩阵。给定这样的 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/6b3e95b7-0207-425e-b217-6a476a631685.png 距离矩阵 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/3972f882-789e-45e6-afcc-d1663436700b.png，具有距离值 ![]，我们可以按以下方式创建数据点的拉普拉斯矩阵：
这里，https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/d9020f81-cd8b-4d22-bb23-c67c1326984c.png 是单位矩阵，https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/04f5ec10-48e3-4051-822a-b701c6660c95.png 是包含 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/a68ab07c-1d27-4a51-ba39-01c32f498aef.png 行和列和的对角矩阵：
数据聚类是通过 L 的特征向量得到的。在只有两类数据点的最简单情况下，第一个特征向量（即对应最大特征值的那个）通常足以将数据分开。
这是一个简单的二类聚类示例。以下代码创建了一些二维数据点，并基于拉普拉斯矩阵的第一个特征向量对其进行聚类：
import scipy.linalg as sl
# create some data points
n = 100
x1 = 1.2 * random.randn(n, 2)
x2 = 0.8 * random.randn(n, 2) + tile([7, 0],(n, 1))
x = vstack((x1, x2))
# pairwise distance matrix
M = array([[ sqrt(sum((x[i] - x[j])**2))
for i in range(2*n)]
for j in range(2 * n)])
# create the Laplacian matrix
D = diag(1 / sqrt( M.sum(axis = 0) ))
L = identity(2 * n) - dot(D, dot(M, D))
# compute eigenvectors of L
S, V = sl.eig(L)
# As L is symmetric the imaginary parts
# in the eigenvalues are only due to negligible numerical errors S=S.real
V=V.real
对应最大特征值的特征向量给出了分组（例如，通过在 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/e6b3a21e-af03-46d9-80e9-3822e2be13d2.png 处进行阈值化）并可以通过以下方式展示：
largest=abs(S).argmax()
plot(V[:,largest])
以下图 (图 19.2) 显示了简单两类数据集的谱聚类结果：
图 19.2：简单两类聚类的结果
对于更难的数据集和更多的类别，通常会采用与最大特征值对应的 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/67e01793-7925-419d-9a24-2f5622e534c5.png 特征向量，然后使用其他方法对数据进行聚类，但使用的是特征向量而不是原始数据点。常见的选择是 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/2e934e12-050c-43a9-b47d-31e0f3116490.png-均值聚类算法，这是下一个示例的主题。
特征向量作为输入用于 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/8ec8c440-a221-41d9-a7e7-5e81d4aa3d23.png-均值聚类，如下所示：
import scipy.linalg as sl
import scipy.cluster.vq as sc
# simple 4 class data
x = random.rand(1000,2)
ndx = ((x[:,0] < 0.4) | (x[:,0] > 0.6)) & \
((x[:,1] < 0.4) | (x[:,1] > 0.6))
x = x[ndx]
n = x.shape[0]
# pairwise distance matrix
M = array([[ sqrt(sum((x[i]-x[j])**2)) for i in range(n) ]
for j in range(n)])
# create the Laplacian matrix
D = diag(1 / sqrt( M.sum(axis=0) ))
L = identity(n) - dot(D, dot(M, D))
# compute eigenvectors of L
_,_,V = sl.svd(L)
k = 4
# take k first eigenvectors
eigv = V[:k,:].T
# k-means
centroids,dist = sc.kmeans(eigv,k)
clust_id = sc.vq(eigv,centroids)[0]
请注意，我们在此使用奇异值分解 sl.svd 计算了特征向量。由于 L 是对称的，结果与使用 sl.eig 得到的结果相同，但 svd 已经按照特征值的顺序给出了特征向量。我们还使用了临时变量。svd 返回一个包含三个数组的列表：左奇异向量 U、右奇异向量 V 和奇异值 S，如下所示：
U, S, V = sl.svd(L)
由于我们在这里不需要 U 和 S，因此可以在解包 svd 返回值时丢弃它们：
_, _, V = sl.svd(L)
结果可以通过以下方式绘制：
for i in range(k):
ndx = where(clust_id == i)[0]
plot(x[ndx, 0], x[ndx, 1],'o')
axis('equal')
以下图显示了简单 多类数据集 的谱聚类结果：
图 19.3：简单四类数据集的谱聚类示例
19.3 求解初值问题
在本节中，我们将考虑数值求解给定初值的常微分方程组的数学任务：
该问题的解是一个函数 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/b2d66ae4-6c74-4bba-996c-d8e96dfe0de4.png。数值方法在离散的通信点 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/86022693-bb97-4372-8e87-d5da6318504d.png 计算近似值，https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/eca4b607-0793-4c58-b3e8-ea1fc795b35d.png，位于感兴趣的区间 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/623b90f9-82a7-4c3d-8b57-3e05a6e609bb.png 内。我们将描述问题的数据收集到一个类中，如下所示：
class IV_Problem:
"""
Initial value problem (IVP) class
"""
def __init__(self, rhs, y0, interval, name='IVP'):
"""
rhs 'right hand side' function of the ordinary differential
equation f(t,y)
y0 array with initial values
interval start and end value of the interval of independent
variables often initial and end time
name descriptive name of the problem
"""
self.rhs = rhs
self.y0 = y0
self.t0, self.tend = interval
self.name = name
微分方程：
描述了一个数学摆；* https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/d2dd5a5d-e940-4681-8d88-8805e1ad95bb.png，初始角速度为零。
摆锤问题成为问题类的一个实例，如下所示：
def rhs(t,y):
g = 9.81
l = 1.
yprime = array([y[1], g / l * sin(y[0])])
return yprime
pendulum = IV_Problem(rhs, array([pi / 2, 0.]), [0., 10.] ,
'mathem. pendulum')
对当前问题可能会有不同的看法，这会导致不同的类设计。例如，有人可能希望将自变量的区间视为解过程的一部分，而不是问题定义的一部分。考虑初值时也是如此。我们在这里将其视为数学问题的一部分，而其他作者可能希望将初值作为解过程的一部分，从而允许初值的变化。
解决过程被建模为另一个类：
class IVPsolver:
"""
IVP solver class for explicit one-step discretization methods
with constant step size
"""
def __init__(self, problem, discretization, stepsize):
self.problem = problem
self.discretization = discretization
self.stepsize = stepsize
def one_stepper(self):
yield self.problem.t0, self.problem.y0
ys = self.problem.y0
ts = self.problem.t0
while ts <= self.problem.tend:
ts, ys = self.discretization(self.problem.rhs, ts, ys,
self.stepsize)
yield ts, ys
def solve(self):
return list(self.one_stepper())
我们接下来首先定义两个离散化方案：
- 显式欧拉方法：
def expliciteuler(rhs, ts, ys, h):
return ts + h, ys + h * rhs(ts, ys)
- 经典龙格-库塔四阶法（RK4）：
def rungekutta4(rhs, ts, ys, h):
k1 = h * rhs(ts, ys)
k2 = h * rhs(ts + h/2., ys + k1/2.)
k3 = h * rhs(ts + h/2., ys + k2/2.)
k4 = h * rhs(ts + h, ys + k3)
return ts + h, ys + (k1 + 2*k2 + 2*k3 + k4)/6.
使用这些工具，我们可以创建实例以获得相应的摆锤常微分方程的离散化版本：
pendulum_Euler = IVPsolver(pendulum, expliciteuler, 0.001)
pendulum_RK4 = IVPsolver(pendulum, rungekutta4, 0.001)
我们可以解决这两个离散模型并绘制解与角度差：
sol_Euler = pendulum_Euler.solve()
sol_RK4 = pendulum_RK4.solve()
tEuler, yEuler = zip(*sol_Euler)
tRK4, yRK4 = zip(*sol_RK4)
subplot(1,2,1), plot(tEuler,yEuler),\
title('Pendulum result with Explicit Euler'),\
xlabel('Time'), ylabel('Angle and angular velocity')
subplot(1,2,2), plot(tRK4,abs(array(yRK4)-array(yEuler))),\
title('Difference between both methods'),\
xlabel('Time'), ylabel('Angle and angular velocity')
图 19.4：使用显式欧拉方法的摆锤模拟，并与更精确的龙格-库塔四阶法的结果进行比较
讨论替代的类设计是值得的。什么应该放在独立的类中，什么应该归并到同一个类中？
-
我们严格将数学问题与数值方法分开。初值应该放在哪里？它们应该是问题的一部分还是解算器的一部分？或者它们应该作为解算器实例的求解方法的输入参数？甚至可以设计程序，允许多种可能性。选择使用这些替代方案之一取决于未来程序的使用。像参数识别中的各种初值循环，留初值作为求解方法的输入参数会更方便。另一方面，模拟具有相同初值的不同模型变体时，将初值与问题绑定起来会更为合理。
-
为简化起见，我们仅展示了使用常定步长的求解器。
IVPsolver类的设计是否适用于未来扩展自适应方法，其中给定的是容差而非步长？ -
我们之前建议使用生成器构造来实现步进机制。自适应方法需要时不时地拒绝某些步长。这种需求是否与
IVPsolver.onestepper中步进机制的设计冲突？ -
我们鼓励你检查两个 SciPy 工具的设计，这些工具用于解决初值问题，分别是
scipy.integrate.ode和scipy.integrate.odeint。
19.4 小结
本书中大部分内容已经整合到本章的三个较长示例中。这些示例模拟了代码开发并提供了原型，鼓励你根据自己的想法进行修改和对比。
你会发现，科学计算中的代码因其与数学定义的算法有着紧密关系，往往具有其独特风格，并且通常明智的做法是保持代码与公式之间的关系可见。Python 提供了实现这一点的技巧，正如你所看到的。
19.5 练习
例 1： 实现一个方法 __add__，通过将两个给定的多项式 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/724805b9-bea3-4a02-b1d4-5d81e8bd4634.png 和 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/e174fe97-c3d6-4bfb-a93a-1025e6028d8a.png 相加，构造一个新的多项式 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/77fa6522-5de3-4e15-b0df-d9a450605771.png。在单项式形式中，多项式的相加只是将系数相加，而在牛顿形式中，系数依赖于插值点的横坐标 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/171118dd-dce4-49d6-8f04-7b4855f2b4f2.png。在相加两个多项式的系数之前，多项式 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/8bdd5965-3db4-42e5-9eae-f37213bda9ef.png 必须获得新的插值点，这些插值点的横坐标 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/40eb87a9-fc00-48f8-a883-b285e80c460c.png 必须与 https://github.com/OpenDocCN/freelearn-ds-pt3-zh/raw/master/docs/sci-comp-py-2e/img/296d4fcd-75bd-4819-8c45-3e084799040e.png 的横坐标一致，并且必须提供方法 __changepoints__ 来实现这一点。该方法应更改插值点，并返回一组新的系数。
例 2： 编写转换方法，将一个多项式从牛顿形式转换为单项式形式，反之亦然。
例 3： 编写一个名为 add_point 的方法，该方法接受一个多项式 q 和一个元组 ![] 作为参数，并返回一个新的多项式，该多项式插值 self.points 和 ![]。
例 4： 编写一个名为 LagrangePolynomial 的类，该类实现拉格朗日形式的多项式，并尽可能多地继承多项式基类。
例 5： 编写多项式类的测试。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
