树莓派超算和科学计算教程-二-

布客飞龙
135 阅读45分钟

树莓派超算和科学计算教程（二）

原文：Raspberry Pi Supercomputing and Scientific Programming

协议：CC BY-NC-SA 4.0

八、Python 3 中的并行编程

在上一章中，我们学习了如何对各种型号的 Raspberry Pi 进行超频，以增加其计算能力。在本章中，我们将学习如何用 Python 和 MPI4PY 编写并行程序。我更喜欢 Python，因为它简单，而且 Python 中的代码不那么吓人。我们将探索 MPI 概念，并用 MPI4PY 在 Python 中实现这些概念。

我们将研究和实现的 MPI 概念如下:

  • MPI 等级和流程
  • 发送和接收数据
  • 数据标记
  • 广播数据
  • 分散和收集数据

MPI4PY 的基础知识

在本书的前面部分，我们学习了一些 MPI 概念。这一章我们再研究几个。

MPI 使用单程序多数据( SPMD )的概念。以下是 SPMD 建筑的要点:

  • 所有进程(称为等级)运行相同的代码，每个进程访问不同的数据部分。
  • 所有进程同时启动。

并行程序被分解成单独的进程，称为等级。每个列都有自己的地址空间，这需要跨列划分数据。每个等级在自己的私有内存中保存一部分程序数据。等级按从 0 到 n-1 的顺序编号。下图(图 8-1 )描述了同时运行的多个队列。

A447085_1_En_8_Fig1_HTML.gif

图 8-1。

Multiple ranks running simeltaneously

在下一节中，我们将看到一个带有等级的基本程序。

MPI4PY 入门

让我们从简单的 Hello World 开始吧！用 MPI4PY 用 Python 编程。

from mpi4py import MPI
import sys

comm = MPI.COMM_WORLD
name = MPI.Get_processor_name()

sys.stdout.write("Hello World!")
sys.stdout.write(" Name: %s, My rank is %d\n" % (name, comm.rank))

Listing 8-1.prog01.py

在上面的代码(清单 8-1 )中，语句from mpi4py import MPI导入所需的 MPI4PY 库。在第六章中，我们学习了 MPI 中通信者的概念。MPI.COMM_WORLD是传播者。它用于运行在集群进程上的进程之间的所有 MPI 通信。Get_processor_name()返回当前进程运行的主机名。comm.rank是当前流程的等级。下图(图 8-2 )描绘了COMM_WORLD。

A447085_1_En_8_Fig2_HTML.gif

图 8-2。

COMM_WORLD

你可能已经注意到我们正在控制台上使用sys.stdout.write()进行打印。这是因为我希望代码兼容 Python 编程语言的两种解释器，python(Python 2 的解释器)和python3。在本书中，我们不会使用任何特定于任一解释器的特性或编程结构。因此，代码可以使用两种解释器运行。

这一章我们已经开始编码了，后面的章节有很多代码样例和练习。将代码和数据组织在单独的目录中是一个好主意。在 lxterminal 中逐一运行以下命令:

mpirun -hostfile myhostfile -np 4 mkdir /home/pi/book
mpirun -hostfile myhostfile -np 4 mkdir /home/pi/book/code
mpirun -hostfile myhostfile -np 4 mkdir /home/pi/book/code/chapter08

这将在小型超级计算机的所有节点上创建相同的目录结构。现在将上面的代码保存在∼/book/code/chapter08目录下一个名为prog01.py的文件中。使用scp as将代码文件复制到所有节点上的那个目录，如下所示:

scp book/code/chapter08/prog01.py 192.168.0.2:/home/pi/book/code/chapter08/
scp book/code/chapter08/prog01.py 192.168.0.3:/home/pi/book/code/chapter08/
scp book/code/chapter08/prog01.py 192.168.0.4:/home/pi/book/code/chapter08/

最后，在pi001上运行mpirun,如下所示:

mpirun -hostfile myhostfile -np 4 python3 ∼/book/code/chapter08/prog01.py

以下是输出:

Hello World! Name: pi001, My rank is 0
Hello World! Name: pi002, My rank is 1
Hello World! Name: pi004, My rank is 3
Hello World! Name: pi003, My rank is 2

对于我们将在本章剩余部分讨论的所有其他代码示例，我们必须遵循相同的步骤。让我简单重复一遍:在chapter08目录下创建一个 Python 代码文件，将该文件复制到集群所有节点的chapter08目录下，最后使用mpirun和 Python 解释器来执行代码。

条件语句

我们可以在 MPI4PY 代码中使用如下条件语句(清单 8-2 ):

from mpi4py import MPI
import sys

comm = MPI.COMM_WORLD
sys.stdout.write("My rank is: %d\n" % (comm.rank))

if comm.rank == 0:
    sys.stdout.write("Doing the task of Rank 0\n")

if comm.rank == 1:
    sys.stdout.write("Doing the task of Rank 1\n")

Listing 8-2.prog02.py

在这段代码中，我们检查进程等级是 0 还是 1，然后将更多消息打印到控制台。用mpiexec运行程序如下:

mpirun -hostfile myhostfile -np 4 python3 ∼/book/code/chapter08/prog02.py

上面的程序(清单 8-2 )的输出如下:

My rank is: 0
Doing the task of Rank 0
My rank is: 1
Doing the task of Rank 1
My rank is: 3
My rank is: 2

检查流程的数量

让我们编写代码(清单 8-3 )来显示 MPI 进程的等级和数量。

from mpi4py import MPI
import sys

comm = MPI.COMM_WORLD
rank = comm.rank
size = comm.size

sys.stdout.write("Rank: %d," % rank)
sys.stdout.write(" Process Count: %d\n" % size)

Listing 8-3.prog03.py

在上面的代码中，comm.size给出了在集群中运行的 MPI 进程的数量。使用mpiexec运行上面的代码，如下所示:

mpirun -hostfile myhostfile -np 4 python3 ∼/book/code/chapter08/prog03.py

输出如下所示:

Rank: 0, Process Count: 4
Rank: 1, Process Count: 4
Rank: 2, Process Count: 4
Rank: 3, Process Count: 4

发送和接收数据

使用send()和receive()进行进程间的数据传输是进程间最简单的通信形式。我们可以用这个实现一对一的交流。下图(图 8-3 )清楚地解释了这一点。

A447085_1_En_8_Fig3_HTML.gif

图 8-3。

One-to-one communication

让我们来看看代码示例(清单 8-4 )的相同之处。

from mpi4py import MPI
import time
import sys

comm = MPI.COMM_WORLD

rank = comm.rank
size = comm.size
name = MPI.Get_processor_name()

shared = 3.14

if rank == 0:
    data = shared
    comm.send(data, dest=1)
    comm.send(data, dest=2)
    sys.stdout.write("From rank %s, we sent %f\n" % (name, data))
    time.sleep(5)

elif rank == 1:
    data = comm.recv(source=0)
    sys.stdout.write("On rank %s, we received %f\n" % (name, data))

elif rank == 2:
    data = comm.recv(source=0)
    sys.stdout.write("On rank %s, we received %f\n" % (name, data))

Listing 8-4.prog04.py

在上面的代码示例中，我们从等级为 0 的流程中发送数据。等级为 1 和 2 的进程正在接收数据。

让我们运行程序。

mpirun -hostfile myhostfile -np 4 python3 ∼/book/code/chapter08/prog04.py

上面的程序(清单 8-4 )的输出如下:

On rank pi002, we received 3.140000
On rank pi003, we received 3.140000
From rank pi001, we sent 3.140000

动态发送和接收数据

到目前为止，我们已经为发送和接收数据的过程编写了条件语句。然而，在大型分布式网络中，由于进程数量的不断变化，这种类型的数据传输并不总是可能的。此外，用户可能不想手工编写条件语句。

