使用 Python 进行科学计算

作者:Joey Bernard

随着计算机变得越来越强大,科学计算正成为深入研究我们世界运作方式的基础研究中越来越重要的组成部分。我们现在能做的事情比十年前甚至无法想象的还要多。

这项工作的大部分传统上都是用更底层的语言完成的,例如 C 或 FORTRAN。最初,这样做是为了最大限度地提高代码的效率,并从计算机中榨取最后一丝性能。随着计算机现在达到多 GHz 的速度,这不再是曾经的瓶颈。其他效率因素开始发挥作用,其中程序员的效率至关重要。考虑到这一点,人们正在考虑其他语言,这些语言有助于最大限度地利用程序员的时间和精力。

本文讨论了其中一个选择:Python。尽管 Python 是一种解释型语言,并且不公正地遭受由此带来的污名,但它因其清晰的风格和许多有用的软件包的可用性而在科学家中越来越受欢迎。我在本文中特别关注的软件包旨在提供快速、强大的数学和科学工具,这些工具的运行速度几乎与 C 或 FORTRAN 代码一样快。

设置

我在这里关注的软件包称为 numpy 和 scipy。它们都可从主要的 SciPy 站点获得(参见“资源”)。但在我们下载它们之前,numpy 和 scipy 到底是什么?

numpy 是一个 Python 包,它提供了扩展的数学功能。这些功能包括新的数据类型,例如无限大小的长整数和复数。它还提供了一种新的数组数据类型,允许构建向量和矩阵。应用于这些新数据类型的所有基本操作也包含在内。有了这个,我们已经可以完成相当多的科学工作了。

scipy 是构建在 numpy 之上的进一步扩展。这个软件包简化了许多需要处理的更常见的任务,包括用于查找多项式根、进行傅里叶变换、进行数值积分和增强 I/O 的工具。有了这些功能,用户可以在相对较短的时间内开发非常复杂的科学应用程序。

现在我们知道了 numpy 和 scipy 是什么,我们如何获得它们并开始使用它们呢?大多数发行版都包含这两个软件包,这使得安装它们变得容易。只需使用您的发行版的包管理器进行安装即可。例如,在 Ubuntu 中,您将在终端窗口中键入以下内容

sudo apt-get install python-scipy

这将安装 scipy 及其所有依赖项。

如果您想使用最新版本,并且不想等待您的发行版更新,则可以通过 Subversion 获得它们。只需执行以下操作

svn co http://svn.scipy.org/svn/numpy/trunk numpy svn co
http://svn.scipy.org/svn/scipy/trunk scipy

构建和安装由源目录中的 setup.py 脚本处理。对于大多数人来说,构建和安装只需要

python setup.py build
python setup.py install    # done as root

如果您没有 root 访问权限,或者不想安装到系统软件包目录中,则可以使用以下命令安装到不同的目录中

python setup.py install --prefix=/path/to/install/dir

还有其他选项可用,您可以通过使用以下命令了解这些选项

python setup.py --help-commands

花时间进行实验,看看您是否可以在特定情况下使用任何额外的选项。

基本数学

现在我们已经安装了 scipy 和 numpy,让我们开始我们的旅程,看看科学计算中经常使用的一些基本函数。最常见的任务之一是矩阵数学。当您使用 numpy 时,这大大简化了。使用 numpy 进行两个矩阵乘法的最基本代码如下所示

import numpy
a1=numpy.empty((500,500))
a2=numpy.empty((500,500))
a3=a1*a2

将此与我们在 C 中编写的代码进行对比

#include <stdlib.h>
int main() {
   double a1[500][500];
   double a2[500][500];
   double a3[500][500];
   int i, j, k;
   for (i=0; i<500; i++) {
      for (j=0; j<500; j++) {
         a3[i][j] = 0;
         for (k=0; k<500; k++) {
            a3[i][j] += a1[i][k] * a2[k][j];
         }
      }
   }
}

Python 代码更短更简洁,代码的意图也更清晰。代码的这种清晰性意味着程序员可以更多地关注算法,而不是实现的细节。有一些 C 库,例如 LAPACK,可以帮助简化 C 中的这项工作。但是,即使这些库也无法与 scipy 的简洁性相媲美。

“但是效率呢?”,我听到你问。好吧,让我们用一些定时运行来看看它。以我们上面的例子为例,我们可以在实际的矩阵乘法部分周围放置一些调用,看看每个调用需要多长时间。有关结果,请参见表 1。

