LWN: 介绍一下SciPy！

关注了就能看到更多这么棒的文章哦～

An introduction to SciPy

January 19, 2021

This article was contributed by Lee Phillips
DeepL assisted translation
https://lwn.net/Articles/842964/

SciPy 是一组用于科学计算和数值计算的 Python 库。几乎所有真正在使用 Python 进行科学研究的用户都在使用 SciPy。由于 Python 在所有科学领域都很流行，并且在一些研究领域比如如数据科学领域，Python 仍然是主流使用的语言，因此 SciPy 拥有庞大的用户群。在新年前夕，SciPy 公布了 scipy 库的 1.6 版本，这是 SciPy 各个组件中的核心。这个版本给了我们一个很好的机会来深入了解这个软件，并给大家展示一些使用实例。

What is SciPy

SciPy 这个名字其实指的是几个相关的含义。SciPy 既用在一些与科学研究中使用 Python 相关的国际会议的标题中，也是 scipy 库，这个 scipy 库包含了用于科学和数值计算各个领域的多个 python module。比如其中包含：

• scipy.integration，包含用于数值积分和求解微分方程的功能。
  

• scipy.fft, 用于傅立叶变换。
  

• scipy.linalg, 用于线性代数。
  

• scipy.stats，用于统计
  

• scipy.special, 提供了许多数学物理学（mathematical physics）中使用的特殊函数，如 airy, elliptic, bessel, 和 hypergeometric 函数。
  

除了 scipy 之外，SciPy 中还包括了一些有用的组件，这些组件被打包在一起用于科学研究。NumPy 是其中的基础库，用来在 Python 中添加了一个数组数据类型，并封装了 C 和 Fortran 的底层实现，从而使得 Python 的数值计算能力变得真正实用了。还有用于数据科学（data science）的 pandas，和 Python 绘图库 Matplotlib。用户界面 IPython 和 Jupyter 也是集合的一部分，还有 Cython，用于构建 C 扩展。还有其他的一些。所以在指向这个领域的时候，SciPy 代表的是一整套软件生态系统：许多可以一起使用的 module，共同构建被研究的问题的解决方案，并与他人分享。

SciPy 背后的理念是希望让科学家们只需一个安装命令就能拥有一套完整的工具。尽管有一些小毛病，但是主要目的还是达到了。确实不时出现一些卸载不干净的旧版本导致的版本冲突，不过这个项目还是在很大程度上将研究者从寻找十几个不同库的兼容版本并分别安装的艰巨任务中解放了出来。此外，SciPy 还提供了一个社区，以及一个寻求文档和帮助的中心场所。

Using SciPy

SciPy 的应用范围非常广泛，因此，这里举的少量例子，肯定无法涵盖它的绝大多数能力。然而，这几个例子就可以让人认识到大家是如何利用这个项目的。这个主要示例与 LWN 最近关于 Pluto 和 Julia 的文章中的问题类似，这样读者可以比较不同项目提供的解决方法。

当然，第一步是完成安装。该项目提供了官方的安装说明，介绍了一些将库安装好的方法，它使用的是比大多数发行版的包管理器中提供的最新版本还要新的版本。我选择在一个全新的虚拟环境（virtual environment）中使用 pip 安装。安装好之后，我的 IPython 无法按 tab 键完成补全，并且经常崩溃。在 Google 上进行了一些搜索后，发现了一个版本冲突的问题，而大家使用虚拟环境（virtual environment）应该有助于预防这种情况。按照文章建议用命令对 completion module 进行降级就解决了这个问题：

$ python -m pip install jedi==0.17.2 

有使用 Python，有三种主要方式，因此 SciPy 的主要使用方式也同样是：使用常用的编辑器创建脚本文件；使用 Jupyter 笔记本；或者使用 IPython 的 read-eval-print loop（REPL）方式。当然，在 shell 中输入 python 就可以使用 "原有的" REPL 功能，但毕竟 IPython 在此基础上提供了显著增强，所以没有什么理由不使用 IPython。

我们的第一个例子将是我们在讲 Julia 的文章中攻克的普通微分方程（ODE）数值求解的问题。这个问题被称为逻辑方程（logistic equation），具体方程是：

f' = f - f²

主要是时间导数。

下面的列表显示了设置和解决这个问题的 IPython 会话的步骤，中间还穿插了代码的解释。

from scipy.integrate import solve_ivp

首先我们 import solve_ivp，这是 SciPy 的基本 ODE 数值求解器。scipy.integration 是一个 module，它包括各种用于逼近积分和寻找 ODE 解的数值逼近的例程。

# Here we define a function that returns the RHS of the ODE:    
def f(t, f):
    return [f - f**2]