下面的例子(清单 8-5 )展示了动态数据传输的概念。

from mpi4py import MPI
import sys

comm = MPI.COMM_WORLD
rank = comm.rank
size = comm.size
name = MPI.Get_processor_name()

shared = (rank+1)*(rank+1)

comm.send(shared, dest=(rank+1) % size)
data = comm.recv(source=(rank-1) % size)

sys.stdout.write("Name: %s\n" % name)
sys.stdout.write("Rank: %d\n" % rank)
sys.stdout.write("Data %d came from rank: %d\n" % (data, (rank-1) % size))

Listing 8-5.prog05.py

在上面的代码(清单 8-5 )中，每个进程都从前面的进程接收数据。这种情况一直持续到结束并循环，这样第一个进程从最后一个进程接收数据。

让我们运行代码。

mpirun -hostfile myhostfile -np 4 python3 ∼/book/code/chapter08/prog05.py

代码的输出如下所示:

Name: pi001
Rank: 0
Data 16 came from rank: 3
Name: pi002
Rank: 1
Data 1 came from rank: 0
Name: pi003
Rank: 2
Data 4 came from rank: 1
Name: pi004
Rank: 3
Data 9 came from rank: 2

如前所述，等级为 0 的进程(第一个进程)从等级为 3 的进程(最后一个进程)接收数据。

数据标记

在前面的例子(清单 8-5 )中，我们研究了如何用 MPI 发送和接收数据。这给好奇的程序员提出了一个基本问题:我们如何在进程间交换多个数据项？我们可以将多个数据项从一个进程发送到另一个进程。然而，在接收端，我们会遇到区分不同数据项的问题。解决这个问题的方法是标记。看看下面的代码示例(清单 8-6 )。

from mpi4py import MPI
import sys

comm = MPI.COMM_WORLD
rank = comm.rank
size = comm.size
name = MPI.Get_processor_name()

if rank == 0:
    shared1 = {'d1': 55, 'd2': 42}
    comm.send(shared1, dest=1, tag=1)

    shared2 = {'d3': 25, 'd4': 22}
    comm.send(shared2, dest=1, tag=2)

if rank == 1:
    receive1 = comm.recv(source=0, tag=1)
    sys.stdout.write("d1: %d, d2: %d\n" % (receive1['d1'], receive1['d2']))
    receive2 = comm.recv(source=0, tag=2)
    sys.stdout.write("d3: %d, d4: %d\n" % (receive2['d3'], receive2['d4']))

Listing 8-6.prog06.py

在上面的例子中，我们将两个不同的字典shared1和shared2从等级为 0 的进程发送到等级为 1 的进程。在源端，shared1被标记为 1，shared2被标记为 2。在目的地，我们可以从与它们相关联的标签中区分不同的数据项。

用下面的命令运行上面的代码(清单 8-6 ):

mpirun -hostfile myhostfile -np 4 python3 ∼/book/code/chapter08/prog06.py

输出如下所示:

d1: 55, d2: 42
d3: 25, d4: 22

数据标记使程序员能够更好地控制数据流。当多个数据在流程之间交换时，数据标记是必须的。

数据广播

当数据从一个进程发送到所有其他进程时，就称为广播。考虑下面的代码(清单 8-7 ):

from mpi4py import MPI
import sys

comm = MPI.COMM_WORLD
rank = comm.rank

if rank == 0:
    data = {'a': 1, 'b': 2, 'c': 3}
else:
    data = None

data = comm.bcast(data, root=0)
sys.stdout.write("Rank: %d, Data: %d, %d, %d\n"
                 % (rank, data['a'], data['b'], data['c']))

Listing 8-7.prog07.py

在上面的代码(清单 8-7 )中，在if语句中，只有当进程的等级为 0 时，我们才给数据分配一个字典。bcast()向所有进程广播数据。

运行程序。

mpirun -hostfile myhostfile -np 4 python3 ∼/book/code/chapter08/prog07.py

输出如下所示:

Rank: 0, Data: 1, 2, 3
Rank: 1, Data: 1, 2, 3
Rank: 2, Data: 1, 2, 3
Rank: 3, Data: 1, 2, 3

数据分散

在广播中，我们向所有进程发送相同的数据。在分散中，我们将不同的数据块发送给所有进程。例如，我们有一个包含四个项目的列表。在广播中，我们将所有这四个项目发送给所有进程，而在分散中，我们将列表中的项目发送给进程，这样每个进程都从列表中接收一个项目。下面的程序(清单 8-8 )演示了这一点。

from mpi4py import MPI
import sys

comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()

if rank == 0:
    data = [x for x in range(0,size)]
    sys.stdout.write("We will be scattering: ")
    sys.stdout.write(" ".join(str(x) for x in data))
    sys.stdout.write("\n")
else:
    data = None

data = comm.scatter(data, root=0)
sys.stdout.write("Rank: %d has data: %d\n" % (rank, data))

Listing 8-8.prog08.py

在上面的代码(清单 8-8 )中，我们创建了一个名为 data 的列表，它的元素数量等于集群中的进程数。scatter()用于将数据分散到所有进程中。

运行代码。

mpirun -hostfile myhostfile -np 4 python3 ∼/book/code/chapter08/prog08.py

以下是输出:

Rank: 1 has data: 1
We will be scattering: 0 1 2 3
Rank: 0 has data: 0
Rank: 2 has data: 2
Rank: 3 has data: 3

正如我们所看到的，每个进程从列表中接收一个项目。scatter()的限制是我们分散的数据列表的大小不能超过进程的数量。

数据采集

收集数据的想法与分散相反。主进程收集由其他进程处理的所有数据。

下面的程序(列表 8-9 )演示了gather()方法。

from mpi4py import MPI
import sys

comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()

if rank == 0:
    data = [x for x in range(0,size)]
    sys.stdout.write("We will be scattering: ")
    sys.stdout.write(" ".join(str(x) for x in data))
    sys.stdout.write("\n")
else:
    data = None

data = comm.scatter(data, root=0)
sys.stdout.write("Rank: %d has data: %d\n" % (rank, data))
data *= data

newData = comm.gather(data, root=0)

if rank == 0:
    sys.stdout.write("We have gathered: ")
    sys.stdout.write(" ".join(str(x) for x in newData))
    sys.stdout.write("\n")

Listing 8-9.prog09.py

在上面的程序(列表 8-9 )中，主进程分散数字列表。所有 MPI 进程都从列表中接收一个元素(列表的大小等于 MPI 进程的数量)。每个进程执行它接收的元素的一个操作。在我们的例子中，它是数字平方的计算。然而，在现实世界的超级计算中，操作可能相当复杂。

一旦操作完成，主进程将所有已处理的元素收集到一个新列表中。

运行代码。

mpirun -hostfile myhostfile -np 4 python3 ∼/book/code/chapter08/prog09.py

输出如下所示:

Rank: 1 has data: 1
Rank: 3 has data: 3
We will be scattering: 0 1 2 3
Rank: 0 has data: 0
We have gathered: 0 1 4 9
Rank: 2 has data: 2

结论

在本章中，我们介绍了 Python 的 MPI4PY 库。我们用 MPI4PY 学习并实验了并行编程中各种有趣的概念。在下一章中，我们将开始使用 Python 3 中带有 Raspberry Pi 的 SciPy 栈。

九、SciPy 栈和符号编程简介

在上一章中，我们学习了如何使用我们为 MPI4PY 和 Python 3 并行编程而构建的 Raspberry Pi 集群。在本章中，我们将介绍 SciPy 栈并将其安装在 Pi 上。我们还将开始使用 SymPy 进行符号编程。

科学 Python 栈

SciPy(Scientific Python 的缩写)是一个用 Python 编写的科学和技术计算的开源库。