表 1. 平均运行时间

语言平均时间(秒)
C1.620
C (-O3)0.010
Python0.250

虽然您的结果可能会有所不同,因为这些时间取决于您的硬件以及您的机器上还在运行的其他程序,但我们可以看到一个总体趋势。Python 代码实际上比没有命令行选项编译的 C 代码快大约八倍。这实际上非常令人惊讶。一旦我们使用优化命令行选项,我们就会看到 C 代码现在更快了,大约快了 25 倍。因此,我们可以使用优化的 C 代码获得更快的代码,但我们需要意识到,在四分之一秒内将两个各有 250,000 个元素的矩阵相乘可能已经足够快了。

此外,当我们使用 Python 时,我们获得了一定程度的保护。如果我们尝试将两个矩阵相乘,而这种乘法在数学上没有意义,会发生什么?当我们尝试将两个不同大小的矩阵相乘时,Python 会给出我们

ValueError: shape mismatch: objects cannot be 
 broadcast to a single shape

在 C 中,我们根本没有收到任何错误。这是因为当我们处理矩阵时,我们实际上是在使用指针算术。因此,我们所做的几乎任何事情对于 C 都是有效的,即使它在问题域中没有意义。

我们也可以同样轻松地处理复数。如果我们想创建一个 64 位复数数组,我们可以这样写

a=zeros((500,500), dtype=complex64)

这将为我们提供一个用零初始化的 500x500 元素矩阵。我们使用以下命令访问每个元素的实部和虚部

a.real[0,0]=1.0 a.imag[0,0]=2.0

这将值 1+2j 设置到 [0,0] 元素中。

还有一些函数可以为我们提供更复杂的结果。这些包括点积、内积、外积、逆矩阵、转置、迹等等。不用说,我们已经掌握了大量的工具来完成相当多的科学工作了。但这足够了吗?当然不够。

开始“真正”的科学

现在我们可以做一些数学运算了,我们如何完成一些“真正”的科学工作呢?这就是我们开始使用我们感兴趣的第二个软件包 scipy 的功能的地方。有了这个软件包,我们有更多的函数可用于完成一些相当复杂的计算科学。让我们看一个简单的数据分析示例,以展示可能完成的工作类型。

假设您收集了一些数据,并且想看看这些数据的形式,是否存在任何周期性。以下代码使我们能够做到这一点

import scipy
inFile = file('input.txt', r)
inArray = scipy.io.read_array(inFile)
outArray = fft(inArray)
outFile = file('output.txt', w)
scipy.io.write_array(outFile, outArray)

如您所见,读取数据是一行代码。在本例中,我们使用 FFT 函数将信号转换为频域。这使我们能够看到数据中频率的分布。等效的 C 或 FORTRAN 代码太长,无法在此处包含。

但是,如果我们想查看这些数据以查看是否有任何有趣的东西呢?幸运的是,还有另一个名为 matplotlib 的软件包,可以用于生成用于此目的的图形。如果我们生成一个正弦波并将其传递给 FFT,我们可以通过绘制它来查看此数据的形式(图 1 和图 2)。

Use Python for Scientific Computing

图 1. 正弦波

Use Python for Scientific Computing

图 2. 正弦波的 FFT

我们看到正弦波看起来很规则,FFT 通过在正弦波频率处有一个单峰来证实这一点。我们刚刚做了一些基本的数据分析。

这向我们展示了进行相当复杂的科学编程是多么容易。而且,如果我们使用交互式 Python 环境,我们可以以探索性的方式进行这种科学分析,从而使我们能够在接近实时的状态下对数据进行实验。

对我们来说幸运的是,SciPy 项目的人们已经考虑到了这一点,并为我们提供了程序 ipython。这也可以在主要的 SciPy 站点上找到。ipython 的编写目的是以非常无缝的方式与 scipy、numpy 和 matplotlib 一起工作。要使用 matplotlib 支持执行它,请键入

ipython -pylab

该界面是一个简单的 ASCII 界面,如图 3 所示。

Use Python for Scientific Computing

图 3. ipython 窗口

如果我们使用它来绘制上面的正弦波,它只会弹出一个显示窗口来绘制绘图(图 4)。

绘图窗口允许您保存您出色的图形和绘图,这样您就可以向全世界展示您的科学突破。本文的所有绘图实际上都是以这种方式生成的。

