Python3 科学计算教程（三）
九、重复
在本章中，我们将介绍使用循环和迭代器的迭代。我们将展示如何将它用于列表和生成器的示例。迭代是计算机有用的基本操作之一。传统上，迭代是通过for循环来实现的。for循环是指令块重复一定次数。在循环内部，可以访问一个循环变量，其中存储了迭代号。
Python 的习惯用法略有不同。Python 中的for循环主要是为了穷尽一个列表，也就是枚举一个列表的元素。如果使用包含第一个 n 个整数的列表，效果类似于刚才描述的重复效果。
一个for循环一次只需要列表中的一个元素。因此，希望对能够按需创建这些元素的对象使用for循环，一次一个。这就是迭代器在 Python 中实现的功能。
for 语句
for语句的主要目的是遍历列表:
for s in ['a', 'b', 'c']:
print(s), # a b c
在本例中，循环变量 s 被连续分配给列表中的一个元素。请注意，循环变量在循环结束后可用。这有时可能有用；例如，参考章节中的例子控制回路内的流量。
for循环最常见的用途之一是使用功能range将给定任务重复定义的次数(参见第 1 章、入门的列表部分)。
for iteration in range(n): # repeat the following code n times
...
如果循环的目的是遍历列表，许多语言(包括 Python)提供以下模式:
for k in range(...):
...
element = my_list[k]
如果该代码的目的是浏览列表my_list，前面的代码不会使它非常清楚。因此，更好的表达方式如下:
for element in my_list:
...
现在乍看之下很清楚，前面的代码通过了my_list列表。请注意，如果您真的需要索引变量 k ，您可以用以下代码替换前面的代码:
for k, element in enumerate(my_list):
...
这段代码的目的是在保持索引变量 k 可用的同时通过my_list。类似的数组结构是命令ndenumerate。
控制回路内的流量
有时需要跳出循环，或者直接进入下一个循环迭代。这两个操作由break和continue命令执行。顾名思义，break关键字打破了这个循环。循环中断时会出现两种情况:
- 循环被完全执行。
- 该循环在完全执行之前被留下(
break)。
对于第一种情况，可以在else块中定义特殊动作，如果遍历了整个列表，则执行该块。如果for循环的目的是寻找某样东西并停止，这通常是有用的。例如，在列表中搜索一个满足某个属性的元素。如果没有找到这样的元素，则执行else块。
这是科学计算中的一个常见用法。通常，我们使用的迭代算法不能保证成功。在这种情况下，最好使用(大的)有限循环，这样程序就不会陷入无限循环。for / else构造允许这样的实现:
maxIteration = 10000
for iteration in range(maxIteration):
residual = compute() # some computation
if residual < tolerance:
break
else: # only executed if the for loop is not broken
raise Exception("The algorithm did not converge")
print("The algorithm converged in {} steps".format(iteration+1))
迭代器
for循环主要用于遍历列表，但它一次只选取列表中的一个元素。特别是，不需要将整个列表存储在内存中，循环就能正常工作。允许for循环在没有列表的情况下工作的机制是迭代器。
可迭代对象产生对象(传递到for循环)。这样的物体obj可以在for回路中使用，如下所示:
for element in obj:
...
迭代器的概念概括了列表的概念。列表给出了可迭代对象的最简单的例子。生成的对象只是存储在列表中的对象:
L = ['A', 'B', 'C']
for element in L:
print(element)
可迭代对象不需要产生现有对象。相反，这些物体可以在飞行中产生。
典型的可迭代对象是函数range返回的对象。这个函数的工作方式就好像它会生成一个整数列表，但是相反，连续的整数是在需要时动态生成的:
for iteration in range(100000000):
# Note: the 100000000 integers are not created at once
if iteration > 10:
break
如果真的需要一个 0 到 100，000，000 之间的所有整数的列表，那么它必须是显式形成的:
l=list(range(100000000))
发电机
您可以使用yield关键字创建自己的迭代器。例如，小于 n 的奇数发生器可以定义为:
def odd_numbers(n):
"generator for odd numbers less than n"
for k in range(n):
if k % 2 == 1:
yield k
那么您可以按如下方式使用它:
g = odd_numbers(10)
for k in g:
... # do something with k
或者甚至像这样:
for k in odd_numbers(10):
... # do something with k
迭代器是一次性的
迭代器的一个显著特点是它们只能使用一次。为了再次使用迭代器，您必须创建一个新的迭代器对象。请注意，可迭代对象可以根据需要多次创建新的迭代器。让我们来研究一下列表的情况:
L = ['a', 'b', 'c']
iterator = iter(L)
list(iterator) # ['a', 'b', 'c']
list(iterator) # [] empty list, because the iterator is exhausted
new_iterator = iter(L) # new iterator, ready to be used
list(new_iterator) # ['a', 'b', 'c']
每次调用生成器对象时，都会创建一个新的迭代器。因此，当迭代器用完时，必须再次调用生成器来获得新的迭代器:
g = odd_numbers(10)
for k in g:
... # do something with k
# now the iterator is exhausted:
for k in g: # nothing will happen!!
...
# to loop through it again, create a new one:
g = odd_numbers(10)
for k in g:.
...
迭代器工具
这里有几个迭代器工具，它们经常会派上用场:
enumerate用于枚举另一个迭代器。它产生一个新的迭代器，产生对(迭代，元素)，其中iteration存储迭代的索引:
A = ['a', 'b', 'c']
for iteration, x in enumerate(A):
print(iteration, x)
# result: (0, 'a') (1, 'b') (2, 'c')
reversed通过向后遍历列表，从列表中创建一个迭代器。请注意，这不同于创建反向列表:
A = [0, 1, 2]
for elt in reversed(A):,
print(elt)
# result: 2 1 0
itertools.count可能是整数的无限迭代器:
for iteration in itertools.count():
if iteration > 100:
break # without this, the loop goes on forever
print("integer {}".format(iteration))
# prints the 100 first integer
intertools.islice使用熟悉的slicing语法截断迭代器；参考第三章、集装箱类型。一个应用正在从无限迭代器创建有限迭代器:
from itertools import count, islice
for iteration in islice(count(), 10):
# same effect as range(10)
...
例如，让我们通过将islice与一个无限发生器相结合来找到一些奇数。首先，我们修改奇数的生成器，使其成为无限生成器:
def odd_numbers():
k=-1
while True:
k+=1
if k%2==1:
yield k
然后，我们用它和islice一起得到一些奇数的列表:
list(itertools.islice(odd_numbers(),10,30,8)) # returns [21, 37, 53]
递归序列的生成器
假设一个序列由一个归纳公式给出。例如，考虑由递归公式定义的斐波那契数列:*un= un*T6】-1+uT11】n-2。
该序列取决于两个初始值，即 u 0 和 u 1 ，尽管对于标准斐波那契序列，这些数字分别取为 0 和 1。编写这种序列生成程序的一种巧妙方法是使用生成器，如下所示:
def fibonacci(u0, u1):
"""
Infinite generator of the Fibonacci sequence.
"""
yield u0
yield u1
while True:
u0, u1 = u1, u0+u1
yield u1
例如，可以这样使用:
# sequence of the 100 first Fibonacci numbers:
list(itertools.islice(fibonacci(0, 1), 100))
算术几何平均值
基于迭代计算算术和几何平均的迭代称为 AGM 迭代(更多信息请参考【1，第 598 页】):
它具有迷人的特性:
右边的积分叫做第一类完全椭圆积分。我们现在开始计算这个椭圆积分。我们使用生成器来描述迭代:
def arithmetic_geometric_mean(a, b):
"""
Generator for the arithmetic and geometric mean
a, b initial values
"""
while True: # infinite loop
a, b = (a+b)/2, sqrt(a*b)
yield a, b
当序列{ a i }收敛时，由{ c i }定义的序列{ c i }收敛到 0-这一事实将用于终止计算椭圆积分的程序中的迭代:
def elliptic_integral(k, tolerance=1e-5):
"""
Compute an elliptic integral of the first kind.
"""
a_0, b_0 = 1., sqrt(1-k**2)
for a, b in arithmetic_geometric_mean(a_0, b_0):
if abs(a-b) < tolerance:
return pi/(2*a)
我们必须确保算法停止。请注意，该代码完全依赖于算术几何平均迭代收敛(快速)的数学陈述。在实际计算中，我们在应用理论结果时必须小心，因为它们在有限精度的算术中可能不再有效。使前面的代码安全的正确方法是使用itertools.islice。安全代码如下(参见控制回路内流量一节下的示例，了解for / else语句的另一个典型用法):
from itertools import islice
def elliptic_integral(k, tolerance=1e-5, maxiter=100):
"""
Compute an elliptic integral of the first kind.
"""
a_0, b_0 = 1., sqrt(1-k**2)
for a, b in islice(arithmetic_geometric_mean(a_0, b_0),
maxiter):
if abs(a-b) < tolerance:
return pi/(2*a)
else:
raise Exception("Algorithm did not converge")
作为一种应用，椭圆积分可用于计算长度为 L 的摆的周期，该摆以角度θ开始(更多信息，请参考【18，第 114 页】),使用:
使用这个公式，钟摆的周期很容易得到:
def pendulum_period(L, theta, g=9.81):
return 4*sqrt(L/g)*elliptic_integral(sin(theta/2))
收敛加速度
我们给出了一个应用发电机加速收敛的例子。本演示紧跟 Python 生成器技巧中 Pramode C.E 给出的示例(更多信息请参考【9】)。
请注意，一个生成器可能会将另一个生成器作为输入参数。例如，假设我们已经定义了一个生成收敛序列元素的生成器。然后，由于欧拉和艾特肯，通常称为艾特肯的δ2-方法(参考【33】*)，有可能通过加速技术来改善收敛。*它通过定义将序列 s i 转换为另一个序列
两个序列具有相同的极限，但是序列收敛得明显更快。一种可能的实现如下:
def Euler_accelerate(sequence):
"""
Accelerate the iterator in the variable `sequence`.
"""
s0 = next(sequence) # Si
s1 = next(sequence) # Si+1
s2 = next(sequence) # Si+2
while True:
yield s0 - ((s1 - s0)**2)/(s2 - 2*s1 + s0)
s0, s1, s2 = s1, s2, next(sequence)
例如，我们使用经典系列:
向π/4 *收敛。*我们在下面的代码中将这个系列实现为一个生成器:
def pi_series():
sum = 0.
j = 1
for i in itertools.cycle([1, -1]):
yield sum
sum += i/j
j += 2
我们现在可以使用这个序列的加速版本:
Euler_accelerate(pi_series())
因此，该加速序列的第一 N 元素通过以下方式获得:
itertools.islice(Euler_accelerate(pi_series()), N)
例如，下图(图 9.1 )显示了由上述公式定义的序列的标准版本及其加速版本的误差对数的收敛速度:
图 9.1:定义的序列与其加速版本之间的比较
列出填充模式
在本节中，我们将比较不同的列表填充方式。它们在计算效率和代码可读性方面是不同的。
列表填充用追加法
无处不在的编程模式是计算元素并将它们存储在列表中:
L = []
for k in range(n):
# call various functions here
# that compute "result"
L.append(result)
这种方法有许多缺点:
- 迭代的次数是预先决定的。如果有
break指令，那么前面的代码负责生成值和决定何时停止。这是不可取的，也缺乏灵活性。 - 它假设用户想要所有迭代的整个计算历史。假设我们只对所有计算值的总和感兴趣。如果有许多计算值，存储它们是没有意义的，因为一次添加一个会更有效。
迭代器列表
迭代器为我们之前讨论的问题提供了一个优雅的解决方案:
def result_iterator():
for k in itertools.count(): # infinite iterator
# call various functions here
# that compute "result"
...
yield result
使用迭代器，我们将生成计算值的任务分开，而不用担心停止条件或存储。如果该代码的用户想要存储 n 个第一值，可以使用list构造函数轻松完成:
L = list(itertools.islice(result_iterator(), n)) # no append needed!
如果用户想要第一个 n 个生成值的总和，建议采用以下结构:
# make sure that you do not use scipy.sum here
s = sum(itertools.islice(result_iterator(), n))
我们在这里做的是一方面分离元素的生成，另一方面存储这些元素。
如果目的确实是建立一个列表，并且当每一步的结果不依赖于先前计算的元素时，可以使用列表理解语法(更多信息，请参考第 3 章、容器类型的列表部分):
L = [some_function(k) for k in range(n)]
当迭代计算依赖于先前计算值的值时，列表理解没有帮助。
存储生成的值
使用迭代器来填充列表在大多数情况下都会很好地工作，但是当计算新值的算法容易抛出异常时，这种模式就会变得复杂；如果迭代器在这个过程中引发了异常，那么列表将不可用！下面的例子说明了这个问题。
假设我们生成由递归定义的序列。如果初始数据 u 0 大于 1，则该序列迅速发散至无穷大。让我们用一个生成器来生成它:
import itertools
def power_sequence(u0):
u = u0
while True:
yield u
u = u**2
如果试图通过执行获得序列的第一个 20 元素(由 u 0 = 2 初始化)，
list(itertools.islice(power_sequence(2.), 20))
将引发异常，并且没有可用的列表，甚至没有引发异常之前的元素列表。目前还没有办法从可能有故障的发电机获得一个部分填充的列表。唯一的办法是使用封装在异常捕获块中的追加方法(更多详细信息，请参考第 10 章、错误处理中的异常一节):
generator = power_sequence(2.)
L = []
for iteration in range(20):
try:
L.append(next(generator))
except Exception:
...
当迭代器表现为列表时
一些列表操作也适用于迭代器。我们现在将检查列表理解和列表压缩的等效项(更多详细信息，请参考第 3 章、容器类型的列表部分)。
生成器表达式
这相当于对生成器的列表理解。这样的构造称为生成器表达式:
g = (n for n in range(1000) if not n % 100)
# generator for 100, 200, ... , 900
这对于计算和或积特别有用，因为那些运算是增量的；他们一次只需要一个元素:
sum(n for n in range(1000) if not n % 100) # returns 4500
在这段代码中，您注意到sum函数有一个参数，这是一个生成器表达式。请注意，当生成器仅用作函数的参数时，Python 语法允许我们省略生成器的括号。
*让我们计算黎曼泽塔函数 ζ ，其表达式为
使用生成器表达式，我们可以在一行中计算这个系列的部分和:
sum(1/n**s for n in itertools.islice(itertools.count(1), N))
请注意，我们还可以定义序列 1nnsT3 的生成器，如下所示:
def generate_zeta(s):
for n in itertools.count(1):
yield 1/n**s
然后我们简单地获得第一个 N 个项的和，使用:
def zeta(N, s):
# make sure that you do not use the scipy.sum here
return sum(itertools.islice(generate_zeta(s), N))
我们指出，我们使用这种计算ζ(ζ)函数的方法，以优雅的方式演示了发电机的使用。这当然不是评估该函数的最准确和计算效率最高的方法。
压缩迭代器
我们在部分列表、第 3 章、容器类型中看到，可以通过将两个容器拉在一起来创建一个列表。迭代器也存在相同的操作:
xg = x_iterator() # some iterator
yg = y_iterator() # another iterator
for x, y in zip(xg, yg):
print(x, y)
一旦其中一个迭代器用完，压缩迭代器就会停止。这与列表上的压缩操作的行为相同。
迭代器对象
正如我们前面提到的，一个for循环只需要一个可迭代的对象。尤其是列表是可重复的。这意味着列表能够从其内容创建迭代器。事实上，这对于任何对象(不仅仅是列表)都是正确的:任何对象都可能是可重复的。
这是通过__iter__方法实现的，该方法应该返回一个迭代器。这里我们举一个例子，其中__iter__方法是一个生成器:
class OdeStore:
"""
Class to store results of ode computations
"""
def __init__(self, data):
"data is a list of the form [[t0, u0], [t1, u1],...]"
self.data = data
def __iter__(self):
"By default, we iterate on the values u0, u1,..."
for t, u in self.data:
yield u
store = OdeStore([[0, 1], [0.1, 1.1], [0.2, 1.3]])
for u in store:
print(u)
# result: 1, 1.1, 1.3
list(store) # [1, 1.1, 1.3]
如果您试图对不可迭代的对象使用迭代器的功能，将会引发异常:
>>> list(3)
TypeError: 'int' object is not iterable
在这个例子中，列表函数试图通过调用__iter__方法迭代对象 3 。但是这个方法不是为整数实现的，因此引发了异常。如果我们试图循环通过一个不可迭代的对象，也会发生同样的情况:
>>> for iteration in 3: pass
TypeError: 'int' object is not iterable
无限迭代
无限迭代或者通过无限迭代器，通过while循环，或者通过递归获得。显然，在实际情况下，某些条件会停止迭代。有限迭代的区别在于，通过粗略地检查代码，不可能说迭代是否会停止。
while 循环
while循环可用于重复一个代码块，直到满足一个条件:
while condition:
<code>
一个while循环相当于以下代码:
for iteration in itertools.count():
if not condition:
break
<code>
因此while循环相当于无限迭代器，如果满足某个条件，可能会停止。这种构造的危险是显而易见的:如果条件永远不满足，代码可能会陷入无限循环。
科学计算的问题在于，人们并不总是确定算法会收敛。例如，牛顿迭代可能根本不会收敛。如果该算法在一个while循环中实现，那么对于一些初始条件的选择，相应的代码将陷入无限循环。
因此，我们建议有限迭代器通常更适合这样的任务。以下结构通常有利地取代了while环的使用:
maxit = 100
for nb_iterations in range(max_it):
...
else:
raise Exception("No convergence in {} iterations".format(maxit))
第一个优点是，无论发生什么，代码都保证在有限的时间内执行。第二个优点是变量nb_iterations包含算法收敛所需的迭代次数。
递归
当一个函数调用自己时，就会发生递归(参见第 7 章、函数中递归函数一节)。
在进行递归时，递归深度，也就是迭代的次数，将你的计算机带到了极限。我们在这里通过考虑一个简单的递归来演示这一点，它实际上根本不包含任何计算。它只给迭代赋值零:
def f(N):
if N == 0:
return 0
return f(N-1)
根据您的系统，该程序可能会阻塞 N ≥ 10000 (使用了太多内存)。结果是 Python 解释器崩溃了，没有进一步的异常。Python 提供了一种机制，可以在检测到过高的递归深度时引发异常。该最大递归深度可以通过执行以下操作来更改:
import sys
sys.setrecursionlimit(1000)
递归极限的实际值可以通过sys.getrecursionlimit()得到。
但是请注意，选择过高的数字可能会危及代码的稳定性，因为 Python 可能会在达到最大深度之前崩溃。因此，保持递归极限不变通常是明智的。
相比之下，以下非递归程序运行千万次迭代没有任何问题:
for iteration in range(10000000):
pass
我们主张，如果可能的话，应该在 Python 中避免递归。这显然只适用于有合适的替代迭代算法可用的情况。第一个原因是深度 N 的递归同时涉及 N 函数调用，这可能会导致很大的开销。第二个原因是它是一个无限迭代，也就是说，在递归结束之前，很难给出必要步骤数的上限。
请注意，在一些非常特殊的情况下(树遍历)，递归是不可避免的。此外，在某些情况下(递归深度很小)，递归程序由于可读性可能更好。
总结
在这一章中，我们研究了迭代器，这是一种非常接近迭代方法的数学描述的编程构造。你看到了yield关键字，遇到了有限和无限迭代器。
我们展示了迭代器可以被耗尽。更多的特殊方面，如迭代器理解和递归迭代器被引入，并通过例子进行了演示。
练习
Ex。1 →计算总和的值:
Ex。2 →创建一个生成器，用于计算由以下关系定义的序列:
Ex。3 →生成所有偶数。
Ex。4 →让。微积分中显示
。通过实验确定最小数量 n ，使得
。为此任务使用生成器。
Ex。5 →生成所有小于给定整数的素数。使用名为厄拉多塞筛的算法。
Ex。6 →应用显式欧拉方法求解微分方程得到递归:
编写一个生成器，计算给定初始值 u 0 和给定时间步长值 h 的解数值 u n 。
Ex。7 →使用公式计算π:
积分可以用复合梯形法则来近似，即通过以下公式:
其中。
为值 y i = f(x i ) 编程一个生成器，并通过对一个又一个项求和来评估公式。将您的结果与 SciPy 的quad功能进行比较。
Ex。8 →让 x = [1，2，3]和 y = [-1，-2，-3]。代码zip(*zip(x, y))有什么作用？解释它是如何工作的。
Ex。9 →完全椭圆积分可以通过函数scipy.special.ellipk计算。编写一个函数，计算 AGM 迭代所需的迭代次数，直到结果符合给定的公差(注意ellipk中的输入参数 m 对应于*算术几何平均)*一节定义中的 k 2 。
Ex。10 →考虑以下定义的顺序:
它单调收敛到零:E1T10】E2T11。。。> 0 。通过分部积分，可以看出序列 E n 实现了以下递归:
使用适当的生成器计算递归的前 20 项，并将结果与通过scipy.integrate.quad数值积分获得的结果进行比较。通过反转递归来执行相同的操作:
使用exp函数评估指数函数。你观察到了什么？你有解释吗？(参考【29】)
图 9.2:逼近 sin(x)的函数的收敛性研究
Ex。11 →正弦函数可由欧拉表示为
编写一个生成函数值 P k (x) 的生成器。设置x=linspace(-1,3.5*pi,200) 并通过图形演示 P k (x) 对增加 k 来说有多好。上图(图 9.2 )显示了可能的结果(参考[【11，Th。5.2，第 65 页]](16.html "Appendix . References") )。*
十、错误处理
在本章中，我们将介绍错误、异常以及如何查找和修复它们。处理异常是编写可靠和可用代码的重要部分。我们将介绍基本的内置异常，并展示如何使用和处理异常。我们将介绍调试，并向您展示如何使用内置的 Python 调试器。
有哪些例外？
程序员(即使是有经验的程序员)发现的一个错误是代码有不正确的语法，这意味着代码指令的格式不正确。
考虑一个语法错误的例子:
>>> for i in range(10)
File “<stdin>”, line 1
for i in range(10)
^
SyntaxError: invalid syntax
出现错误是因为for声明末尾缺少冒号。这是引发异常的一个示例。在SyntaxError的情况下，它告诉程序员代码有不正确的语法，并打印错误发生的行，箭头指向该行中的问题所在。
Python 中的异常是从名为Exception的基类派生(继承)而来的。Python 有许多内置的例外。一些常见的异常类型列于表 10.1、(内置异常的完整列表参见*【38】*)。
以下是两个常见的例外示例。如你所料，当你试图除以 0 时ZeroDivisionError被提升。
def f(x):
return 1/x
>>> f(2.5)
0.4
>>> f(0)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "exception_tests.py", line 3, in f
return 1/x
ZeroDivisionError: integer division or modulo by zero
表 10.1:一些常用的内置异常及其含义
零除会引发ZeroDivisionError并打印出发生错误的文件、行和函数名。
正如我们之前看到的，数组只能包含相同数据类型的元素。如果您试图分配一个不兼容类型的值，则会引发ValueError。值错误的示例:
>>> a = arange(8.0)
>>> a
array([ 0., 1., 2., 3., 4., 5., 6., 7.])
>>> a[3] = 'string'
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ValueError: could not convert string to float: string
这里，引发ValueError是因为数组包含浮点，并且不能为元素分配字符串值。
基本原则
让我们看看如何使用异常的基本原则，用raise引发异常，用try语句捕捉异常。
引发异常
创建错误被称为引发异常。在前一节中，您看到了一些异常示例。您也可以定义自己的预定义类型或无类型的例外。使用如下命令可以引发异常:
raise Exception("Something went wrong")
当出现问题时，很容易打印出错误消息，例如:
print("The algorithm did not converge.")
出于多种原因，不建议这样做。首先，打印输出很容易被遗漏，尤其是当信息被隐藏在打印到控制台的许多其他信息中时。其次，更重要的是，它会使您的代码无法被其他代码使用。调用代码将无法知道发生了错误，因此无法处理它。
由于这些原因，最好还是引发异常。异常应始终包含描述性消息，例如:
raise Exception("The algorithm did not converge.")
这条信息对用户来说会很明显。它还为调用代码提供了知道发生了错误并可能找到补救方法的机会。
下面是一个典型的例子，检查函数内部的输入，以确保它在继续之前是可用的。例如，对负值和正确数据类型的简单检查确保了计算阶乘的函数的预期输入:
def factorial(n):
if not (n >=0 and isinstance(n,(int,int32,int64))):
raise ValueError("A positive integer is expected")
...
如果给出了不正确的输入，函数的用户将立即知道错误是什么，处理异常是用户的责任。请注意在引发预定义异常类型时异常名称的使用，在这种情况下ValueError后跟消息。通过指定异常的类型，调用代码可以根据引发的错误类型决定以不同的方式处理错误。
总而言之，引发异常总是比打印错误消息更好。
捕捉异常
处理异常被称为捕获异常。使用try和except命令检查异常。
异常会停止程序执行流程，并寻找最近的try封闭块。如果没有捕获到异常，程序单元被留下，它继续在调用栈中较高的程序单元中搜索下一个封闭的try块。如果没有找到任何块，并且没有处理异常，则执行完全停止；将显示标准追溯信息。
让我们看一个try语句的例子:
try:
<some code that might raise an exception>
except ValueError:
print("Oops, a ValueError occurred")
在这种情况下，如果try块中的代码引发了类型为ValueError的错误，将会捕获异常并打印except块中的消息。如果在try块内没有发生异常，except块被完全跳过，继续执行。
except语句可以捕捉多个异常。这是通过简单地将它们分组到一个元组中来完成的，如下所示:
except (RuntimeError, ValueError, IOError):
try块也可以有多个except语句。这使得根据类型以不同的方式处理异常成为可能。让我们看一个多种异常类型的例子:
try:
f = open('data.txt', 'r')
data = f.readline()
value = float(data)
except OSError as oe:
print("{}: {}".format(oe.strerror, oe.filename))
except ValueError:
print("Could not convert data to float.")
例如，如果文件不存在，这里将捕捉到一个OSError；例如，如果文件第一行中的数据与浮点数据类型不兼容，则会捕获一个ValueError。
在本例中，我们通过关键字as将OSError分配给变量oe。这允许在处理该异常时访问更多细节。这里我们打印了错误字符串oe.strerror和相关文件的名称oe.filename。根据类型的不同，每种错误类型都可以有自己的一组变量。如果文件不存在，在前面的示例中，消息将是:
I/O error(2): No such file or directory
另一方面，如果文件存在，但您没有打开它的权限，消息将是:
I/O error(13): Permission denied
这是捕获异常时格式化输出的有用方法。
try - except组合可以通过可选的else和finally模块进行扩展。使用else的例子可以在第 13 章测试的测试等分算法部分看到。当最终需要进行清理工作时，将try和finally结合起来会给出一个有用的结构:
确保文件正确关闭的示例:
try:
f = open('data.txt', 'r')
# some function that does something with the file
process_file_data(f)
except:
...
finally:
f.close()
这将确保无论在处理文件数据时引发什么异常，文件最终都是关闭的。在try语句中没有处理的异常在finally块之后被保存和引发。这个组合用在with语句中；请参见上下文管理器-附带声明一节。
用户定义的异常
除了内置的 Python 异常，还可以定义自己的异常。这种用户定义的异常应该继承自Exception基类。当您定义自己的类时，如第 14 章综合示例的多项式一节中的多项式类，这可能会很有用。
看看这个简单的用户定义异常的小例子:
class MyError(Exception):
def __init__(self, expr):
self.expr = expr
def __str__(self):
return str(self.expr)
try:
x = random.rand()
if x < 0.5:
raise MyError(x)
except MyError as e:
print("Random number too small", e.expr)
else:
print(x)
生成一个随机数。如果该数字低于 0.5，将引发异常，并打印一条消息，指出该值太小。如果没有引发异常，则打印该数字。
在本例中，您还看到了在try语句中使用else的情况。如果没有异常发生，else下的区块将被执行。
建议您使用以Error结尾的名称定义异常，如标准内置异常的命名。
上下文管理器 with 语句
Python 中有一个非常有用的构造，用于在处理上下文(如文件或数据库)时简化异常处理。该语句将try ... finally结构封装在一个简单的命令中。下面是一个使用with读取文件的例子:
with open('data.txt', 'r') as f:
process_file_data(f)
这将尝试打开文件，对文件运行指定的操作(例如，读取)，并关闭文件。如果在执行process_file_data的过程中出现任何问题，文件会被正确关闭，然后引发异常。这相当于:
f = open('data.txt', 'r')
try:
# some function that does something with the file
process_file_data(f)
except:
...
finally:
f.close()
我们将在第 12 章、输入输出的文件处理一节中，在读写文件时使用这个选项。
前面的文件读取示例是使用上下文管理器的示例。上下文管理器是 Python 对象，有两种特殊的方法，_ _enter_ _和_ _exit_ _。实现这两种方法的类的任何对象都可以用作上下文管理器。在这个例子中，文件对象f是一个上下文管理器，因为有f._ _enter_ _和f._ _exit_ _方法。
_ _enter_ _方法应该实现初始化指令，比如打开文件或者数据库连接。如果此方法有一个返回语句，则使用as构造访问返回的对象。否则省略as关键字。_ _exit_ _方法包含清理指令，例如，关闭文件或提交事务并关闭数据库连接。如需了解更多解释和自写上下文管理器示例，请参阅第 13 章、测试的使用上下文管理器计时一节。
有一些 NumPy 函数可以用作上下文管理器。例如load功能支持某些文件格式的上下文管理器。NumPy 的函数errstate可以用作上下文管理器，指定代码块内的浮点错误处理行为。
下面是使用 errstate 和上下文管理器的一个示例:
import numpy as np # note, sqrt in NumPy and SciPy
# behave differently in that example
with errstate(invalid='ignore'):
print(np.sqrt(-1)) # prints 'nan'
with errstate(invalid='warn'):
print(np.sqrt(-1)) # prints 'nan' and
# 'RuntimeWarning: invalid value encountered in sqrt'
with errstate(invalid='raise'):
print(np.sqrt(-1)) # prints nothing and raises FloatingPointError
参考第 2 章、变量和基本类型的无限和非数字部分，了解本例的更多细节，参考第 13 章、测试的部分，了解另一个示例。
查找错误:调试
软件代码中的错误有时被称为bug。调试是发现和修复代码中的错误的过程。这个过程可以在不同的复杂程度下进行。最有效的方法是使用一种叫做调试器的工具。进行单元测试是早期识别错误的好方法，请参考第 13 章、测试的使用单元测试一节。当不清楚问题在哪里或者是什么的时候，调试器是非常有用的。
虫子
通常有两种错误:
- 会引发异常，但不会被捕获。
- 代码运行不正常。
第一种情况通常更容易解决。第二个可能更难，因为问题可能是一个错误的想法或解决方案，一个错误的实现，或者两者的结合。
我们只关心接下来的第一种情况，但是同样的工具可以用来帮助找到为什么代码没有做它应该做的事情。
堆栈
当引发异常时，您会看到调用堆栈。调用堆栈包含调用引发异常的代码的所有函数的跟踪。
一个简单的堆栈示例:
def f():
g()
def g():
h()
def h():
1//0
f()
这种情况下的堆栈是f、 g和h。运行这段代码生成的输出如下所示:
Traceback (most recent call last):
File "stack_example.py", line 11, in <module>
f()
File "stack_example.py", line 3, in f
g()
File "stack_example.py", line 6, in g
h() File "stack_example.py", line 9, in h
1//0
ZeroDivisionError: integer division or modulo by zero
打印错误。显示了导致错误的函数序列。第 11 行的函数f被调用，依次调用g和h。这就造成了ZeroDivisionError。
堆栈跟踪报告程序执行中某一点的活动堆栈。堆栈跟踪允许您跟踪调用到给定点的函数序列。通常这是在引发了未捕获的异常之后。这有时被称为事后分析，堆栈跟踪点就是异常发生的地方。另一种选择是手动调用堆栈跟踪来分析您怀疑有错误的代码段，可能是在异常发生之前。
Python 调试器
Python 自带一个名为 pdb 的内置调试器。一些开发环境集成了调试器。以下过程在大多数情况下仍然适用。
使用调试器最简单的方法是在您想要调查的代码点启用堆栈跟踪。下面是基于第七章、函数的返回值一节中提到的例子触发调试器的一个简单例子:
import pdb
def complex_to_polar(z):
pdb.set_trace()
r = sqrt(z.real ** 2 + z.imag ** 2)
phi = arctan2(z.imag, z.real)
return (r,phi)
z = 3 + 5j
r,phi = complex_to_polar(z)
print(r,phi)
pdb.set_trace()命令启动调试器并启用后续命令的跟踪。前面的代码将显示这一点:
> debugging_example.py(7)complex_to_polar()
-> r = sqrt(z.real ** 2 + z.imag ** 2)
(Pdb)
调试器提示用(Pdb)表示。调试器停止程序执行，并给你一个提示，让你检查变量，修改变量，单步执行命令，等等。
每一步都会打印当前行，因此您可以跟踪自己的位置以及接下来会发生什么。通过命令n(下一步)来执行命令，如下所示:
> debugging_example.py(7)complex_to_polar()
-> r = sqrt(z.real ** 2 + z.imag ** 2)
(Pdb) n
> debugging_example.py(8)complex_to_polar()
-> phi = arctan2(z.imag, z.real)
(Pdb) n
> debugging_example.py(9)complex_to_polar()
-> return (r,phi)
(Pdb)
...
命令n(下一步)将继续到下一行并打印该行。如果需要同时查看多行，列表命令l(列表)显示当前行及周围代码:
在调试器中列出周围的代码:
> debugging_example.py(7)complex_to_polar()
-> r = sqrt(z.real ** 2 + z.imag ** 2)
(Pdb) l
2
3 import pdb
4
5 def complex_to_polar(z):
6 pdb.set_trace()
7 -> r = sqrt(z.real ** 2 + z.imag ** 2)
8 phi = arctan2(z.imag, z.real)
9 return (r,phi)
10
11 z = 3 + 5j
12 r,phi = complex_to_polar(z)
(Pdb)
变量的检查可以通过使用命令p(打印)后跟变量名将其值打印到控制台来完成。打印变量的示例:
> debugging_example.py(7)complex_to_polar()
-> r = sqrt(z.real ** 2 + z.imag ** 2)
(Pdb) p z
(3+5j) (Pdb) n
> debugging_example.py(8)complex_to_polar()
-> phi = arctan2(z.imag, z.real)
(Pdb) p r
5.8309518948453007
(Pdb) c
(5.8309518948453007, 1.0303768265243125)
p(打印)命令将打印变量；命令c(继续)继续执行。
在执行过程中更改变量是有用的。只需在调试器提示符下分配新值，然后单步执行或继续执行。
一个改变变量的例子:
> debugging_example.py(7)complex_to_polar()
-> r = sqrt(z.real ** 2 + z.imag ** 2)
(Pdb) z = 2j
(Pdb) z
2j
(Pdb) c
(2.0, 1.5707963267948966)
这里，变量z被赋予一个新值，在剩余代码中使用。请注意，最终打印结果已经更改。
概述-调试命令
在表 10.2 中，显示了最常见的调试命令。有关命令的完整列表和描述，请参见文档【25】了解更多信息。请注意，任何 Python 命令也可以工作，例如，为变量赋值。
型式
短变量名
如果要检查名称与调试器的任何短命令一致的变量，例如h，必须使用!h显示该变量。
表 10.2:调试器最常见的调试命令。
IPython 中的调试
IPython 附带了名为ipdb的调试器版本。在撰写本书时，差异非常小，但这可能会改变。
IPython 中有一个命令，在出现异常时自动打开调试器。这在尝试新想法或代码时非常有用。如何在 IPython 中自动打开调试器的示例:
In [1]: %pdb # this is a so - called IPython magic command
Automatic pdb calling has been turned ON
In [2]: a = 10
In [3]: b = 0
In [4]: c = a/b
___________________________________________________________________
ZeroDivisionError Traceback (most recent call last)
<ipython-input-4-72278c42f391> in <module>()
—-> 1 c = a/b
ZeroDivisionError: integer division or modulo by zero
> <ipython-input-4-72278c42f391>(1)<module>()
-1 c = a/b
ipdb>
IPython 提示符下的 IPython 魔法命令%pdb在引发异常时自动启用调试器。这里调试器提示显示ipdb来指示调试器正在运行。
总结
本章中的关键概念是异常和错误。我们展示了一个异常是如何在另一个程序单元中被捕获的。您可以定义自己的异常，并为它们配备消息和给定变量的当前值。
代码可能会返回意外结果，但不会引发异常。定位错误结果来源的技术称为调试。我们介绍了调试方法，并希望鼓励您对它们进行培训，以便在需要时随时可用。认真调试的需求来得比你预期的要早。
十一、名称空间、范围和模块
在本章中，我们将介绍 Python 模块。模块是包含函数和类定义的文件。本章还解释了命名空间的概念以及跨函数和模块的变量范围。
命名空间
Python 对象的名称，如变量、类、函数和模块的名称，都收集在名称空间中。模块和类有自己的命名空间，与这些对象同名。这些命名空间是在导入模块或实例化类时创建的。模块命名空间的生存期与当前 Python 会话一样长。类实例的命名空间的生存期到该实例被删除为止。
函数在执行(调用)时会创建一个本地命名空间。当函数通过常规返回或异常停止执行时，它被删除。本地命名空间未命名。
名称空间的概念将变量名放在其上下文中。例如，有几个名为sin的函数，它们通过它们所属的命名空间来区分，如下面的代码所示:
import math
import scipy
math.sin
scipy.sin
它们确实不同，因为scipy.sin是一个接受列表或数组作为输入的通用函数，其中math.sin只接受浮点数。命令dir(<name of the namespace>)可以获得一个包含特定名称空间中所有名称的列表。它包含两个特殊的名字__name__和__doc__。前者指的是模块的名称，后者指的是它的文档字符串:
math.__name__ # returns math
math.__doc__ # returns 'This module is always ...'
有一个特殊的命名空间__builtin__，它包含 Python 中可用的名称，没有任何import。它是一个命名空间，但是在引用内置对象时不需要给出它的名称:
'float' in dir(__builtin__) # returns True
float is __builtin__.float # returns True
变量的范围
在程序的一个部分定义的变量不需要在其他部分知道。已知某个变量的所有程序单元称为该变量的作用域。我们先举个例子；让我们考虑两个嵌套函数:
e = 3
def my_function(in1):
a = 2 * e
b = 3
in1 = 5
def other_function():
c = a
d = e
return dir()
print("""
my_function's namespace: {}
other_function's namespace: {}
""".format(dir(),other_function()))
return a
执行my_function(3)会导致:
my_function's namespace: ['a', 'b', 'in1', 'other_function']
other_function's namespace: ['a', 'c', 'd']
变量e位于包含函数my_function 的程序单元的名称空间中。变量a在这个函数的名称空间中，它本身包含了最里面的函数other_function。对于这两个函数，e是一个全局变量。
一个好的做法是只通过参数列表向函数传递信息，而不使用前面例子中的构造。在第七章、函数的匿名函数一节中可以找到一个例外，其中全局变量用于闭包。通过给变量赋值，变量会自动变成局部变量:
e = 3
def my_function():
e = 4
a = 2
print("my_function's namespace: {}".format(dir()))
执行
e = 3
my_function()
e # has the value 3
给出:
my_function's namespace: ['a', 'e']
其中e变成了局部变量。事实上，这段代码现在有两个属于不同名称空间的变量e。
通过使用global声明语句，函数中定义的变量可以成为全局变量，也就是说，它的值甚至可以在该函数之外访问。global声明的使用演示如下:
def fun():
def fun1():
global a
a = 3
def fun2():
global b
b = 2
print(a)
fun1()
fun2() # prints a
print(b)
型式
避免使用全局
建议避免使用这种构造和global的使用。这类代码很难调试和维护。类的使用(详见第八章、类，使得global主要过时。
模块
在 Python 中，模块只是一个包含类和函数的文件。通过在会话或脚本中导入文件，函数和类变得可用。
简介
默认情况下，Python 附带许多不同的库。您可能还想为特定目的安装更多这样的软件，例如优化、绘图、读/写文件格式、图像处理等。NumPy 和 SciPy 是这种库的两个重要例子，matplotlib 用于绘图是另一个例子。在本章的最后，我们将列出一些有用的库。
要使用库，您可以:
-
仅从库中加载某些对象，例如从 NumPy:
from numpy import array, vander -
或者加载整个库:
from numpy import * -
Or give access to an entire library by creating a namespace with the library name:
import numpy ... numpy.array(...)在库中的函数前面加上名称空间，就可以访问该函数，并将该函数与其他同名对象区分开来。
此外，命名空间的名称可以与import命令一起指定:
import numpy as np
...
np.array(...)
您使用的选项会影响代码的可读性以及出错的可能性。一个常见的错误是阴影:
from scipy.linalg import eig
A = array([[1,2],[3,4]])
(eig, eigvec) = eig(A)
...
(c, d) = eig(B) # raises an error
避免这种意外影响的方法是使用import:
import scipy.linalg as sl
A = array([[1,2],[3,4]])
(eig, eigvec) = sl.eig(A) # eig and sl.eig are different objects
...
(c, d) = sl.eig(B)
在本书中，我们使用了许多命令、对象和函数。这些是通过以下语句导入本地命名空间的:
from scipy import *
以这种方式导入对象不会使从中导入对象的模块变得明显。下表给出了一些例子(表 11.1 ):
| **图书馆** | **方法** | | `numpy` | `array`、`arange`、`linspace`、`vstack`、`hstack`、`dot`、`eye`、`identity`和`zeros`。 | | `numpy.linalg` | `solve`、`lstsq`、`eig`和`det`。 | | `matplotlib.pyplot` | `plot`、`legend`和`cla`。 | | `scipy.integrate` | `quad`。 | | `copy` | `copy`和`deepcopy`。 |表 11.1:导入对象的示例
IPython 中的模块
IPython 在代码开发下使用。一个典型的场景是，您处理一个文件，其中包含一些您在开发周期中更改的函数或类定义。要将此类文件的内容加载到 shell 中，可以使用import，但文件只加载一次。更改文件对以后的导入没有影响。这就是 IPyhthon 的魔法指令run登场的地方。
IPython 魔法命令
IPython 有一个名为run的特殊魔法命令，它会执行一个文件，就像你在 Python 中直接运行它一样。这意味着文件的执行独立于 IPython 中已经定义的内容。当您想要测试作为独立程序的脚本时，这是从 IPython 中执行文件的推荐方法。您必须以与从命令行执行相同的方式在执行的文件中导入所有需要的内容。在myfile.py中运行代码的一个典型例子是:
from numpy import array
...
a = array(...)
该脚本文件由exec(open('myfile.py').read())在 Python 中执行。或者，在 IPython 中，如果您想确保脚本独立于之前的导入运行，可以使用魔法命令run myfile。文件中定义的所有内容都被导入到 IPython 工作区。
变量 name
在任何模块中，特殊变量__name__被定义为当前模块的名称。在命令行中(在 IPython 中)，该变量被设置为__main__，这允许以下技巧:
# module
import ...
class ...
if __name__ == "__main__":
# perform some tests here
测试将仅在文件直接运行时运行，而不是在文件导入时运行。
一些有用的模块
有用的 Python 模块有很多。在下表中，我们给出了一个非常简短的列表，重点是与数学和工程应用相关的模块(表 11.2) :
| **模块** | **描述** | | `scipy` | 科学计算中使用的函数 | | `numpy` | 支持数组和相关方法 | | `matplotlib` | 使用导入子模块 pyplot 进行绘图和可视化 | | `functools` | 函数的部分应用 | | `itertools` | 迭代器工具提供特殊的功能，比如对生成器进行切片 | | `re` | 高级字符串处理的正则表达式 | | `sys` | 系统特定功能 | | `os` | 操作系统界面，如目录列表和文件处理 | | `datetime` | 表示日期和日期增量 | | `time` | 返回挂钟时间 | | `timeit` | 测量执行时间 | | `sympy` | 计算机算术包(符号计算) | | `pickle` | 酸洗，特殊文件输入和输出格式 | | `shelves` | 架子，特殊文件输入和输出格式 | | `contextlib` | 上下文管理器工具 |表 11.2:工程应用中有用的 Python 包的非穷举列表
总结
我们在书的开头告诉你，你必须导入 SciPy 和其他有用的模块。现在你完全明白进口意味着什么。我们介绍了名称空间，并讨论了import和from ... import *之间的区别。变量的范围已经在早期的第 7 章、函数中介绍过了，但是现在你对这个概念的重要性有了更完整的了解。
十二、输入和输出
在本章中，我们将介绍一些处理数据文件的选项。根据数据和所需的格式，有几种读写选项。我们将展示一些最有用的替代方案。
文件处理
在许多情况下，文件输入/输出(输入和输出)是必不可少的。例如:
- 使用测量或扫描的数据。测量值存储在需要读取以进行分析的文件中。
- 与其他程序交互。将结果保存到文件中，以便在其他应用中导入，反之亦然。
- 存储信息供将来参考或比较。
- 与他人共享数据和结果，可能在其他平台上使用其他软件。
在本节中，我们将介绍如何用 Python 处理文件输入/输出。
与文件交互
在 Python 中，file类型的对象表示存储在磁盘上的物理文件的内容。可以使用以下语法创建新的file对象:
myfile = open('measurement.dat','r') # creating a new file object from an existing file
例如，可以通过以下方式访问文件的内容:
print(myfile.read())
文件对象的使用需要小心。问题是文件必须先关闭，然后才能被其他应用重新读取或使用，这是使用以下语法完成的:
myfile.close() # closes the file object
然而，事情并没有那么简单，因为在执行对close的调用之前可能会触发异常，这将跳过结束代码(考虑以下示例)。确保文件被正确关闭的一个简单方法是使用上下文管理器。在第 10 章、错误处理的例外一节中，使用with关键字对该结构进行了更详细的解释。以下是它如何用于文件:
with open('measurement.dat','r') as myfile:
... # use myfile here
这确保了当一个人退出with块时文件是关闭的，即使在块内部引发了异常。该命令适用于上下文管理器对象。我们建议您在第 10 章、错误处理的例外部分阅读更多关于上下文管理器的内容。这里有一个例子说明了为什么with结构是可取的:
myfile = open(name,'w')
myfile.write('some data')
a = 1/0
myfile.write('other data')
myfile.close()
文件关闭前会引发异常。文件保持打开状态，并且不能保证文件中写入了什么数据或何时写入。因此，实现相同结果的正确方法是:
with open(name,'w') as myfile:
myfile.write('some data')
a = 1/0
myfile.write('other data')
在这种情况下，文件会在引发异常(这里是ZeroDivisionError)后干净地关闭。还要注意，没有必要显式关闭文件。
文件是可重复的
特别是，文件是可迭代的(参见第 9 章、迭代的迭代器部分)。文件重复它们的行:
with open(name,'r') as myfile:
for line in myfile:
data = line.split(';')
print('time {} sec temperature {} C'.format(data[0],data[1]))
文件的行作为字符串返回。字符串方法split是将字符串转换为字符串列表的可能工具。例如:
data = 'aa;bb;cc;dd;ee;ff;gg'
data.split(';') # ['aa', 'bb', 'cc', 'dd', 'ee', 'ff', 'gg']
data = 'aa bb cc dd ee ff gg'
data.split(' ') # ['aa', 'bb', 'cc', 'dd', 'ee', 'ff', 'gg']
由于myfile对象是可迭代的，我们也可以直接提取到列表中，如下所示:
data = list(myfile)
文件模式
从这些文件处理的例子中可以看出，open函数至少有两个参数。第一个显然是文件名，第二个是描述文件使用方式的字符串。打开文件有几种这样的模式；基本的有:
with open('file1.dat','r') as ... # read only
with open('file2.dat','r+') as ... # read/write
with open('file3.dat','rb') as ... # read in byte mode
with open('file4.dat','a') as ... # append (write to the end of the file)
with open('file5.dat','w') as ... # (over-)write the file
with open('file6.dat','wb') as ... # (over-)write the file in byte mode
'r'、'r+'和'a'模式要求文件存在，而'w'如果不存在同名文件，将创建一个新文件。用'r'和'w'阅读和写作是最常见的，正如你在前面的例子中看到的。
考虑一个例子，使用追加'a'模式打开一个文件并在文件末尾添加数据，而不修改已经存在的内容。注意断线，\n:
with open('file3.dat','a') as myfile:
myfile.write('something new\n')
NumPy 方法
NumPy 内置了将 NumPy 数组数据读写到文本文件的方法。这些是numpy.loadtxt和numpy.savetxt。
savetxt
将数组写入文本文件很简单:
savetxt(filename,data)
有两个有用的参数作为字符串给出，fmt和delimiter，它们控制格式和列之间的分隔符。默认值是分隔符为空格，格式为%.18e，对应于所有数字的指数格式。格式参数的使用如下:
x = range(100) # 100 integers
savetxt('test.txt',x,delimiter=',') # use comma instead of space
savetxt('test.txt',x,fmt='%d') # integer format instead of float with e
loadtxt
从文本文件读取数组是在以下语法的帮助下完成的:
filename = 'test.txt'
data = loadtxt(filename)
由于数组中的每一行都必须具有相同的长度，因此文本文件中的每一行都必须具有相同数量的元素。与savetxt类似，默认值为float，分隔符为space。可以使用dtype和delimiter参数进行设置。另一个有用的参数是comments，可以用来标记数据文件中的注释使用了什么符号。使用格式化参数的示例如下:
data = loadtxt('test.txt',delimiter=';') # data separated by semicolons
data = loadtxt('test.txt',dtype=int,comments='#') # read to integer type,
#comments in file begin with a hash character
腌制
您刚才看到的读写方法在写入之前将数据转换为字符串。复杂类型(如对象和类)不能用这种方式编写。使用 Python 的 pickle 模块，您可以将任何对象以及多个对象保存到文件中。
数据可以以明文(ASCII)格式保存，也可以使用稍微高效的二进制格式保存。主要有两种方法:dump，将 Python 对象的腌制表示保存到文件中；以及load，从文件中检索腌制对象。基本用法是这样的:
import pickle
with open('file.dat','wb') as myfile:
a = random.rand(20,20)
b = 'hello world'
pickle.dump(a,myfile) # first call: first object
pickle.dump(b,myfile) # second call: second object
import pickle
with open('file.dat','rb') as myfile:
numbers = pickle.load(myfile) # restores the array
text = pickle.load(myfile) # restores the string
请注意这两个对象的返回顺序。除了两种主要方法之外，将 Python 对象序列化为字符串而不是文件有时也很有用。这是通过dumps和load完成的。考虑一个序列化数组和字典的示例:
a = [1,2,3,4]
pickle.dumps(a) # returns a bytes object
b = {'a':1,'b':2}
pickle.dumps(b) # returns a bytes object
使用dumps的一个很好的例子是当你需要将 Python 对象或 NumPy 数组写入数据库时。这些通常支持存储字符串，这使得无需任何特殊模块就可以轻松地写入和读取复杂的数据和对象。除了 pickle 模块，还有一个优化版叫做cPickle。它是用 C 写的，如果你需要快速阅读和写作，它是一个选项。pickle 和 cPickle 产生的数据相同，可以互换。
货架
字典中的对象可以通过键来访问。访问文件中的特定数据也有类似的方法，首先为其分配一个密钥。这可以通过使用模块架来实现:
from contextlib import closing
import shelve as sv
# opens a data file (creates it before if necessary)
with closing(sv.open('datafile')) as data:
A = array([[1,2,3],[4,5,6]])
data['my_matrix'] = A # here we created a key
在文件处理一节中，我们看到内置的open命令生成了一个上下文管理器，我们看到了为什么这对处理外部资源(如文件)很重要。与此命令相反，sv.open本身并不创建上下文管理器。需要来自contextlib模块的closing命令将其转换为合适的上下文管理器。考虑以下恢复文件的示例:
from contextlib import closing
import shelve as sv
with closing(sv.open('datafile')) as data: # opens a data file
A = data['my_matrix'] # here we used the key
...
搁置对象具有所有字典方法，例如键和值，并且可以像字典一样使用。请注意，只有在调用了close或sync方法后，更改才会写入文件。
读写 Matlab 数据文件
SciPy 能够使用该模块以 Matlab 的.mat文件格式读写数据。命令是loadmat和savemat。要加载数据，请使用以下语法:
import scipy.io
data = scipy.io.loadmat('datafile.mat')
变量数据现在包含一个字典，键对应于保存在.mat文件中的变量名。变量是 NumPy 数组格式。保存到.mat文件包括创建一个包含所有要保存的变量(变量名和值)的字典。接下来的命令是savemat:
data = {}
data['x'] = x
data['y'] = y
scipy.io.savemat('datafile.mat',data)
这将在读入 Matlab 时保存名称相同的 NumPy 数组x和y。
读写图像
SciPy 附带了一些处理图像的基本功能。模块功能将图像读取到 NumPy 数组。该函数将数组保存为图像。下面将把一个 JPEG 图像读入一个数组，打印形状和类型，然后用一个调整大小的图像创建一个新的数组，并将新的图像写入文件:
import scipy.misc as sm
# read image to array
im = sm.imread("test.jpg")
print(im.shape) # (128, 128, 3)
print(im.dtype) # uint8
# resize image
im_small = sm.imresize(im, (64,64))
print(im_small.shape) # (64, 64, 3)
# write result to new image file
sm.imsave("test_small.jpg", im_small)
请注意数据类型。图像几乎总是以范围 0 的像素值存储...255 作为 8 位无符号整数。第三个形状值显示图像有多少颜色通道。在这种情况下， 3 表示它是一个彩色图像，其值按以下顺序存储:红色im[0]、绿色im[1]、蓝色im[2]。灰度图像只有一个通道。
为了处理图像，SciPy 模块scipy.misc包含许多有用的基本图像处理功能，如滤波、变换和测量。
总结
当处理大量数据的测量和其他来源时，文件处理是不可避免的。与其他程序和工具的通信也是通过文件处理来完成的。
你学会了像其他人一样用readlines、write等重要方法将文件视为 Python 对象。我们展示了如何通过特殊属性保护文件，这些属性可能只允许读或写访问。
写入文件的方式通常会影响进程的速度。我们看到了如何通过酸洗或使用shelve方法存储数据。
十三、测试
在这一章中，我们将集中讨论科学编程测试的两个方面。第一个方面是在科学计算中测试什么这个经常很难的话题。第二个方面涉及如何测试的问题。我们将区分手动测试和自动测试。手动测试是每个程序员为了快速检查一个实现是否工作而做的事情。自动化测试是这个想法的精炼的、自动化的变体。我们将介绍一些可用于一般自动测试的工具，着眼于科学计算的特殊情况。
手动测试
在代码开发过程中，为了测试它的功能，您会做很多小测试。这可以称为手动测试。通常，您会通过在交互环境中手动测试函数来测试给定的函数是否完成了它应该做的事情。例如，假设您实现了等分算法。这是一种寻找标量非线性函数的零点(根)的算法。要启动算法，必须给定一个区间，其属性是函数在区间边界上取不同的符号，更多信息请参见练习 4 、第 7 章、函数。
然后，您将测试该算法的实现，通常是通过检查:
- 当函数在区间边界具有相反符号时，找到了一个解
- 当函数在区间边界具有相同的符号时，会引发异常
尽管看起来很有必要，但手动测试并不令人满意。一旦你确信代码做了它应该做的事情，你就可以用相对较少的演示例子来说服其他人相信代码的质量。在那个阶段，人们经常对开发过程中进行的测试失去兴趣，它们被遗忘甚至删除。一旦你改变了一个细节，事情就不再正常工作，你可能会后悔你以前的测试不再可用。
自动测试
开发任何代码的正确方法是使用自动测试。优点是:
- 在每次代码重构之后和任何新版本发布之前，大量测试的自动重复。
- 代码使用的静默文档。
- 代码测试覆盖范围的文档:事情是在变更前工作还是某个方面从未测试过？
注
程序中的变化，特别是不影响其功能的结构变化，叫做代码重构。
我们建议与代码并行开发测试。好的测试设计本身就是一门艺术，很少有投资能像好的测试投资那样保证开发时间节约的良好回报。
现在，我们将考虑自动化测试方法来实现一个简单的算法。
测试等分算法
让我们检查一下等分算法的自动化测试。利用该算法，可以找到实值函数的零点。在第七章、功能中练习 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 {}'.format(tol))
我们假设这存储在bisection.py文件中。作为第一个测试用例，我们测试发现函数 f ( x ) = x 的零点为:
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()线手动运行测试。
有许多工具可以自动化这种调用。
现在让我们设置一个测试，当函数在区间两端具有相同的符号时，该测试检查bisect是否引发异常。现在，我们假设引发的异常是ValueError异常。在以下示例中，我们将检查初始间隔[ a 、 b ]。对于等分算法，它应该满足一个符号条件:
def test_badinput():
try:
bisect(lambda x: x,0.5,1)
except ValueError:
pass
else:
raise AssertionError()
test_badinput()
在这种情况下，如果异常不是ValueError类型，则会引发AssertionError。有一些工具可以简化前面的构造，以检查是否引发了异常。
另一个有用的测试是边缘案例测试。在这里，我们测试参数或用户输入，这可能会产生程序员无法预见的数学上未定义的情况或程序状态。例如，如果两个界限相等会发生什么？如果 a > b 会发生什么？
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 interval bounds failed'
test_equal_boundaries()
test_reverse_boundaries()
使用单元测试包
标准的unittest Python 包极大地方便了自动化测试。这个包要求我们重写测试以兼容。第一个测试必须在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开始执行这些方法。
测试设置和拆卸方法
类unittest.TestCase提供了两个特殊的方法，setUp和tearDown，它们在每次调用测试方法之前和之后运行。测试发电机时需要这样做，每次测试后发电机都会耗尽。我们通过测试一个程序来演示这一点，该程序检查给定字符串首次出现的文件中的行:
class NotFoundError(Exception):
pass
def find_string(file, string):
for i,lines in enumerate(file.readlines()):
if string in lines:
return i
raise NotFoundError(
'String {} not found in File {}'.format(string,file.name))
我们假设该代码保存在find_in_file.py文件中。测试必须准备一个文件，并在测试后打开和删除它，如下例所示:
import unittest
import os # used for, for example, deleting files
from find_in_file import find_string, NotFoundError
class TestFindInFile(unittest.TestCase):
def setUp(self):
file = open('test_file.txt', 'w')
file.write('aha')
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, 'aha')
self.assertEqual(line_no, 0)
def test_not_exists(self):
self.assertRaises(NotFoundError, find_string,
self.file, 'bha')
if __name__=='__main__':
unittest.main()
每次测试前setUp运行，然后tearDown执行。
参数化测试
人们经常想用不同的数据集重复同样的测试。当使用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, "test_{name}".format(name=data['name']),
make_test_function(data))
if __name__=='__main__':
unittest.main()
在本例中，数据以字典列表的形式提供。make_test_function函数动态生成一个测试函数，该函数使用一个特定的数据字典，用之前定义的方法checkifzero执行测试。最后，命令setattr用于制作类Tests的这些测试函数方法。
断言工具
在这一部分，我们收集了最重要的工具来提高一个AssertionError。我们从unittest看到了assert命令和两个工具，即assertAlmostEqual。下表(表 13.1 )总结了最重要的断言工具和相关模块:
表 13.1:Python、unittest 和 NumPy 中的断言工具
浮动比较
两个浮点数不应该与==比较，因为一次计算的结果往往会因为舍入误差而略有偏差。出于测试目的，有许多工具可以测试浮点数的相等性。首先，allclose检查两个数组是否几乎相等。它可以用于测试功能，如图所示:
self.assertTrue(allclose(computed, expected))
这里，self指的是一个unittest.Testcase实例。numpy包testing中也有测试工具。这些是通过使用以下方式导入的:
import numpy.testing
使用numpy.testing.assert_array_allmost_equal或numpy.testing.assert_allclose测试两个标量或两个数组是否相等。如上表所示，这些方法描述所需精度的方式不同。
QR 因式分解将给定矩阵分解为正交矩阵 Q 和上三角矩阵 R 的乘积，如下例所示:
import scipy.linalg as sl
A=rand(10,10)
[Q,R]=sl.qr(A)
方法应用正确吗？我们可以通过验证 Q 确实是一个正交矩阵来检验这一点:
import numpy.testing as npt
npt.assert_allclose(
dot(Q.T,self.Q),identity(Q.shape[0]),atol=1.e-12)
此外，我们可以通过检查 A = QR 来执行健全性测试:
import numpy.testing as npt
npt.assert_allclose(dot(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(
dot(self.Q.T,self.Q),identity(self.Q.shape[0]),
atol=1.e-12)
def test_sanity(self):
npt.assert_allclose(dot(self.Q,self.R),self.A)
if __name__=='__main__':
unittest.main()
注意在assert_allclose中，参数atol默认为零，当处理具有小元素的矩阵时，这通常会导致问题。
单元和功能测试
到目前为止，我们只使用了功能测试。功能测试检查功能是否正确。对于二等分算法，这种算法实际上在有一个零的时候会找到一个零。在这个简单的例子中，单元测试是什么并不清楚。虽然这看起来有点做作，但是仍然可以对二等分算法进行单元测试。它将展示单元测试如何经常导致更加条块分割的实现。
因此，在二等分法中，我们想要检查，例如，在每一步中，间隔的选择是否正确。怎么做？请注意，用当前的实现是绝对不可能的，因为算法隐藏在函数内部。一种可能的补救方法是只运行二等分算法的一个步骤。由于所有的步骤都是相似的，我们可能会说我们已经测试了所有可能的步骤。我们还需要能够在算法的当前步骤检查当前界限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.)
调试
在测试时，调试有时是必要的，特别是如果不能立即弄清楚为什么给定的测试没有通过。在这种情况下，能够在交互会话中调试给定的测试是非常有用的。然而，unittest.TestCase类的设计使这变得困难，这阻碍了测试用例对象的简单实例化。解决方案是创建一个仅用于调试目的的特殊实例。
假设，在上面TestIdentity类的例子中，我们想要测试test_functionality方法。这将通过以下方式实现:
test_case = TestIdentity(methodName='test_functionality')
现在，该测试可以通过以下方式单独运行:
test_case.debug()
这将运行这个单独的测试，并允许调试。
测试发现
如果您编写一个 Python 包，各种测试可能会在包中展开。discover模块查找、导入并运行这些测试用例。命令行的基本调用是:
python -m unittest discover
它开始在当前目录中寻找测试用例，并向下递归目录树来寻找名称中包含'test'字符串的 Python 对象。该命令接受可选参数。最重要的是-s修改开始目录和-p定义识别测试的模式:
python -m unittest discover -s '.' -p 'Test*.py'
测量执行时间
为了在代码优化上做出决定，人们通常必须比较几种代码选择，并根据执行时间来决定应该优先选择哪种代码。此外，在比较不同算法时，讨论执行时间是一个问题。在这一节中，我们提供了一种简单易行的方法来测量执行时间。
带魔法功能的计时
测量单个语句执行时间最简单的方法是使用 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参数控制重复次数。
用 Python 模块计时
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
与前面例子中的方法相反，我们获得了所有获得的测量值的列表。由于计算时间可能会根据计算机的总负载而变化，因此可以将此列表中的最小值视为执行语句所需计算时间的良好近似值。
使用上下文管理器计时
最后，我们提出了第三种方法。它展示了上下文管理器的另一个应用。我们首先构建一个上下文管理器对象来测量经过的时间，如图所示:
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('Time elapsed {} seconds'.format(self.elapsed))
return False
回想一下_ _enter_ _和_ _exit_ _方法使这个类成为一个上下文管理器。_ _exit_ _方法的参数ty、val和tb在正常情况下None。如果在执行过程中引发异常，它们会获取异常类型、异常值和追溯信息。return 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
总结
没有测试就没有程序开发！我们展示了良好组织和记录测试的重要性。一些专业人员甚至通过首先指定测试来开始开发。自动测试的一个有用工具是模块unittest，我们已经详细解释过了。虽然测试提高了代码的可靠性，但是需要分析来提高性能。替代的编码方式可能会导致很大的性能差异。我们展示了如何测量计算时间以及如何定位代码中的瓶颈。
练习
Ex。1 →两个矩阵 A 、 B 叫做相似，如果存在一个矩阵 S ，这样 B = S -1 A S 。 A 和 B 具有相同的特征值。写一个测试，通过比较两个矩阵的特征值来检查它们是否相似。是功能测试还是单元测试？
Ex。2 →创建两个高维向量。比较各种方式计算dot乘积的执行时间:
- SciPy 功能:
dot(v,w) - 生成器和总和:
sum((x*y for x,y in zip(v,w))) - 综合清单及合计:
sum([x*y for x,y in zip(v,w)])
Ex。3 →让 u 成为向量。矢量 v 带分量
叫做 u 的移动平均线。确定计算 v 的两种选择中哪一种更快:
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)])
十四、综合示例
在这一章中，我们提出了一些全面和更长的例子，并简要介绍了理论背景和它们的完整实现。通过这一点，我们想向你展示如何在实践中使用本书中定义的概念。
多项式
首先，我们将通过设计一个多项式类来展示 Python 构造的强大功能。我们将给出一些理论背景，这将引导我们得到一个需求列表，然后我们将给出代码，并给出一些注释。
注意，这个类在概念上不同于类numpy.poly1d。
理论背景
一个多项式:p(x)= anxn+an-1xn-1T13】+…+aT16】1T18】x+aT20】0 由其度、其表示及其系数定义。上式所示的多项式表示称为单项表示。在该表示中，多项式被写成单项式的线性组合，xIT25。或者，多项式可以写成:
-
Newton representation with the coefficients ci and n points, x0, …, xn-1:
p(x)= c**+【c】*x【x】-我...。+c【n】t30(x0)(*xn-)
-
Lagrange representation with the coefficients *yi*and n+1 points, x0, … , xn:
p(x)=yT6】0T8】lT10】0(x)+yT16】1T18】lT20】1(x)+……+ynln
具有基本功能:
有无限多种表现形式，但我们在这里仅限于这三种典型的表现形式。
多项式可以由插值条件确定:
p(xIT5)=yT8】IT10】I= 0，…， n
给定不同的值 x i 和任意值 y i 作为输入。在拉格朗日公式中，插值多项式是直接可用的，因为它的系数是插值数据。牛顿表示中的插值多项式的系数可以通过一个递归公式获得，称为除差公式:
c i ，0 = y i、 和
。
最后，设置。
单项表示中的插值多项式的系数是通过求解一个线性系统获得的:
以给定多项式 p (或其倍数)作为特征多项式的矩阵称为伴随矩阵。伴随矩阵的特征值是多项式的零点(根)。计算 p 零点的算法可以通过首先建立其伴随矩阵，然后用eig计算特征值来构造。牛顿表示的多项式的伴随矩阵如下:
任务
我们现在可以制定一些编程任务:
- 用
points、degree、coeff和basis属性写一个名为PolyNomial的类，其中:points是元组列表( x i ，yIT6)degree是对应插值多项式的次数coeff包含多项式系数basis是一个字符串，说明使用了哪种表示法
- 为类提供一种在给定点计算多项式的方法。
- 为该类提供一个名为
plot的方法，在给定的时间间隔内绘制多项式。 - 写一个名为
__add__的方法，返回两个多项式之和的多项式。请注意，只有在单项情况下，总和才可以通过对系数求和来计算。 - 写一个计算以单项形式表示的多项式系数的方法。
- 写一个计算多项式伴随矩阵的方法。
- 写一个通过计算伴随矩阵的特征值来计算多项式零点的方法。
- 写一个计算多项式的方法，该多项式是给定多项式的第I次导数。
- 写一个检查两个多项式是否相等的方法。可以通过比较所有系数来检查相等性(零前导系数应该无关紧要)。
多项式类
现在让我们基于多项式的单项公式设计一个多项式基类。多项式可以通过给出其关于单项式基的系数或给出插值点列表来初始化，如下所示:
import scipy.linalg as sl
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__方法使用了第 7 章、函数中参数和参数一节中讨论的**args构造。如果没有给定参数，则假设多项式为零。如果多项式由插值点给出，则通过求解范德蒙系统来计算系数的方法如下:
def point_2_coeff(self):
return sl.solve(vander(self.x),self.y)
如果给定了 k 系数，则 k 插值点通过以下方式构造:
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 章、类中特殊方法一节的例子。)这里，该方法使用命令polyval。下一步，为了方便起见，我们只是增加了两种方法，我们用property装饰器来装饰(参见第 7 章、功能中作为装饰器的功能部分】 :
@property
def x(self):
return self.points[:,0]
@property
def y(self):
return self.points[:,1]
让我们解释一下这是怎么回事。我们定义了一种提取数据的 x 值的方法，用于定义多项式。类似地，定义了提取数据的 y 值的方法。使用property装饰器，调用方法的结果就好像它只是多项式的一个属性一样。有两种编码选择:
-
We use a method call:
def x(self): return self.interppoints[:,0]这允许通过调用
p.x()访问 x 值。 -
我们使用
property装饰机。它允许我们通过以下语句访问 x 值:p.x
我们选择第二种变体。定义一个__repr__方法(参考第八章类中属性一节)总是一个好的做法。至少对于结果的快速检查，这种方法是有用的:
def __repr__(self):
txt = 'Polynomial of degree {degree} \n'
txt += 'with coefficients {coeff} \n in {base} basis.'
return txt.format(coeff=self.coeff, degree=self.degree,
base=self.base)
我们现在提供一种绘制多项式的方法，如下所示:
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)
plot(x,y)
xlabel('$x$')
ylabel('$p(x)$')
if plotinterp:
plot(self.x, self.y, 'ro')
注意vectorize命令的使用(参见第 4 章、线性代数-数组中函数作用于数组一节。__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)
为此，功能eigvals必须首先从scipy.linalg导入。我们来举一些用法的例子。
首先，我们从给定的插值点创建一个多项式实例:
p = PolyNomial(points=[(1,0),(2,3),(3,8)])
多项式关于单项式基的系数可作为p的属性获得:
p.coeff # returns array([ 1., 0., -1.])
这对应于多项式。由
p.plot(-3.5,3.5)得到的多项式默认图，结果如下图(图 14.1 ):
图 14.1:多项式作图法的结果
最后，我们计算多项式的零点，在这种情况下是两个实数:
pz = p.zeros() # returns array([-1.+0.j, 1.+0.j])
可以通过在以下几点评估多项式来验证结果:
p(pz) # returns array([0.+0.j, 0.+0.j])
牛顿多项式
NewtonPolyNomial类定义了一个关于牛顿基描述的多项式。我们通过使用super命令(参见第 8 章、类中的子类化和继承一节)，让它从多项式基类继承一些常见的方法，例如polynomial.plot、polynomial.zeros，甚至__init__方法的部分内容:
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()))
这里，我们使用除差来计算多项式的牛顿表示，它在这里被编程为生成器:
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]
让我们简单检查一下这是如何工作的:
pts = array([[0.,0],[.5,1],[1.,0],[2,0.]]) # here we define the
interpolation data: (x,y) pairs
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.
多项式求值不同于基类的相应方法。Newton.PolyNomial.__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 dot(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
注意布尔数组的使用。这些练习将进一步建立在这个基础上。
光谱聚类
特征向量的一个有趣的应用是对数据进行聚类。使用从距离矩阵导出的矩阵的特征向量，可以将未标记的数据分成组。谱聚类方法的名字来源于这个矩阵的谱的使用。 n 元素的距离矩阵(例如，数据点之间的成对距离)是 n × n 对称矩阵。给定这样一个 n × n 距离矩阵 M 和距离值MijT7】，我们可以如下创建数据点的拉普拉斯矩阵:
这里，I 是单位矩阵， D 是包含 M 行和的对角矩阵，
数据聚类是从 L 的特征向量中获得的。在只有两类数据点的最简单情况下，第一个特征向量(即对应于最大特征值的特征向量)通常足以分离数据。
下面是一个简单的两类聚类的例子。下面的代码创建一些 2D 数据点，并基于拉普拉斯矩阵的第一特征向量对它们进行聚类:
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
对应于最大特征值的特征向量给出了分组(例如，通过在 0 处设置阈值)，并且可以用以下公式表示:
largest=abs(S).argmax()
plot(V[:,largest])
下图(图 14.2 )显示了简单两类数据集的光谱聚类结果:
图 14.2:显示了简单的两类聚类的结果
对于更难的数据集和更多的类，通常取 k 最大特征值对应的 k 特征向量，然后用其他方法对数据进行聚类，但使用特征向量代替原始数据点。一个常见的选择是k-意味着聚类算法，这是下一个例子的主题:
特征向量用作 k 的输入-意味着聚类，如下所示:
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返回一个包含三个数组的列表，左右奇异向量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')
下图显示了简单多类数据集的光谱聚类结果:
图 14.3:简单四类数据集的光谱聚类示例。
求解初值问题
在这一节中，我们将考虑对给定初始值数值求解一个常微分方程组的数学任务:
和 ( t ) = f ( t，y ) 和 ( t 0】
这个问题的解决方案是一个函数 y 。一种数值方法旨在计算好的近似值，yIT5】∑y*(tI)在离散点，即通信点 t i ，在感兴趣的区间内[ t 0 ，t e 。我们在一个类中收集描述问题的数据，如下所示:*
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
微分方程:
描述一个数学钟摆； y 1 描述其相对于垂直轴的角度， g 为引力常数， l 为其长度。初始角度为π/2，初始角速度为零。
钟摆问题成为问题类的一个实例，如下所示:
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.
有了这些，我们可以创建实例来获得摆式 ODE 的相应离散化版本:
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')
图 14.4:用显式欧拉方法模拟摆，并与更精确的龙格-库塔 4 方法的结果进行比较
讨论替代的班级设计是值得的。什么应该放在单独的班级里，什么应该捆绑到同一个班级里？
- 我们把数学问题和数值方法严格分开了。初始值应该放在哪里？他们应该是问题的一部分还是解决者的一部分？还是应该将它们作为求解器实例的求解方法的输入参数？人们甚至可以设计程序，使其允许多种可能性。使用其中一种替代品的决定取决于该程序的未来用途。通过将初始值作为求解方法的输入参数，可以简化参数识别中对各种初始值的循环。另一方面，用相同的初始值模拟不同的模型变量会促使将初始值与问题联系起来。
- 为了简单起见，我们只介绍了具有恒定给定步长的解算器。
IVPsolver类的设计是否适合自适应方法的未来扩展，在这种扩展中给出的是公差而不是步长？ - 我们之前建议使用发电机结构作为步进机构。适应性方法需要不时地拒绝步骤。这种需求是否与
IVPsolver.onestepper中步进机构的设计相冲突？ - 我们鼓励您检查用于求解初始值的两个 SciPy 工具的设计，即
scipy.integrate.ode和scipy.integrate.odeint。
总结
我们在这本书里解释的大部分内容被捆绑到本章的三个较长的例子中。它们模仿代码开发并给出原型，鼓励你改变和面对自己的想法。
您看到科学计算中的代码可以有自己的味道，因为它与数学定义的算法有很强的关系，并且保持代码和公式之间的关系可见通常是明智的。正如您所看到的，Python 在这方面有技巧。
练习
Ex。1 →实现一种方法__add__，通过将两个给定的多项式 p 和 q 相加，构造一个新的多项式 p+q 。在单项式中，多项式只需将系数相加即可，而在牛顿式中，系数取决于插值点的横坐标 x i 。在将两个多项式的系数相加之前，多项式 q 必须获得新的插值点，其横坐标 x i 与 p 的横坐标一致，并且必须为此提供方法__changepoints__。它应该改变插值点并返回一组新的系数。
Ex。2 →编写转换方法，将一个多项式从牛顿形式转换为单项式形式，反之亦然。
Ex。3 →编写一个名为add_point的方法，该方法将一个多项式 q 和一个元组 (x，y) 作为参数，并返回一个新的插值self.points和 (x，y) 的多项式。
Ex。4 →编写一个名为LagrangePolynomial的类，以拉格朗日形式实现多项式，并尽可能从多项式基类继承。
Ex。5 →为多项式类编写测试。