描述方程右侧的 f() 函数需要返回一个列表。在这种情况下，它只包含单个元素，因为方程是一阶的，这意味着它只有一个一阶导数。我们必须将高阶方程重写为一阶方程，在这种情况下，返回的数组将为系统中的每个方程提供一个元素：

# We'll store the time interval for the solution here:
tint = [0, 8]

import numpy as np

# linspace creates an array of equally-spaced values:
t = np.linspace(0, 8, 100)

变量 t 将包含一个从 0 到 8 的数组，以 100 个等距的步长计算。这些将是计算求解的时间。如果我只想知道一个时间间隔结束时的解的值，我就不需要建立这样的数组，但是我打算把求解结果绘制出来，所以我想知道它在很多点的值。

# Now we can solve, storing the solution in sol:
sol = solve_ivp(f, tint, [0.05], t_eval=t)

调用求解器(solver)的第三个参数是第一次取值时的初始条件，在本例中是 0.05。它和函数本身一样，是在一个数组中，这是用来处理多方程系统的情况。返回值存储在 sol 中，是一个对象，有少数几个字段包含了关于解的信息。解本身存储在 sol.y 中，对应的独立变量值存储在 sol.t 中。在这种情况下，sol.t 应该与 t 完全相同，但如果求解器没有提供输入时间数组（input time array），它将决定在哪些时间提供解值，这就需要查询 sol.t：

import matplotlib.pyplot as plt

# A no-frills plot, using the array stored in fields of sol:
plt.plot(sol.t, sol.y[0])

# Here is the command to display the plot:
plt.show()

结果如下图所示。如果我们有一组方程，而不仅仅是一个方程，就我们可以通过向绘图函数传递 sol.y 数组的额外元素来绘制感兴趣的变量。

[Logistic fn plot]

SciPy 的 scipy 组件包含了一些有趣的模块，这些模块集中在不同的领域。例如，让我们来探索一下图像处理库 scipy.ndimage 的部分功能。下面的命令来自另一个 IPython 会话，其中我读取了一张图像并对其进行了四次变换。

from scipy import ndimage
import matplotlib.pyplot as plt

# Read an image file from disk
bogy = plt.imread('bogart.jpg')

ndimage 库中的变换是基于矩阵操作的。imread()函数会读取图像，并将其作为一个 1200×1200×3 的数字矩阵存储在变量 bogy 中，这个矩阵是按原始图像尺寸，针对每个原色（primary color）来生成的。下面使用的每一个图像变换函数都使用这个输入图像作为第一个参数，第二参数就是变换选项。调用 gaussian_filter()和 gaussian_laplace()时使用的 sigma 参数决定了变换的强度。gaussian_laplace()滤波器通过使用二次导数来检测强度的变化，从而对图像进行平滑。

下图是原始图像。每次进行变换（transformation）后，显示的是结果图像。

[Original Bogart]

# The first transformation: Gaussian blurring:
bogyBlurred = ndimage.gaussian_filter(bogy, sigma=17)

[Blurred]

# The rotation function takes degrees:
bogyTilted = ndimage.rotate(bogy, 45)

[tilted]

# This transformation can create sketch-like images
# from photographs:
bogySketched = ndimage.gaussian_laplace(bogy, sigma=3)

[Sketched]

# The geometric_transform routine uses a mapping function
# to transform coordinates:
def mapping(t):
    return (t[0], t[1]**1.3, t[2])

bogySqueezed = ndimage.geometric_transform(bogy, mapping)

[Squeezed]

# Each processed image is stored in its own variable. We could
# apply additional transformations to them or save them
# to disk with the following:
plt.imshow(bogyTilted)
plt.savefig('bogyTilted.png')

ndimage 库还有更多的功能，可不仅仅是用来变换黑白电影明星图像。它可以有效地应用于诸如自动计算显微镜载玻片上的细菌数量，或者测量照片中物体的尺寸等问题。

Changes in 1.6

在 1.6 版本的发布公告中，它被描述为 "6 个月努力工作的结晶。它包含了许多新的功能和大量的 bug 修复，更高的测试覆盖率，以及更好的文档"。在这样一个大型项目中，通常情况下，长长的改动列表中的每一项都只有那些碰巧正在使用这一部分的人感兴趣。但值得注意的是，许多改动是针对上一节探讨的 ndimage 库的，包括对这个库的文档进行了全面的修改。

Help and documentation

SciPy 的广泛采用和悠久的历史，意味着有大量的帮助和参考资料可用。官方文档的起始页是一个索引，指向 SciPy 的每个主要组件的文档。每一个组件都有自己的项目，有自己的团队，这个索引会引导到各个项目网站。例如，其中一个组件是 SymPy，我最近在 LWN 这里介绍过，这个该索引就直接通向 SymPy 网站上维护的文档。