Use Python for Scientific Computing

图 4. ipython 绘图

因此,我们已经开始进行一些真正的计算科学和一些基本的数据分析。我们接下来做什么?当然,我们会做得更大。

并行化

到目前为止,我们已经研究了相对较小的数据集和相对简单的计算。但是,如果我们有大量的数据,或者我们有更复杂的分析想要运行怎么办?我们可以利用并行性,并在高性能计算集群上运行我们的代码。

SciPy 站点的优秀人员编写了另一个名为 mpi4py 的模块。此模块提供了 MPI 标准的 Python 实现。有了它,我们可以编写消息传递程序。但是,安装它确实需要一些工作。

第一步是在您的机器上安装 MPI 实现(例如 MPICH、OpenMPI 或 LAM)。大多数发行版都有 MPI 的软件包,因此这是安装它的最简单方法。然后,您可以按照通常的方式使用以下命令构建和安装 mpi4py

python setup.py build python setup.py install

要测试它,请执行

mpirun -np 5 python tests/helloworld.py

这将运行一个五进程 Python 并运行测试脚本。

现在,如果我们的瓶颈在于划分大型数据集,我们可以编写程序在可用处理器之间划分大型数据集。或者,如果我们想进行大型模拟,我们可以将模拟空间划分到所有可用处理器之间。不幸的是,对 MPI 编程的有用讨论将是另一篇或两篇文章的内容。但是,我鼓励您获得一本关于 MPI 的好教科书,并亲自进行一些实验。

结论

尽管任何解释型语言都很难与编译优化的语言的速度相媲美,但我们已经看到,这不再像以前那样是一个巨大的障碍。现代机器的运行速度足够快,足以弥补解释的开销。这为使用像 Python 这样的语言打开了复杂应用程序的世界。

本文只能介绍最基本的可用的功能。幸运的是,已经编写了许多非常好的教程,可以从主要的 SciPy 站点获得。因此,走出去,以 Python 的方式做更多的科学研究吧。

ScientificPython

numpy 和 scipy 不是 Python 程序员可用的唯一选择。另一个流行的软件包是 ScientificPython。它包括几何类型(例如向量、张量和四元数)、多项式、基本统计信息、导数、插值等等。这与 scipy 中提供的功能类型相同。主要区别在于 ScientificPython 具有内置的并行编程能力,而 scipy 则需要额外的模块。这是通过 MPI 的部分实现和 Bulk Synchronous Parallel 库 (BSPlib) 的实现来完成的。

LAPACK 和 BLAS

可以提出的论点是,将 C 和 FORTRAN 的复杂性与 Python 的复杂性进行比较是不公平的,因为我们实际上是在 Python 中使用附加软件包。等效的库可以在 C 和 FORTRAN 中使用,其中 LAPACK 和 BLAS 是一些更流行的库。BLAS 提供基本的线性代数函数,而 LAPACK 基于这些函数提供更复杂的科学函数。虽然这些库提供了优化的例程,可以从您的硬件中提取每个有用的周期,并且比直接编写 C 或 FORTRAN 简单得多,但它们仍然比 Python 中的等效例程复杂几个数量级。但是,如果您真的需要从您的机器中榨取最后一丝性能,那么没有什么能比得上这些类型的库。

并行编程的类型

一般来说,并行程序可以分为两大类:共享内存和消息传递。在共享内存并行编程中,代码在一台物理机器上运行,并使用多个处理器。这种类型的并行编程的示例包括 POSIX 线程和 OpenMP。这种类型的并行代码受限于您可以构建的机器的大小。

为了绕过此限制,您可以使用消息传递并行代码。在这种形式中,独立的执行单元通过来回传递消息进行通信。这意味着它们可以位于不同的机器上,只要它们具有某种通信方式即可。这种类型的并行编程的示例包括 MPICH 和 OpenMPI。大多数科学应用程序都使用消息传递来实现并行性。

资源

Python 编程语言 — 官方网站:www.python.org

SciPy: www.scipy.org

ScientificPython — 理论生物物理学、分子模拟和数值密集型计算:dirac.cnrs-orleans.fr/plone/software/scientificpython

Joey Bernard 具有物理学和计算机科学的背景。最后,他在 ACEnet 的最新工作使他有机会同时使用这两个学位,帮助研究人员进行 HPC 工作。

加载 Disqus 评论