SciPy 具有用于数值运算、常微分方程解算器、快速傅立叶变换、最优化、线性代数、积分、插值、信号处理和图像处理的模块。SciPy 被全世界的科学、数学和工程社区广泛使用。还有许多其他库使用 SciPy 和 NumPy 的核心模块进行操作。OpenCV 和 SciKit 是其他主要库使用 NumPy 和/或 SciPy 的最好例子。

SciPy 栈包含以下组件:

  • NumPy 是一个用于数值计算的库。它提供了数字和科学计算所需的所有基本数据类型。
  • SciPy 库有许多用于科学编程的模块。
  • Matplotlib 用于数据可视化。
  • SymPy 用于符号编程。
  • IPython 是一个高级 Python 解释器，增加了一些特性。
  • Pandas 用于数据分析。
  • Nose 用于自动化测试用例。

下图(图 9-1 )总结了 Python SciPy 栈在科学计算世界中的作用。

A447085_1_En_9_Fig1_HTML.gif

图 9-1。

The components of the SciPy stack

在本书中，我们将学习 NymPy、SciPy 库、Matplotlib 和 SymPy。

SciPy 栈的安装

在 Raspberry Pi 上安装 SciPy 栈的最佳方式是使用apt-get和pip。

首先，使用以下命令升级pip:

sudo python3 -m pip install --upgrade pip

使用以下命令安装 SymPy:

sudo pip3 install sympy

SciPy 栈的其余组件可以使用apt-get实用程序方便地安装，如下所示:

sudo apt-get install python3-matplotlib -y
sudo apt-get install python3-scipy -y
sudo apt-get install python3-numpy -y

这将安装所需的 SciPy 栈组件。

-很好

SymPy 网站上说:

SymPy is a Python library of symbolic mathematics. Its goal is to become a fully functional computer algebra system (CAS), while keeping the code as simple as possible so as to be easy to understand and expand. SymPy is written entirely in Python.

SymPy 是 BSD 许可的，是免费的。它完全是用 Python 编写的。它依赖于mpmath并且本质上是轻量级的，因为没有其他依赖项。

入门指南

让我们从 SymPy 开始。运行以下命令创建并导航到该章节的目录:

cd ∼
cd book
cd code
mkdir chapter09
cd chapter09

我们将把这一章的代码保存在目录chapter09中，并将在本书的其余部分继续这一做法，即按章创建目录来组织代码。

我假设您有一些数学和微积分的基础知识，所以我在解释代码时不必解释这些基础知识。让我们看一个符号计算概念的简单例子(列表 9-1 )。

import math
import sympy
print(math.sqrt(9))
print(math.sqrt(8))
print(sympy.sqrt(9))
print(sympy.sqrt(8))
Listing 9-1.prog01.py

运行上面的代码(清单 9-1 )。以下是输出:

3.0
2.8284271247461903
3
2*sqrt(2)

在上面的输出中，我们可以看到math.sqrt()方法直接产生数字格式的结果，而sympy.sqrt()方法只在结果是整数的情况下才产生数字格式的结果。它不会产生一个小数值，而是保持sqrt(2)不变。这样，我们就可以象征性地计算很多数学方程。有了符号数学概念的这些知识，让我们用 Python 3 更深入地研究 SymPy。

标志

我们来研究一下符号的概念。符号类似于数学中的变量。我们可以用它们来计算方程式和表达式。它们也可以用来解方程。sympy.symbols()方法将一个字符串转换成符号变量，如下所示(清单 9-2 ):

from sympy import *

x = symbols('x')
sigma, y = symbols('sigma y')
print(x, y, sigma)

Listing 9-2.prog02.py

输出如下所示:

x y sigma

上面的代码(清单 9-2 )展示了symbols()方法可以接受一个字符串作为参数，其中的标记由空格分隔。

让我们再看一个例子(清单 9-3 )，它演示了用符号对表达式求值。

from sympy import *

x = symbols('x')
expr = x + 1
print(expr.subs(x, 3))

Listing 9-3.prog03.py

这将在表达式中用 3 代替x,并对其求值。代码(清单 9-3 )产生4作为输出。

我们可以替换多个符号如下(列表 9-4 ):

from sympy import *

x, y = symbols('x y')
expr = x + y
print(expr.subs({x:3, y:2}))

Listing 9-4.prog04.py

这里，我们一次替换多个符号。运行代码并检查输出。

将字符串转换为 SymPy 表达式

我们可以将字符串转换成 SymPy 表达式。就像在 Python 中一样，我们可以对指数使用**操作符。下面的代码(清单 9-5 )显示了这一点:

from sympy import *

str = "x**3 + 4*y - 1/5"
expr = sympify(str)
print(expr)

Listing 9-5.prog05.py

方法将一个字符串转换成一个 SymPy 表达式。

我们还可以使用evalf()方法将表达式求值为浮点数。它的默认精度是小数点后 15 位数字；但是，我们可以将精度作为参数传递。下面(列表 9-6 显示了evalf()方法的用例示例:

from sympy import *

expr = sqrt(10)
print(expr.evalf())

print(pi.evalf(20))

x = symbols('x')
expr = sin(2*x)
print(expr.evalf(subs={x: 2.4}))
The output is as follows,
3.16227766016838
3.1415926535897932385
-0.996164608835841
Printing in SymPy

Listing 9-6.prog06.py

Sympy 的打印功能

SymPy 有很多打印机。在任何环境中，在命令提示符下使用init_session()方法都会启动一个交互式会话。下面是一个交互式会话的示例。我在控制台中输入的命令以粗体突出显示。

pi@pi001:∼/book/code/chapter09 $ python3

Python 3.4.2 (default, Oct 19 2014, 13:31:11)
[GCC 4.9.1] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from sympy import *

>>> init_session()

Python console for SymPy 1.0 (Python 3.4.2-32-bit) (ground types: python)

这些命令被执行:

>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.0/

>>> Integral(sqrt(1/x), x)

⌠
⎮        ___
⎮      ╱ 1
⎮    ╱   ─  dx
⎮╲╱     x
⌡
>>> sqrt(x)

√x
>>> (sqrt(x) + sqrt(y))**2

         2
(√x + √y)
>>>

这就是我们如何在交互式控制台中以良好的格式打印表达式。

症状的简化

我们可以使用simplify()方法尽可能简化数学表达式。这下涵盖了大量的表达式。以下(清单 9-8 )就是一个例子:

from sympy import *
x = symbols('x')
print(simplify(sin(x)**2 + cos(x)**2))
print(simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1)))
print(simplify(gamma(x)/gamma(x - 2)))
Listing 9-8.prog08.py

简化的输出如下:

1
x - 1
(x - 2)*(x - 1)

SymPy 中有更多的简化方法。我们可以使用expand()来展开一个多项式表达式，如下所示(列表 9-9 )。

from sympy import *
x, y = symbols('x y')
print(expand((x + y)**2))
print(expand((x + 3)*(y + 5)))
Listing 9-9.prog09.py

以下是扩展后的输出:

x**2 + 2*x*y + y**2
x*y + 5*x + 3*y + 15

类似地，我们可以使用factor()方法(列表 9-10 )来寻找多项式的不可约因子。

from sympy import *
x = symbols('x')
print(factor(x**3 - x**2 + x))
Listing 9-10.prog10.py

输出如下所示:

x*(x**2 - x + 1)

结石

我们甚至可以用 SymPy 来做微积分。我们可以使用diff()方法计算导数如下(列表 9-11 ):

from sympy import *
x = symbols('x')
print(diff(x**3 - x**2 + x, x))
print(diff(x**5, x))
print(diff(sin(x), x))
Listing 9-11.prog11.py

输出是:

3*x**2 - 2*x + 1
5*x**4
cos(x)