SciPy 团队维护着 Planet SciPy，这是网络上关于 SciPy 的文章和帖子的聚合。其中很多是关于各种主题的教程，无论是初学者还是专家都能找到感兴趣的部分。SciPy Cookbook 也不应该被忽视。这是一个有很多细节、相当广泛的代码样例的集合，给出了使用 SciPy 库来解决各种问题的例子。这些例子从初学者上手到非常专业的高级应用都有。

除了官方来源之外，还有 SciPy Lecture Notes。这是一个全面的文档，它从 Python 语言和 SciPy 生态系统的概述开始，一直深入到许多领域的高级主题。它有许多作者贡献，并通过 GitHub 来寻求社区的贡献。这些图文并茂的讲义（Lecture），除了作为自学的阅读材料外，还被设计成可供教师使用的投影幻灯片。为此，代码被显示在了带有 widget 的 box 中，这些 widget 将每个代码样例折叠成更简洁的形式，去除了提示和一些输出。这使得幻灯片不那么紧凑，对于帮助讲解更加有效率，同时也让在线读者在喜欢的时候可以随时了解所有细节。

物理学家 Robert Johansson 在 GitHub 上维护了一系列关于用 Python 进行科学计算的讲座，可以作为使用 SciPy 的非常有用的信息来源。讲座以 Jupyter 笔记本的形式供下载，但作者也提供了 PDF 版本。

对于那些喜欢读纸质书的人来说，Johansson 的《Numerical Python》一书涵盖了很多这方面的内容。《Python Data Science Handbook》受到数据科学家的推崇，作者 Jake VanderPlas 以传统形式和一套 Jupyter 笔记本的形式免费提供全文阅读。另一本物理学家的书是《Learning Scientific Programming with Python》。它的作者 Christian Hill 维护了一个 GitHub 页面，里面有一些有趣的案例研究，都是以 Jupyter 的形式，与他书中提到的主题相关。

当然也有丰富的参考文档和教程文档，各种形式都有。而当一个人被卡住的时候，SciPy 被人们广泛采用的一个好处就是很容易在 Stack Overflow 这样的网站上找到解决方法。

SciPy or not

Python、NumPy 以及其他那些最终组合成 SciPy 的组件，共同为科学领域的免费软件应用铺平了道路。尽管其他一些语言和库现在可能更胜一筹，但 Python 和 SciPy 为采用它们奠定了基础，其中部分原因是 IPython/Jupyter 笔记本的吸引力。Python 生态系统让科学家们看到了自由软件作为研究工具的优势。使用自由软件不仅意味着避免了许可费（licensing fee，这在世界许多地方都是真正的障碍），而且还避免了黑箱，从而使得结果的开放性、可靠性、可重复性和信心都得到了更好的支撑。这对科学的好处是巨大的。

如果一个人有理由要使用 Python 生态系统——也许是需要一个特定库的特殊功能，或者是已经用 Python 编程的成员之一了，那么 SciPy 将是一个显然的选择。但 SciPy 可能不是每个寻找免费软件工具的科学家的最佳选择。那些只需要一个免费、快速、复杂的程序来进行矩阵计算、方程求解和类似的任务，并且喜欢安装一个不依赖于其他部件的简单的程序的人，应该考虑使用 Octave，在 12 月份的时候 LWN 里就已经有过介绍。希望作为自己的纸笔计算的补充的数学家，可能不会用到 SciPy 中的大部分内容。他们最好安装 SymPy 或 Maxima，这两个软件除了可以进行符号操作外，还可以进行一些数值计算和绘图。

如果是一个正在开始新项目的科学家，并且不需要与使用 Python 的同事进行交互，应该认真考虑把 Julia 作为高级语言。之前介绍 Julia 的文章概述了 Julia 在科学界崛起的一些原因。Python 的主要吸引力——特别是在数据科学等领域——是它的一些被广泛使用且久经考验的函数库，但这并不应该使得我们放弃 Julia，毕竟 Julia 访问 Python 库几乎和从 Python 代码中使用这些库一样容易。

一个科学家只要有一台笔记本电脑和一个互联网连接，就可以安装这些工具，并逐个试用。如今可用于科学研究的免费软件的深度和广度，在十几年前是无法想象的。此外，可以将所有这些软件都安装上，完全可行，这样我们就可以保持总是使用最合适的工具来辅助我们的研究。

全文完
LWN 文章遵循 CC BY-SA 4.0 许可协议。

欢迎分享、转载及基于现有协议再创作～

长按下面二维码关注，关注 LWN 深度文章以及开源社区的各种新近言论～