我们还可以找到如下表达式的高阶导数(列表 9-12 ):

from sympy import *
x = symbols('x')
print(diff(10*x**4, x, x, x))
print(diff(10*x**4, x, 3))
Listing 9-12.prog12.py

输出如下所示:

240*x
240*x

我们也可以使用integrate()方法计算带有 SymPy 的积分。下面的代码(清单 9-13 )演示了这一点:

from sympy import *
x = symbols('x')
print(integrate(sin(x), x))
Listing 9-13.prog13.py

输出如下所示:

-cos(x)

我们还可以集成如下限制(列表 9-14 ):

from sympy import *
x = symbols('x')
print(integrate(exp(-x), (x, 0, oo)))
Listing 9-14.prog14.py

这里，我们对-x的指数从零到无穷大进行积分(用oo表示)。运行这个并检查输出。我们也可以计算多重极限的多重积分如下(列表 9-15 ):

from sympy import *
x, y = symbols('x y')
print(integrate(exp(-x)*exp(-y), (x, 0, oo), (y, 0, oo)))
Listing 9-15.prog15.py

运行代码(列表 9-15 )并通过手动执行集成来验证输出。

Note

SymPy 真是个大话题。这不可能在一章中完全涵盖。我建议读者在 http://docs.sympy.org 和 www.sympy.org 网站上多探索。

结论

在这一章中，我们从 SymPy 开始，学习了如何在 Python 中执行符号计算。在下一章，我们将从 NumPy 和 Matplotlib 开始。

十、NumPy 简介

在上一章中，我们学习了如何安装 SciPy 栈，以及如何在 Python 3 中使用 SymPy 进行符号计算。在这一章中，我们将被介绍到 NumPy 库，我们将学习 NumPy 的基础知识。我们还将学习使用 Matplotlib 绘制和可视化数据的基础知识。因此，让我们通过学习 NumPy 的基础知识，开始进入科学计算世界的激动人心的旅程。

NumPy 基础

NumPy 是数字(al) Python 的缩写。它的网站( http://www.numpy.org )上说:

努比切望 python(python 语言)

进行科学计算的基础包

其特点如下:

它有一个强大的自定义 N 维数组对象，可以高效方便地表示数据。

它拥有与其他用于科学编程的编程语言(如 C/C++和 FORTRAN)集成的工具。

它用于数学运算，如线性代数、矩阵运算、图像处理和信号处理。

朱皮特

到目前为止，我们一直将代码保存在.py文件中，并用 Python 3 解释器运行它。在这一章中，我们将使用一个叫做 Jupyter 的工具，这是一个基于网络的高级工具，用于在编程语言 Julia、Python 和 R .

Ju lia + pyt hon + r = Jupyter

它将 Python 3(或任何其他支持的语言，如 R 和 Julia)代码和结果保存在一种叫做笔记本的交互格式中。Jupyter 使用 Python 2 和 Python 3 的 IPython 内核。IPython 是 Python 的高级交互式 shell，具有可视化功能。Jupyter 项目是 IPython 的副产品。

Jupyter 和 IPython 具有以下特性:

  • 基于终端和 Qt 的交互式 shells
  • 基于浏览器的笔记本，支持代码和交互式可视化
  • 支持并行计算

可以使用以下命令轻松安装它:

sudo pip3 install --upgrade pip
sudo pip3 install jupyter

这将安装 Jupyter 及其对 Raspberry Pi 的所有依赖。

jupyter 笔记型电脑

Jupyter Notebook 是由 Jupyter Notebook 应用程序生成的文档，其中包含 Python 代码和丰富的文本元素，如段落、等式、图形、链接和交互式可视化。笔记本有人类可读组件和机器可读(可执行)组件。

让我们现在就开始使用 NumPy 和 Jupyter 笔记本吧。打开 lxterminal 并运行以下命令序列:

cd ∼
cd book
cd code
mkdir chapter10
cd chapter10

这将创建并导航到对应于当前章节chapter10的目录。

现在，这是我们的笔记本启动文件夹，让我们从这里使用以下命令启动笔记本:

jupyter notebook

它将启动 Jupyter 笔记本应用程序，并打开一个浏览器窗口(在 Raspbian 的最新版本中为 Chromium 浏览器)。

下面(图 10-1 )是 Jupyter 笔记本启动时的控制台截图:

A447085_1_En_10_Fig1_HTML.jpg

图 10-1。

Jupyter Notebook App console

下面(图 10-2 )是 Chromium 浏览器窗口标签运行笔记本 app 的截图:

A447085_1_En_10_Fig2_HTML.jpg

图 10-2。

Jupyter Notebook App running in Chromium

在浏览器窗口的右上方，单击 New，然后在随后的下拉菜单中选择 Python 3。见下图(图 10-3 )截图:

A447085_1_En_10_Fig3_HTML.jpg

图 10-3。

New Python 3 notebook

它将在同一浏览器窗口中打开一个新的笔记本标签(图 10-4 )。

A447085_1_En_10_Fig4_HTML.jpg

图 10-4。

Python 3 notebook tab

将 Jupyter 笔记本的名称改为Chapter10_Practice，如下图截图所示(图 10-5 )。

A447085_1_En_10_Fig5_HTML.jpg

图 10-5。

Renaming the notebook

笔记本应用程序将显示一个新笔记本的实例，其名称已更新，状态为“正在运行”，如下图截图所示(图 10-6 )。

A447085_1_En_10_Fig6_HTML.jpg

图 10-6。

Notebook running

现在，如果你查看Chapter10目录，你会发现一个与笔记本对应的文件Chapter10_Practice.ipynb。

在窗口顶部的菜单栏中，有您在任何其他 IDE 中都有的选项，例如保存、复制、粘贴和运行。

在第一个单元格中键入import numpy as np,然后单击 Run 按钮。该控件将自动创建下一个文本单元格，并将光标置于其上，如下所示(图 10-7 )。

A447085_1_En_10_Fig7_HTML.jpg

图 10-7。

Working with Python 3 code

我们刚刚将 NumPy 导入到我们的笔记本中，因此我们不必再次导入。此外，我们也可以编辑笔记本的前一个单元格。在执行类型中，如果解释器突出了语法中的错误，我们可以通过编辑任何单元格来修复它。随着科学计算的发展，我们将对 Jupyter 有更多的了解。

N 维数组(ndarray)

NumPy 最强大的构造是 N 维数组(ndarray)。ndarray为多维同质数据提供通用容器。同构意味着ndarray中的数据项是相同的数据类型。让我们看一个 Python 中各种ndarray类型变量的例子。在笔记本中键入以下代码:

x = np.array([1, 2, 3], np.int16)
y = np.array([[0, 1, 2], [3, 4, 5]], np.int32)
z = np.array([[[0, 1, 2], [2, 3, 4], [4, 5, 6]],[[1, 1, 1], [0, 0, 0], [1, 1, 1]]], np.float16)

恭喜，我们刚刚创建了一维、二维和三维的ndarray对象。我们可以通过在 Jupyter notebook 中运行以下代码来验证这一点:

print(x.shape)
print(y.shape)
print(z.shape)

输出是:

(3,)
(2, 3)
(3, 3, 3)

NumPy 中的ndarray的索引方案与 C 语言中的相同，即第一个元素从 0 开始索引。下面一行打印了一个ndarray中第二行第三列的元素值。

我们可以按如下方式对数组进行切片:

print(z[:,1])

以下是输出:

[[ 2\.  3\.  4.]
 [ 0\.  0\.  0.]]

n 数组属性

以下是对ndarray重要属性的演示:

print(z.shape)
print(z.ndim)
print(z.size)
print(z.itemsize)
print(z.nbytes)
print(z.dtype)

以下是输出:

(2, 3, 3)
3
18
2
36
float16

下面是正在发生的事情:

  • ndarray.shape返回数组维数的元组。
  • ndarray.ndim返回数组维数。
  • ndarray.size返回数组中元素的个数。
  • ndarray.itemsize以字节为单位返回一个数组元素的长度。
  • ndarray.nbytes返回数组元素消耗的总字节数。它是属性ndarray.size和ndarray.itemsize相乘的结果。
  • ndarray.dtype返回数组中数据项的数据类型。

数据类型

数据类型用于表示数据的类型和大小。例如，int16表示 16 位有符号整数，uint8表示 8 位无符号整数。NumPy 支持广泛的数据类型。有关数据类型的更多详细信息，请浏览数据类型对象网页( https://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html )。

Note

要查看 NumPy 支持的各种数据类型的详细列表，请访问 https://docs.scipy.org/doc/numpy/user/basics.types.html 。

具备了 NumPy 中的ndarray构造的基础，我们将在下一节开始数组创建例程。

数组创建例程

让我们用内置数组创建例程来创建ndarrays。我们还将使用 matplotlib 来可视化创建的数组。我们将使用 matplotlib 在需要时可视化数据。这本书的最后一章是专门介绍 matplotlib 的，在这里我们会学到更多的内容。

创建ndarrays的第一种方法是ones()。您一定已经猜到了，它创建了所有值都为 1 的数组。在笔记本中运行以下语句，查看该方法的演示:

a = np.ones(5)
print(a)
b = np.ones((3,3))
print(b)
c = np.ones((5, 5), dtype=np.int16)
print(c)

类似地，我们有生成多维零数组的np.zeros()方法。只需通过向方法np.zeros()传递相同的一组参数来运行上面的例子。np.eye()方法用于创建对角矩阵。我们也可以指定对角线的指数。在笔记本中运行以下示例进行演示:

a = np.eye(5, dtype=np.int16)
print(a)
b = np.eye(5, k=1)
print(b)
c = np.eye(5, k=-1)
print(c)

np.arange()返回指定范围内等距数字的列表。尝试笔记本中的以下示例:

a = np.arange(5)
print(a)
b = np.arange(3, 25, 2)
print(b)

输出如下所示:

[0 1 2 3 4]
[ 3  5  7  9 11 13 15 17 19 21 23]

现在让我们用 matplotlib 找点乐子吧。在笔记本中运行以下示例:

import matplotlib.pyplot as plt
x = np.arange(10)
y = np.arange(10)
print(x)
print(y)
plt.plot(x, y, 'o')
plt.show()

在上面的例子中，我们用图形表示了ndarrays。文本结果如下:

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

图形输出(图 10-8 )如下:

A447085_1_En_10_Fig8_HTML.jpg

图 10-8。

Graphical Output of arange()

通过语句import matplotlib.pyplot as plt，我们正在导入myplotlib库的 pyplot 模块。plt.plot()准备显示的标绘图，plt.show()在屏幕上显示。

我们有一个类似的方法，np.linspace(),来生成一个线性数组。linspace()和arange()的主要区别在于linspace()接受要生成的样本数作为参数，而不是步长。下面的代码片段展示了linspace()生成数据的方式:

N = 8
y = np.zeros(N)
x1 = np.linspace(0, 10, N)
x2 = np.linspace(0, 10, N)
plt.plot(x1, y - 0.5, 'o')
plt.plot(x2, y + 0.5, 'o')
plt.ylim([-1, 1])
plt.show()

plt.ylim()指定图形上 Y 坐标的极限。以下(图 10-9 )为输出:

A447085_1_En_10_Fig9_HTML.jpg

图 10-9。

linspace() graph

类似地，我们有logspace()方法，它生成对数标度的数组值。下面的代码演示了:

N = 64
x = np.linspace(1, N, N)
y = np.logspace(0.1, 1, N)
plt.plot(x, y, 'o')
plt.show()

运行笔记本中的语句，它会生成以下(图 10-10 )输出:

A447085_1_En_10_Fig10_HTML.jpg

图 10-10。

linspace() versus logspace() graph

矩阵和线性代数

NumPy 有创建矩阵的例程。我们来看几个例子。np.matrix()将给定数据解释为矩阵。以下是一些例子:

a = np.matrix('1 2; 3 4')
b = np.matrix([[1, 2], [3, 4]])
print(a)
print(b)

运行笔记本中的代码，您会看到两个示例都返回了矩阵，即使输入的格式不同。

[[1 2]
 [3 4]]
[[1 2]
 [3 4]]

从数组或序列构建一个矩阵。

A = np.mat('1 1; 1 1')
B = np.mat('2 2; 2 2')
C = np.mat('3 4; 5 6')
D = np.mat('7 8; 9 0')
a = np.bmat([[A, B], [C, D]])
print(a)

上面的代码返回一个由所有序列组合而成的矩阵。运行笔记本中的代码，查看如下输出:

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

np.matlib.zeros()和np.matlib.ones()分别返回 0 和 1 的矩阵。np.matlib.eye()返回一个对角矩阵。np.matlib.identity()返回给定大小的正方形单位矩阵。下面的代码示例演示了这些方法:

from numpy.matlib import *
a = zeros((3, 3))
print(a)
b = ones((3, 3))
print(b)
c = eye(3)
print(c)
d = eye(5, k=1)
print(d)
e = eye(5, k=-1)
print(e)
f = identity(4)
print(f)

rand()和randn()方法返回带有随机数的矩阵。

a = rand((3, 3))
b = randn((4, 4))
print(a)
print(b)

我们来学习几个与线性代数(矩阵运算)相关的方法。我们有dot()方法，它计算两个数组的点积，而vdot()计算两个向量的点积。inner()计算两个数组的内积。outer()计算两个向量的外积。下面的代码示例演示了所有这些方法:

a = [[1, 0], [0, 1]]
b = [[4, 1], [2, 2]]
print(np.dot(a, b))
print(np.vdot(a, b))
print(np.inner(a, b))
print(np.outer(a, b))

输出如下所示:

[[4 1]
 [2 2]]
6
[[4 2]
 [1 2]]
[[4 1 2 2]
 [0 0 0 0]
 [0 0 0 0]
 [4 1 2 2]]

三角测量法

让我们把三角方法形象化。我们将使用 matplotlib 可视化sin()、cos()、tan()、sinh()和cosh()方法。下面的示例演示了这些方法的用法:

x = np.linspace(-np.pi, np.pi, 201)
plt.plot(x, np.sin(x))
plt.xlabel('Angle in radians')
plt.ylabel('sin(x)')
plt.show()

输出(图 10-11 )如下:

A447085_1_En_10_Fig11_HTML.jpg

图 10-11。

Graph for sin(x)

以下代码示例演示了cos()的用法:

x = np.linspace(-np.pi, 2*np.pi, 401)
plt.plot(x, np.cos(x))
plt.xlabel('Angle in radians')
plt.ylabel('cos(x)')
plt.show()

以下(图 10-12 )为输出:

A447085_1_En_10_Fig12_HTML.jpg

图 10-12。

graph for cos(x)

让我们继续讨论双曲余弦和正弦波，如下所示:

x = np.linspace(-5, 5, 2000)
plt.plot(x, np.cosh(x))
plt.show()
plt.plot(x, np.sinh(x))
plt.show()

以下(图 10-13 )是cosh()的输出:

A447085_1_En_10_Fig13_HTML.jpg

图 10-13。

Graph of cosh(x)

下面(图 10-14 )是sinh(x)的曲线图:

A447085_1_En_10_Fig14_HTML.jpg

图 10-14。

Graph of sinh(x)

随机数和统计

rand()方法生成给定维数的随机矩阵。randn()方法使用从正态分布中抽取的数字生成一个矩阵。randint()在给定的界限内生成一个数，不包括这个界限。random_integers()生成包含给定限值的随机整数。下面的代码演示了上述方法的前三种:

import numpy as np
a = np.random.rand(3, 3)
b = np.random.randn(3, 3)
c = np.random.randint(4, 15)
print(a)
print(b)
print(c)

我们可以计算 a ndarray的中值、平均值、均值、标准差和方差，如下所示:

a = np.array([[10, 7, 4], [1, 2, 3], [1, 1, 1]])
print(median(a))
print(average(a))
print(mean(a))
print(std(a))
print(var(a))

傅立叶变换

NumPy 有一个用于基本信号处理的模块。fft模块具有fft()方法，用于计算一维离散傅立叶变换，如下所示:

t = np.arange(256)
sp = np.fft.fft(np.sin(t))
freq = np.fft.fftfreq(t.shape[-1])
plt.plot(freq, sp.real, freq, sp.imag)
plt.show()

输出(图 10-15 )如下:

A447085_1_En_10_Fig15_HTML.jpg

图 10-15。

Fast Fourier Transform

fft2()和fftn()分别用于计算二维和 n 维离散傅立叶变换。试着为这些写代码。

结论

在本章中，我们从 NumPy、matplotlib 和 Jupyter 开始，并学习了如何在 Python 中执行数值计算和基本的可视化。在下一章，我们将从 SciPy 库开始。

十一、SciPy 简介

在上一章中，我们学习了如何使用 NumPy 进行数值计算。我们还学习了如何方便地使用 Jupyter，以及如何使用 matplotlib 进行可视化。在本章中，我们将介绍 SciPy 库。然而，我们与 NumPy 和 matplotlib 的旅程还远未结束。在本书的剩余部分，我们将学习 NumPy 和 matplotlib 的新功能。因此，让我们和 SciPy 一起踏上科学计算的激动人心的旅程。

SciPy 中的科学和数学常数

在我们开始之前，让我们完成为本章创建一个新目录的仪式。请为本章创建一个目录chapter 11。

cd ∼
cd book
cd code
mkdir chapter11
cd chapter11

现在，使用以下命令启动 Jupyter 笔记本应用程序:

jupyter notebook

打开一个新笔记本，重命名为Chapter11_Practice。这个笔记本将保存本章的代码。

SciPy 库有一个名为scipy.constants的模块，其中有许多重要的数学和科学常数的值。让我们试几个。在笔记本中运行以下代码:

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import *
print("Pi = " + str(pi))
print("The golden ratio = " + str(golden))
print("The speed of light = " + str(c))
print("The Planck constant = " + str(h))
print("The standard acceleration of gravity = " + str(g))
print("The Universal constant of gravity = " + str(G))

输出如下所示:

Pi = 3.141592653589793
The golden ratio = 1.618033988749895
The speed of light = 299792458.0
The Planck constant = 6.62606957e-34,
The standard acceleration of gravity = 9.80665
The Universal constant of gravity = 6.67384e-11

请注意，SciPy 常量不包括度量单位，仅包括常量的数值。这些在数值计算中非常有用。

Note

还有更多这样的常数。请访问 https://docs.scipy.org/doc/scipy/reference/constants.html 查看更多。

线性代数

现在我们将学习几个与线性代数相关的方法。让我们从矩阵的逆矩阵开始:

import numpy as np
from scipy import linalg
a = np.array([[1, 4], [9, 16]])
b = linalg.inv(a)
print(b)

以下是输出:

[[-0.8   0.2 ]
 [ 0.45 -0.05]]

我们也可以求解矩阵方程ax = b如下:

a = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]])
b = np.array([2, 4, -1])
from scipy import linalg
x = linalg.solve(a, b)
print(x)
print(np.dot(a, x))

以下是输出:

[ 2\. -2\.  9.]
[ 2\.  4\. -1.]

我们可以如下计算矩阵的行列式:

a = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
print(linalg.det(a))

在笔记本中运行上面的代码，自己看看结果。

我们还可以计算矩阵的范数，如下所示:

a = np.arange(16).reshape((4, 4))
print(linalg.norm(a))
print(linalg.norm(a, np.inf))
print(linalg.norm(a, 1))
print(linalg.norm(a, -1))

以下是结果:

35.2136337233
54
36
24

我们也可以计算 QR 和 RQ 分解如下:

from numpy import random
from scipy import linalg
a = random.randn(3, 3)
q, r = linalg.qr(a)
print(a)
print(q)
print(r)
r, q = linalg.rq(a)
print(r)
print(q)

综合

SciPy 有用于各种集成操作的integrate模块，所以让我们看看它的一些方法。第一个是quad()。它接受要积分的函数以及积分的极限作为参数，然后返回数值和近似误差。以下是几个例子:

from scipy import integrate
f1 = lambda x: x**4
print(integrate.quad(f1, 0, 3))
import numpy as np
f2 = lambda x: np.exp(-x)
print(integrate.quad(f2, 0, np.inf))

以下是输出结果:

(48.599999999999994, 5.39568389967826e-13)
(1.0000000000000002, 5.842606742906004e-11)

trapz()使用梯形法则沿给定轴积分:

print(integrate.trapz([1, 2, 3, 4, 5]))

以下是输出:

12.0

让我们看一个使用梯形法则的累积积分的例子。

import matplotlib.pyplot as plt
x = np.linspace(-2, 2, num=30)
y = x
y_int = integrate.cumtrapz(y, x, initial=0)
plt.plot(x, y_int, 'ro', x, y[0] + 0.5 * x**2, 'b-')
plt.show()

以下(图 11-1 )为输出:

A447085_1_En_11_Fig1_HTML.jpg

图 11-1。

Cumulative integration using the trapezoidal rule

插入文字

这个模块有插值的方法。让我们借助 matplotlib 的图形表示来研究其中的一些。interp1d()用于一维插值，如下所示。

from scipy import interpolate
x = np.arange(0, 15)
y = np.exp(x/3.0)
f = interpolate.interp1d(x, y)
xnew = np.arange(0, 14, 0.1)
ynew = f(xnew)
plt.plot(x, y, 'o', xnew, ynew, '-')
plt.show()

结果如下(图 11-2 ):

A447085_1_En_11_Fig2_HTML.jpg

图 11-2。

1-D interpolation

interp1d()适用于y=f(x)类型的功能。类似地，interp2d()处理z=f(x, y)类型的函数。它用于二维插值。

x = np.arange(-5.01, 5.01, 0.25)
y = np.arange(-5.01, 5.01, 0.25)
xx, yy = np.meshgrid(x, y)
z = np.sin(xx**3 + yy**3)
f = interpolate.interp2d(x, y, z, kind='cubic')
xnew = np.arange(-5.01, 5.01, 1e-2)
ynew = np.arange(-5.01, 5.01, 1e-2)
znew = f(xnew, ynew)
plt.plot(x, z[0, :], 'ro-', xnew, znew[0, :], 'b-')
plt.show()

结果如下(图 11-3 ):

A447085_1_En_11_Fig3_HTML.jpg

图 11-3。

2-D interpolation

接下来，我们要学习splev()、splprep()和splrep()。splev()方法用于计算 B 样条或其导数。我们将把这种方法与用于表示 B 样条的splprep()和splrep(),一起使用。

splrep()用于表示一维曲线，如下所示:

from scipy.interpolate import splev, splrep
x = np.linspace(-10, 10, 50)
y = np.sinh(x)
spl = splrep(x, y)
x2 = np.linspace(-10, 10, 50)
y2 = splev(x2, spl)
plt.plot(x, y, 'o', x2, y2)
plt.show()

结果如下(图 11-4 ):

A447085_1_En_11_Fig4_HTML.jpg

图 11-4。

Representation of a 1-D curve

splprep()用于表示一条 N 维曲线。

from scipy.interpolate import splprep
theta = np.linspace(0, 2*np.pi, 80)
r = 0.5 + np.cosh(theta)
x = r * np.cos(theta)
y = r * np.sin(theta)
tck, u = splprep([x, y], s=0)
new_points = splev(u, tck)
plt.plot(x, y, 'ro')
plt.plot(new_points[0], new_points[1], 'r-')
plt.show()

下面(图 11-5 )是结果。

A447085_1_En_11_Fig5_HTML.jpg

图 11-5。

Representation of an N-D curve

结论

在本章中，我们介绍了 SciPy 库中几个重要且常用的模块。接下来的两章将着重向读者介绍信号和图像处理的专业科学领域。在下一章，我们将学习几个与信号处理相关的模块。

十二、SciPy 信号处理

在上一章中，我们学习了如何用 SciPy 进行科学计算。我们了解了 SciPy 库的几个模块。在本章中，我们将探索一个重要的科学领域，即信号处理，并且我们将学习SciPy.signal模块中的方法。让我们从 SciPy 中的信号处理开始吧。这将是一个简短的章节，只有几个代码示例，让您对信号处理领域的一些基础知识有所了解。

波形

让我们从波形发生器功能开始。在∼/book/code目录下新建一个目录chapter 12。运行以下命令启动 Jupyter 笔记本应用程序:

jupyter notebook

将当前笔记本重命名为Chapter12_Practice。和前面的章节一样，在同一个笔记本中运行本章的所有代码。先导入 NumPy 和 matplotlib。

import numpy as np
import matplotlib.pyplot as plt

我们要研究的第一个例子是锯齿波发生器功能。

from scipy import signal
t = np.linspace(0, 2, 5000)
plt.plot(t, signal.sawtooth(2 * np.pi * 4 * t))
plt.show()

该函数接受信号的时序和宽度，并生成三角形或锯齿形连续信号。以下(图 12-1 )为输出，

A447085_1_En_12_Fig1_HTML.jpg

图 12-1。

Sawtooth wave signal

让我们来看一个方波发生器，它接受时间序列和占空比作为输入。

t = np.linspace(0, 1, 400)
plt.plot(t, signal.square(2 * np.pi * 4 * t))
plt.ylim(-2, 2)
plt.title('Square Wave')
plt.show()

输出如下所示:

A447085_1_En_12_Fig2_HTML.jpg

图 12-2。

Square wave signal

脉宽调制方波正弦波可以如下所示:

sig = np.sin(2 * np.pi * t)
pwm = signal.square(2 * np.pi * 30 * t, duty=(sig +1)/2)
plt.subplot(2, 1, 1)
plt.plot(t, sig)
plt.title('Sine Wave')
plt.subplot(2, 1, 2)
plt.plot(t, pwm)
plt.title('PWM')
plt.ylim(-1.5, 1.5)
plt.show()

输出(图 12-3 )如下:

A447085_1_En_12_Fig3_HTML.jpg

图 12-3。

Modulated wave

窗口功能

窗口函数是一种数学函数，在特定区间外为零。我们现在来看看三种不同的窗口函数。第一个是汉明窗函数。我们必须将输出窗口中的点数作为参数传递给所有函数。

window = signal.hamming(101)
plt.plot(window)
plt.title('Hamming Window Function')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.show()

输出(图 12-4 )如下:

A447085_1_En_12_Fig4_HTML.jpg

图 12-4。

Hamming window demo

汉宁窗函数如下:

window = signal.hanning(101)
plt.plot(window)
plt.title('Hanning Window Function')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.show()

输出(图 12-5 )如下:

A447085_1_En_12_Fig5_HTML.jpg

图 12-5。

Hanning window demo

凯泽窗函数如下:

window = signal.kaiser(51, beta=20)
plt.plot(window)
plt.title('Kaiser Window Function Beta = 20')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.show()

输出(图 12-6 )如下:

A447085_1_En_12_Fig6_HTML.jpg

图 12-6。

Kaiser window demo

墨西哥帽小波

通过将点数和振幅作为参数传递，我们可以使用 Ricker 函数生成一个墨西哥帽小波，如下所示:

plt.plot(signal.ricker(200, 6.0))
plt.show()

墨西哥帽小波是连续小波族中的一个特例。它用于过滤和平均光谱信号。输出如下所示:

A447085_1_En_12_Fig7_HTML.jpg

图 12-7。

Mexican hat wavelet

盘旋

我们可以用convolve()方法卷积两个 N 维数组，如下所示:

sig = np.repeat([0., 1., 0.], 100)
win = signal.hann(50)
filtered = signal.convolve(sig, win, mode='same') / sum(win)
plt.subplot(3, 1, 1)
plt.plot(sig)
plt.ylim(-0.2, 1.2)
plt.title('Original Pulse')
plt.subplot(3, 1, 2)
plt.plot(win)
plt.ylim(-0.2, 1.2)
plt.title('Filter Impulse Response')
plt.subplot(3, 1, 3)
plt.plot(filtered)
plt.ylim(-0.2, 1.2)
plt.title('Filtered Signal')
plt.show()

信号、窗口及其卷积如图 12-8 所示。两个信号的卷积提供了两个信号的组合以形成第三个信号。它是信号处理中一种非常重要的信号合并技术。如果信号代表图像或音频数据，我们可以基于卷积度获得改进的图像或音频。如果您还没有注意到，我们正在使用 NumPy 的repeat()方法，它把要重复的模式和重复次数作为参数。在本例中，我们生成一个样本大小为 300 的信号。

A447085_1_En_12_Fig8_HTML.jpg

图 12-8。

Convolution

结论

在这简短的一章中，我们介绍了 SciPy 的scipy.signal模块中的几个重要的方法类。在下一章，我们将探索图像处理领域。

十三、SciPy 图像处理

在上一章中，我们用 SciPy 研究了信号处理。我们研究了 SciPy 提供的几个主要函数类，合并在scipy.signal下。在这一章中，我们将学习一个 SciPy 包，scipy.misc，以及一些使用它进行图像处理的例子。

第一图像处理程序

scipy.misc模块用于基本的图像处理操作。创建一个名为Dataset的目录来存储我们将要使用的所有样本图像。

cd ∼/book
mkdir Dataset

样本图像数据集可以在本书的页面上的在线下载部分获得。同样，为图书代码创建一个目录chapter 13，如下所示:

cd ∼/book/code
mkdir chapter13
cd chapter13

让我们看一个读取和显示图像的基本例子(清单 13-1 )。

from scipy import misc

img = misc.imread('/home/pi/book/Dataset/4.2.01.tiff')

misc.imshow(img)

Listing 13-1.prog01.py

代码(清单 13-1 将从imread()方法提供的路径中读取一个图像，并且imshow()方法将使用xlib显示它。

scipy.misc有三个内置图像。这些可以如下使用(列表 13-2 ):

from scipy import misc

img1 = misc.face()
img2 = misc.lena()
img3 = misc.ascent()

misc.imshow(img1)
misc.imshow(img2)
misc.imshow(img3)

Listing 13-2.prog02.py

face()是浣熊的脸，lena()是标准测试图像，ascent()是灰度图像。

简单的图像处理

scipy.misc有三种简单操作的方法。scipy.imfilter()对图像应用各种滤镜。以下(清单 13-3 )是一个例子:

from scipy import misc

misc.imshow(misc.imfilter(misc.face(), 'edge_enhance_more'))

Listing 13-3.prog03.py

在上面的代码中(清单 13-3 ，我没有使用任何中间变量来存储图像。我通过传递给imshow()方法来直接显示它。方法imfilter()接受两个参数。

  • 第一个是要过滤的图像。
  • 第二个是要应用的预定义过滤器的类型。

过滤器类型的允许值为'blur', 'contour', 'detail', 'edge_enhance', 'edge_enhance_more', 'emboss', 'find_edges', 'smooth', 'smooth_more', 'sharpen'。

我们可以将图像的大小调整为 50%，如下所示(清单 13-4 ):

from scipy import misc
misc.imshow(misc.imresize(misc.face(), 50))
Listing 13-4.prog04.py

我们也可以将图像旋转一定角度，如下所示(列表 13-5 ):

from scipy import misc
misc.imshow(misc.imrotate(misc.face(), 45))
Listing 13-5.prog05.py

用于图像处理的 NumPy 简介

让我们从使用 NumPy 库进行图像处理的基础开始。考虑以下(列表 13-6 )代码:

from scipy import misc

img = misc.face()

print(type(img))

Listing 13-6.prog06.py

上述程序(清单 13-6 )的输出如下:

<class 'numpy.ndarray'>

这意味着图像的数据类型是 NumPy 中的ndarray。我们需要理解一些重要的ndarray属性，这将帮助我们理解它所代表的图像的重要属性。

考虑下面的代码(清单 13-7 ):

from scipy import misc

img = misc.face()

print(img.dtype)
print(img.shape)
print(img.ndim)
print(img.size)

Listing 13-7.prog07.py

输出如下所示:

uint8
(768, 1024, 3)
3
2359296

让我们来了解一下其中每一项的含义。

  • dtype属性用于表示图像的元素的数据类型。在这种情况下，它是uint8，表示无符号的 8 位整数。这意味着它可以有 256 个不同的值。
  • shape表示图像的尺寸。在这种情况下，它是一个彩色图像。它的分辨率为 1024x768，有三个颜色通道，分别对应于红色、绿色和蓝色。每个像素的每个通道可以有 256 个可能值中的一个。因此，这些的组合可以为每个像素产生 256256256 种不同的颜色。您可以将彩色图像想象为三个二维平面的排列。灰度图像是灰度值的单一平面。
  • ndim代表尺寸。彩色图像具有三维，而灰度图像具有二维。
  • size代表数组中元素的总数。它可以通过乘以尺寸值来计算。在本例中，它是 76810243=2359296。

我们可以看到对应于单个像素的 RGB 值如下(列表 13-8 ):

from scipy import misc

img = misc.face()

print(img[10, 10]))

Listing 13-8.prog08.py

在上面的代码中(清单 13-8 ，我们正在访问位于(10，10)的像素值。输出为[172 169 188]。

这些是 NumPy 用于图像处理的基础。在本章中，我们将在需要时学习更多关于 NumPy 的知识。

用于图像处理的 Matplotlib

我们使用了misc.imshow()方法来显示图像。虽然该方法对于简单的应用程序很有用，但它是原始的。我们需要使用更先进的科学应用框架。我们知道 matplotlib 就是为此服务的。这是一个用于 Python 的 MATLAB 风格的绘图和数据可视化库，我们在安装 SciPy 栈时安装了它。我们在前面的章节中也使用过它。在这一章和下一章，我们将使用它来显示图像。它是 SciPy 栈不可或缺的一部分。就像 NumPy 一样，matplotlib 对于一本书来说是一个太大的主题，值得再写一本书。我们将只使用 matplotlib 中的pyplot模块来满足我们的图像处理需求。让我们看一个简单的图像处理程序(列表 13-9 )如下:

import scipy.misc as misc
import matplotlib.pyplot as plt

img = misc.face()

plt.imshow(img)
plt.show()

Listing 13-9.prog09.py

在上面的代码中(清单 13-9 ，我们正在导入pyplot模块。imshow()方法用于将图像添加到绘图窗口。show()方法显示绘图窗口。输出(图 13-1 )如下:

A447085_1_En_13_Fig1_HTML.jpg

图 13-1。

pyplot imshow() demo for color image

我们还可以关闭轴(或标尺)并给图像添加一个标题，如下所示(清单 13-10 ):

import scipy.misc as misc
import matplotlib.pyplot as plt

img = misc.lena()

plt.imshow(img, cmap='gray')
plt.axis('off')
plt.title('face')
plt.show()

Listing 13-10.prog10.py

由于图像是灰度图像，我们必须在imshow()方法中选择一个gray色彩映射表，以便图像的色彩空间正确显示在绘图窗口中。axis('off')用于关闭斧子。title()方法用于指定图像的标题。输出(图 13-2 )如下:

A447085_1_En_13_Fig2_HTML.jpg

图 13-2。

Lena image with title and axis off

我们可以使用imshow()将多个图像推送到绘图窗口中的图像网格(参见清单 13-11 )，如下所示:

import scipy.misc as misc
import matplotlib.pyplot as plt

img1 = misc.face()
img2 = misc.ascent()
img3 = misc.lena()

titles = ['face', 'ascent', 'lena']
images = [img1, img2, img3]

plt.subplot(1, 3, 1)
plt.imshow(images[0])
plt.axis('off')
plt.title(titles[0])

plt.subplot(1, 3, 2)
plt.imshow(images[1], cmap='gray')
plt.axis('off')
plt.title(titles[1])

plt.subplot(1, 3, 3)
plt.imshow(images[2], cmap='gray')
plt.axis('off')
plt.title(titles[2])

plt.show()

Listing 13-11.prog11.py

在imshow()之前我们已经用过subplot()的方法了。subplot()方法的前两个参数指定网格的尺寸，第三个参数指定图像在网格中的位置。网格中图像的位置编号从左上角开始。左上位置是第一个位置，下一个位置是第二个位置，依此类推。让我们看看输出(图 13-3 ):

A447085_1_En_13_Fig3_HTML.jpg

图 13-3。

Multiple image grid

图像通道

我们可以分离多通道图像的图像通道。代码(列表 13-12 )如下:

import scipy.misc as misc
import matplotlib.pyplot as plt

img = misc.face()

r = img[:, :, 0]
g = img[:, :, 1]
b = img[:, :, 2]

titles = ['face', 'Red', 'Green', 'Blue']
images = [img, r, g, b]

plt.subplot(2, 2, 1)
plt.imshow(images[0])
plt.axis('off')
plt.title(titles[0])

plt.subplot(2, 2, 2)
plt.imshow(images[1], cmap='gray')
plt.axis('off')
plt.title(titles[1])

plt.subplot(2, 2, 3)
plt.imshow(images[2], cmap='gray')
plt.axis('off')
plt.title(titles[2])

plt.subplot(2, 2, 4)
plt.imshow(images[3], cmap='gray')
plt.axis('off')
plt.title(titles[3])

plt.show()

Listing 13-12.Prog12.py

代码(清单 13-12 )的输出(图 13-4 )如下:

A447085_1_En_13_Fig4_HTML.jpg

图 13-4。

Separate image channels

我们可以使用np.dstack()方法合并所有通道来创建原始图像，如下所示(列表 13-13 ):

import scipy.misc as misc
import matplotlib.pyplot as plt
import numpy as np

img = misc.face()

r = img[:, :, 0]
g = img[:, :, 1]
b = img[:, :, 2]

output = np.dstack((r, g, b))

plt.imshow(output)
plt.axis('off')
plt.title('Combined')
plt.show()

Listing 13-13.prog13.py

运行上面的代码(清单 13-13 ，亲自看看np.dstack()是如何工作的。

Note

我写过一本关于用 Raspberry Pi 进行图像处理的详细书籍。可以在www . a press . com/us/book/9781484227305购买。

结论

在本章中，我们了解了使用 SciPy 进行图像处理的世界。我们还学习了如何使用 NumPy ndarray表示图像。我们学习了用scipy.misc包对图像执行一些基本操作。下一章，我们将学习 matplotlib 的数据表示和处理。

avatar
文章
阅读
粉丝