[译]学习IPython进行交互式计算和数据可视化（六）

第五章：高性能并行计算

一个反复被提及的反对使用Python进行高性能数值计算的言论是这种语言是动态解释型的，速度太慢。一种编译型低级语言，如C，能提供比它快几个数量级的运算速度。我们在第三章——使用IPython进行数值计算中已经引入了向量化这一概念表示了对这种观点的反对。NumPy 数组的运算速度甚至可以和C一样快，因为低速的Python循环可以使用快速的C循环替代。尽管有时会出现一些复杂的算法不能进行向量化或很难向量化，幸运的是我们还有其他的解决方案而不用丢弃所有的Python代码用C重写。我们在本章中将会介绍一些这样的解决方案。

首先，一个能够被很好利用的是现代电脑都具有的多核心CPU。一个标准的Python进程正常情况下值运行在单个核心上，但是通过并行处理将任务分发到多个核心甚至多台计算机上是有可能实现的。在IPython中要实现这一功能非常简单。很少的几行代码就可以使用MPI。

另外一种流行的解决方案是先检查出Python代码中的运算时间瓶颈的那一部分，接着用C代码将其替换掉。一般情况下只有非常少的Python代码占用了算法的大部分运算时间，所以可以用C替换掉这一小部分而保留其余的大部分Python代码。Cython是一个让这件任务比它听起来要简单的多的外部包：它提供了Python的编译型超集，可以很容易的和Python代码集成在一起。在IPython使用尤其方便。

到了本章的末尾，我们将已讨论了：

  • 如何在IPython中把独立的函数分布到多个核心上
  • 怎样在IPython中简单地使用MPI
  • 怎样使用一个cell 魔法通过Cython将Python代码转换为C代码
  • 怎样在Cython中使用NumPy来对你代码的速度进行数量级的提升

交互式任务并行化

在这一节，我们将看到怎样在IPython中把独立的函数分布到多个CPU核心上。

使用Python进行并行化计算

Python的原生并行计算功能功能相当不尽人意。一个长期存在的问题是CPython实现了一个全局解释器锁（GIL），而这个GIL就像CPython文档中所说的：

"...一种阻止同时有多个线程执行Python字节码的互斥机制"

由于CPython的内存管理不是线程安全的，所以GIL是必须的，但是它的一个主要缺点是它会阻止CPython的多线程程序充分利用多个处理器的资源。

Python的GIL:感兴趣的读者可以在下面的参考文献找到关于Python GIL的更多信息：

如果NumPy是由ATLAS、MKL、等库编译的，那么其中的一些线性代数函数也许可以通过释放GIL来利用多核处理器的资源。除此之外，使用Python一般用进程替代线程进行分布式任务。由于进程并不共享相同的内存空间，所以一些进程内部通信需要实现，例如，使用Python的原生multiprocessing模块。一个更加强大也更复杂的解决方案是使用消息传递接口Message Passing Interface（MPI）。

对于这两种解决方案IPython都能很好兼容，这一节我们就会讨论一些有关于这方面的问题。它提供了一个强大而通用的并行化计算架构。多个IPython引擎可以同时运行在不同的核心或不同的计算机上。得助于平衡加载load balancing，独立的任务可以非常简单地均匀分发。数据可以在引擎间进行传递，这使得在IPython中实现复杂分布式算法成为可能。

并行化计算是一个非常有难度的主题，我们只会涉及到一些基础的层面。

将任务分发到多个核心

IPython的并行化计算功能是可扩展和高度定制的，但我们只会展示使用它们的一些简单方式。除此之外，我们将关注并行化计算的交互式应用，因为这才是IPython的精髓所在。

将代码分发到一台主机的多个核心需要进行下面的几个步骤：

  1. 启动多个IPython引擎（一般情况下每个处理器对应一个IPython引擎）
  2. 创建一个Client对象作为这些引擎的代理
  3. 使用client（客户端）来启动这些引擎上的任务并获取结果。

任务既可以同步启动，也可以异步启动：

  1. 若采用同步任务（阻塞式），客户端在任务启动之后会立即阻塞，并在任务完成后返回任务结果。
  2. 若采用异步任务（非阻塞式），任务启动后会立即返回一个ASyncResult对象。这个对象会异步地记录任务的状态，任务完成后可以在任意时刻获取任务结果。

开启引擎

开启引擎的最简单方法是使用系统shell命令——ipcluster start进行启动调用。默认状态下，这个命令将会在本地的每一个核心上启动一个引擎。引擎的数目可以通过-n选项进行指定，例如，ipcluster start -n 2将会启动两个引擎。你可以通过ipcluster -h和ipcluster -start -h来查看所有的可用选项。除此之外，notebook有一个名为Cluster的表盘，提供了启动和停止引擎的web接口。

创建一个客户端实例

客户端是用来想引擎发送任务的。在一个IPython控制台或notebook中，我们首先需要从 parallel子包中导入Client类：

In [1]: from IPython.parallel import Client

接着就是创建Client实例：

In [2]: rc = Client()

IPython会自动的探测正在运行的引擎。我们可以向下面这样检查正在运行的引擎数目：

In [3]: rc.ids
Out[3]: [0, 1]

客户端的ids属性可以提供正在运行的引擎的标识符。在这里我们可以看到在本地主机上有两个正在运行的引擎（主机有双核处理单元）。

使用并行魔法

从IPython中想引擎发送任务的最简单方法是使用%px 魔法命令。

他会在引擎上运行一条命令：

In [4]: import os
In [5]: %px print(os.getpid())
[stdout:0] 6224
[stdout:1] 3736

默认状态下，这个命令将会以同步模式在所有正在运行的引擎上执行。有多种方式可以用来指定哪些引擎来执行。

第一种方式是使用%pxconfig魔法命令：

In [6]: %pxconfig --targets 1
In [7]: %px print(os.getpid())
3736

--target 选项接受一个索引或切片对象，例如：::2表示所有偶数索引指示的引擎。这里我们只是制定了第二个引擎（因为总计就只有两个引擎）。随后所有对%px的调用将会在指定的目标引擎上执行。

与其等价的一种方法是使用%%px小格魔法命令：

In [8]: %%px --targets :-1
print(os.getpid())
[stdout:0] 6224

%%px选项将会应用在整个小格，这在notebook环境中用起来非常方便。

另外一个可用的选项是阻塞模式(blocking mode)。默认状态下%px魔法命令使用的是阻塞模式。我们可以使用--noblock选项来使费阻塞模式(non-blocking mode)生效。

In [9]: %%px --noblock
import time
time.sleep(1)
os.getpid()
Out[9]: <AsyncResult: execute>

任务接着就会异步执行。%pxresult魔法命令会阻塞解释器直到任务结束返回结果。

In [10]: %pxresult
Out[1:12]: 3736

并行映射

Python自带的map函数会将元素逐个排列起来。IPython提供一个语义等价但能将不同任务分发给不同引擎的并行化map函数。使用这个函数试讲任务分布到多个核心的最简单的方法。

创建视图

要使用它，我们首先需要使用Client实例获得一个引擎视图。一个视图（view）表示一个或多个引擎，同时还带有客户端的索引语法。例如，我们可以使用下面的命令获得一个包含所有引擎的视图：


In [11]: v = rc[:]

这个视图接着就可以用来在引擎上启动任务。而且我们可以使用sync_imports()方法在引擎上导入包：

In [12]: with v.sync_imports():
import time
importing time on engine(s)

同步映射

让我们来定义一个下面这样的函数：

In [13]: def f(x):
			time.sleep(1)
			return x * x

这个函数接受一个数字参数并在等待一秒后返回它的数字的平方。我们可以使用v.map_sync()这个方法将0-9的数字分别传递给这个函数并使用我们的两个引擎（也就是两个CPU）同步执行：

In [14]: v.map_sync(f, range(10))
Out[14]: [0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

几秒钟之后我们获得了一个结果列表。这里每一个引擎都处理了5个任务，总计处理了10个任务：


In [15]: %timeit -n 1 -r 1 v.map_sync(f, range(10))
1 loops, best of 1: 5.02 s per loop
In [16]: %timeit -n 1 -r 1 map(f, range(10))
1 loops, best of 1: 10 s per loop

异步映射

我们可以使用v.map()方法根据参数列表来异步执行函数。：

In [17]: r = v.map(f, range(10))
In [18]: r.ready(), r.elapsed
Out[18]: False, 2.135
In [19]: r.get()
Out[19]: [0, 1, 4, 9, 16, 25, 36, 49, 64, 81]
In [20]: r.elapsed, r.serial_time
Out[20]: (5.023, 10.008)

r变量是一个ASyncResult对象，该对象包含多种可以用来储蓄进度信息、已用时间和获取任务结果的属性和方法。elapsed
属性可以在任务开始后的任意时间返回从任务开始已经使用的时间。serial_time属性只有在任务完成后才能使用，能返回所有隐情在每一个任务上累计花费的时间。ready()方法在任意时刻都可以返回一个表示任务是否已经完成的数值。get()方法会一直处于阻塞状态直到任务结束返回任务结果。

一个实际示例——蒙特卡罗模拟

为了展示IPython提供的并行计算功能，我们考虑进行一个新的示例。我们来使用蒙特卡罗模拟来估测P一下Pi常量值（π）。原理是如果在一个边为1的正方内有n个随机取样的点，在n趋向于无穷大时，距离正方形一个指定顶点小于1的点所占所有点的比例接近于π/4。下面的图像描绘了实际效果：
【图】
这是一个蒙特卡罗模拟的详细实例，我们将多次设置大量的随机数，and takes an average at the end to estimate
some quantity of interest that would be difficult to obtain with a deterministic
method. 蒙特卡罗模拟被广泛应用于科学、工程和金融领域。由于只是多次运行相同的独立函数，所以很容易进行并行化。

这里我们将使用这个随机试验来估计π的值。虽然我们已经知道通过这种方法获得的π值精度非常低，有很多方法比它更高效且精度更高。但是用这个示例来介绍IPython的并行化计算功能已经足够了。

首先，我们先编写执行这个模拟过程的Python代码。sample函数将会生成n个横纵坐标均在0到1之间随机点，并返回在四分之范围内的点的数目。

In [1]: def sample(n):
return (rand(n) ** 2 + rand(n) ** 2 <= 1).sum()

由于圆括号内的长度为n的响亮是一个掩码数组（也就是它包含的是布尔值），所以他的和是真（True）值的数目，也就是说耳机里的距离在0到1之间的点的数目。

In [2]: n = 1000000.
In [3]: 4 * sample(n) / n
Out[3]: 3.142184

由于π的真实值是3.1415926535...，我们可以看到使用10W个点进行模拟时小数点后只有两位是正确的。现在我们将这些任务分布到多个核心。假设已经启动了多个引擎（比如使用ipcluster start），下面就是怎样使这些代码并行化：

In [4]: from IPython.parallel import Client
			rc = Client()
			v = rc[:]
			with v.sync_imports():
				from numpy.random import rand
In [5]: 4 * sum(v.map_sync(sample, [n] * len(v))) / (n * len(v))
Out[5]: 3.141353

这里,len(v)表示引擎的数目。我们使用相同的参数n调用sample函数len(v)次。所有结果的和是红色点的总数，而总的点数是n*len(v)。最后我们使用先前的公式得到了π的估计值。

在IPython中使用MPI

MPI是一个有名的标准化消息传递系统，对于并行化计算来说效率非常高。我们假设你的系统已经安装好了MPI实现（如Open-MPI,http://www.open-mpi.org）和相应的Python封装——mpi4py(http://mpi4py.scipy.org）。关于如何安装MPI的信息可以在上面的网站中找到。

Windows上的MPI:如果你使用的是Windows操作系统，可以使用HPC Pack中微软的MPI实现（http://
www.microsoft.com/en-us/download/details.
aspx?id=36045）。你也许还会对Visual Studio中的Python工具感兴趣（http://pytools.codeplex.com），其中的工具可以将Visual Studio设置成一个Python IDE。这些工具为IPython提供了原生支持，且已经针对使用MPI进行高性能计算进行了专门的设计。

首先，我们需要为MPI创建一个专门的Ipython配置文件。在shell中输入以下命令：

python profile create --parallel --profile=mpi

接着，编辑IPYTHONDIR/profile_mpi/ipcluster_config.py文件（IPYTHONDIR一般情况下是~/.ipython）将下面的这一行添加进去：

c.IPClusterEngines.engine_launcher_class = 'MPIEngineSetLauncher'

现在，输入下面的命令来启动四个引擎：

ipcluster start -n 4 --profile=mpi

为了在IPython中使用MPI，我们首先需要编写一个使用mpi4py封装的MPI的函数。在这个示例中，我们将在四个核心上并行化计算1到16之间所有整数的和。让我们在一个名为psum.py的文件中写入下面的代码：

from mpi4py import MPI
import numpy as np
# This function will be executed on all processes.
def psum(a):
	# "a" only contains a subset of all integers.
	# They are summed locally on this process. 
	locsum = np.sum(a)
	# We allocate a variable that will contain the final result, 
	that 
# is the sum of all our integers.
	rcvBuf = np.array(0.0,'d')
# We use a MPI reduce operation:
# * locsum is combined from all processes
# * these local sums are summed with the MPI.SUM operation
# * the result (total sum) is distributed back to all processes 
in
# the rcvBuf variable
	MPI.COMM_WORLD.Allreduce([locsum, MPI.DOUBLE],
		[rcvBuf, MPI.DOUBLE],
		op=MPI.SUM)
	return rcvBuf

最后，我们就可以像下面这样在IPython中交互式地使用这个函数了：


In [1]: from IPython.parallel import Client
In [2]: c = Client(profile='mpi')
In [3]: view = c[:]
In [4]: view.activate() # enable magics
In [5]: view.run('psum.py') # the script is run on all processes
In [6]: view.scatter('a', np.arange(16)) # this array is scattered across 
processes
In [7]: %px totalsum = psum(a) # psum is executed on all processes
Parallel execution on engines: [0, 1, 2, 3]
In [8]: view['totalsum']
Out[8]: [120.0, 120.0, 120.0, 120.0]

更多关于如何在IPython中使用MPI的细节信息可以在下面IPython官方文档网页上找到(我们这个示例就来自这里)：http://ipython.org/ipython-doc/stable/parallel/parallel_mpi.html。

IPython的更多高级并行计算功能

我们只涉及了IPython中可用的并行计算功能中非常基础的部分。更多的高级功能如下：

  • 动态平衡加载
  • 引擎间的pushing和pulling对象
  • 使用SSH隧道将引擎运行在不同的主机上
  • 在带有StarCluster的Amazon EC2集群上使用IPython（译者注：Amazon EC2有免费的一年使用权限，前提是你有一张支持美元货币的信用卡）
  • 将所有的请求和结果存储到数据库中
  • 使用有向无环图（DAG）管理任务依赖

这些功能远远超出了本书的范围。感兴趣的读者可以在IPython的官方文档上找到所有这些功能的细节。

基于Cython在IPython中使用C

将独立的任务分布到多个核心的最简单方式是充分利用现代计算机的并行能力从而减少两倍或更多的运行时间。然而，有些算法不能很容易地分解为独立的子任务。除此之外，也许会遇到这样的情况：由于牵涉到了无法矢量化的嵌套循环，python代码编写的算法本身速度非常慢。在这种情况下，一个非常有趣的选择是可以将其中少量的关键代码用C替换掉从而大大地降低Python消耗的时间。这个方案不牵涉任何的并行计算功能，但依然能极大的提高Python脚本的运行效率。除此之外，两种技术是可以同时使用的：局部代码使用C编写，并在IPython中进行并行计算。

Cython可以将一部分设置没有明确以C的形式进行转换的Python代码进行编译；它提供了一个Python扩展语法来调用C函数和定义C数据类型。影响性能的问题代码接着就会被自动的转换为C代码、编译，然后就可以在Python中显式使用了。某些情况下仅仅只用纯Python代码是可以的；有些情况下由于算法的实际特性，速度的提升是十分猛烈的，甚至能达到及格数量级的提示，使用umPy进行向量化也达不到。

在这一节中，我们将看到怎样在IPython中交互式地使用Cython。回回看到一个势力，在这个示例中我们用纯Python代码实现一个算法，对其进行简单地修改，运算速度就有了超过300倍的提升。

安装和配置Cython

相较于其它包，Cython安装起来稍微难一些。原因是使用Cython意味着编译C代码，显然这需要一个C语言编译器（例如非常流行的GNU C编译器 gcc）。在Linux系统上，gcc要么已经默认安装了，要么通过包管理器也能很容易地安装好，例如Ubuntu和Debian发行版的sudo apt-get install构建工具。在OS X系统上，我们可以安装Apple XCode。在Windows上你可以安装开源的gcc发型包——MinGW(http://www.mingw.
org)。接着Cython就可以像其他包一样安装了（查看第一章——开始使用IPython）。更多相关的信息可以在这里找到：http://wiki.cython.org/Installing。

在Windows上配置MinGW 和 Cython：在Windows上，由于MinGW的版本原因，在编译Cython的时候可能会有错误消息出现。要解决这个BUG，你也许需要打开C:\Python27\Lib\distutils\cygwinccompiler.py（或者相似的路径，这取决于你自己的配置）将其中的-mno-cygwin对应的资源指引用空字符串""替换掉。
同事，你要确认C:\MinGW\binis 在你的系统环境PATH变量内。最后，你也许需要编辑（或创建）C:\Python27\Lib\distutils\distutils.cfg这个文件并将下面的几行代码添加进去：

[build]
compiler = mingw32

你可以在下面的网站找到更多的相关信息：http://wiki.
cython.org/InstallingOnWindows。

在IPython中使用Cython

在使用Cython时，一般将代码脚本以.pyx作为后缀，这表示这些文件要使用Cython编译为C语言。接着，生成的C程序会被C编译器编译成.so文件（在Linux上）或.pyd文件（在Windows上），这些文件可以正常被Python代码导入。

这个过程一般会涉及到一个用来设置要编译的文件以及不同的编译选项的distutils setup.py脚本。因为这一步不是非常难，我们在这里就不做叙述了。接着我们将展示怎样在IPython中使用Cython。这里有一个有利条件是Cython进行C语言编译时在后台自动进行的，而不需要一个手工编写的setup.py脚本。IPython notebook在这个地方非常有用，因为相较于控制台，在它上面可以很方便的编写多行代码。

这里我们将展示怎样使用%%cython小格魔法来从IPython执行Cython。第一步是加载cythonmagic扩展。

In [1]: %load_ext cythonmagic

接着，%%cython小格魔法就允许编写可以被自动编译的Cython代码了。小格内定义的函数在交互式会话中成为可用的。Python代码就可以正常的调用他了。

In [2]: %%cython
	def square(x):
		return x * x
In [3]: square(10)
Out[3]: 100

在这里，对`suqare(10)的调用牵涉到了对已编译的用来求数字平方的C语言函数的调用。

通过Cython对纯Python算法进行加速

这里我们将会看到一个牵涉到嵌套循环的纯Python算法怎样在Cython中尽心转化来达到10倍的速度提升。这个算法是爱拉托逊斯筛法，该算法用来寻找小于某一给定数值的所有素数。这是一个非常经典的算法，算法从2到n之间的任意一个整数开始，并逐渐移除已经找到的素数的倍数。算法运行到最后，只有素数会被留下来。我们将用Python实现这个算法，并展示怎样才能在Cython中将其转化。

纯Python版本

纯Python编写的算法有10几行长。有很多方法可以对这个实现进行改进和缩短（存在一个只有一行代码的算法！），但是因为我们只关心纯Python版本和Cython版本之间的相对执行时间，所以这个版本已经足够用了。

In [1]: def primes1(n):
			primes = [False, False] + [True] * (n - 2)
			i = 2
# The exact code from here to the end of the function
# will be referred as #SIEVE# in the next examples.
			while i < n:
# we do not deal with composite numbers
				if not primes[i]:
					i += 1
					continue
				k = i * i
# mark multiples of i as composite numbers
				while k < n:
					primes[k] = False
					k += i
				i += 1
			return [i for i in xrange(2, n) if primes[i]]
In [2]: primes(20)
Out[2]: [2, 3, 5, 7, 11, 13, 17, 19]

变量primes包含用来表示对应索引的数是不是素数的布尔值。根据“若一个正整数只能被两个数整除，则这个正整数为素数”这一定义，我们只用0和1这两个合数来对它进行初始化。接着，在每一次迭代中，我们将在不改变素数的情况下标记越来越多的合数。每一个i表示一个素数，k的迭代用来把所有i的倍数标记为合数。最后我们返回结果为真的索引列表，也就是说列表包含了所有小于n的素数。

现在，让我们卡一下这个函数的执行时间：

In [3]: n = 10000
In [4]: %timeit primes1(n)
100 loops, best of 3: 5.54 ms per loop

我们将尝试使用Cython对这个函数进行加速。

幼稚的Cython转换

作为第一次尝试，我们将在Cython中使用完全相同的代码：

In [5]: %load_ext cythonmagic
In [6]: %%cython
		def primes2(n):
			primes = [False, False] + [True] * (n - 2)
			i = 2
#SIEVE#: see full code above
In [7]: timeit primes2(n)
100 loops, best of 3: 3.25 ms per loop

我们仅仅是在小格的头部添加了了%%cython就获得了70%的速度提升。但是我们可以通过给Cython提供类型信息来做的更好。

添加C语言类型

由于本地变量是动态类型的Python变量，先前实例中的速度提升实在是太逊了。这说明Python和纯C语言代码之间的性能差异与Python的动态特性有很大关系。我们可以使用cdef官架子将Python变量转换成C语言变量来提升性能：

In [8]: %%cython
	def primes3(int n):
		primes = [False, False] + [True] * (n - 2)
		cdef int i = 2
		cdef int k = 0
	#SIEVE#: see full code above

与先前的幼稚版本相比，共有三处改动：参数n被固定的声明为一个整数，局部变量i和k被声明为C语言中的整数变量。这样做有多少速度提升呢：

In [9]: timeit primes3(n)
1000 loops, best of 3: 538 us per loop

这个函数现在比纯Python版本快了十倍，而我们做的只是使用了一个%%cython魔法和一些类型声明。这个结果通过更为合适的数据结构甚至还能改进。

一般情况下，要想知道可以用Cython进行转换作为主要速度提升的那一部分代码需要一些Python内部构造相关的知识，更重要的是要做广泛的性能分析。Python循环（尤其是潜逃循环），Python函数调用、密集循环中的高级数据结构处理是Cython一般的优化目标。

（同时)使用NumPy和Cython

在这一节中，我们见展示怎样把NumPy数组和Cython代码集成到一起。我们还会看到密集循环中的Python函数调用怎样才能通过将Python函数转换成C语言函数来对性能进行极大的优化。

Python版本

在这里我们见使用一个随机过程模拟（或者可以称为布朗运动）的示例。这个过程描述了一个质点的运动轨迹，质点从x=0处开始运动，每一个极短的dx时间间隔内进行随机的+dx或-dx的步长运动。这种类型的过程会频繁地出现在金融、经济、物理、生物等领域。

NumPy的cumsun()和rand()函数可以对这种特定的过程进行高效的模拟。然而，还有很多更为复杂的过程需要模拟，比如，某些模型需要当位置到达某一临界值时出现瞬时跳跃。在这种情况下，向量化不不能实现，因此手动循环是不可避免的。

In [1]: def step():
			return sign(rand(1) - .5)
		def sim1(n):
			x = zeros(n)
			dx = 1./n
			for i in xrange(n - 1):
				x[i+1] = x[i] + dx * step()
			return x

step函数通过NumPy的sign()函数和rand()函数来随机返回+1或-1。在sim1()函数中，轨迹被初始化为一个元素全为零的NumPy向量。接着，每一次迭代，一个写的步长就会被添加到古籍中。then()函数返回全部的轨迹。下面是一个轨迹的示例：

In [2]: plot(sim1(10000))

让我们来看一下这个函数的执行时间。


In [3]: n = 10000
In [4]: timeit sim1(n)
1 loops, best of 3: 249 ms per loop

Cython版本

对于Cython版本，我们将做两件事。第一，我们将对所有的局部变量连同包含运动轨迹的NumPy 数组添加C语言类型。第二，我们将把step()函数转换成不许调用任何NumPy函数的纯C函数。我们将调用以C标准库形式定义的纯C函数。

In [4]: %%cython
	import numpy as np
	cimport numpy as np
	DTYPE = np.double
	ctypedef np.double_t DTYPE_t
	# We redefine step() as a pure C function, using only
# the C standard library.
	from libc.stdlib cimport rand, RAND_MAX
	from libc.math cimport round
	cdef double step():
		return 2 * round(float(rand()) / RAND_MAX) - 1
	def sim2(int n):
# Local variables should be defined as C variables.
		cdef int i
		cdef double dx = 1. / n
		cdef np.ndarray[DTYPE_t, ndim=1] x = np.zeros(n, 
		dtype=DTYPE)
		for i in range(n - 1):
			x[i+1] = x[i] + dx * step()
		return x

我们首先需要导入标准的NumPy库和一个特殊的也称作NumPy的C语言库。这个特殊的C语言库是Cython包的一部分，用cimport进行导入，我们使用ctypedef将NumPy数据类型double定义成对应的C语言数据类型double_t。它在编译时定义x数组的确切类型而不是执行时，这就带来了极大的速度提示。x的维度数也在sim2()函数中被指定。所有的局部变量均定义为C数据类型的C语言变量。

step()函数已经被完全重写。它现在是一个纯C函数（以cdef定义的）。它使用C标准库中的rand()函数，这个函数返回一个介于0到RAND_MAX之间的随机数。math库中的round()函数被用来生成+1或-1这两个随机数。

让我们来检查一下sim2()函数的执行时间：


In [5]: timeit sim2(n)
1000 loops, best of 3: 670 us per loop

Cython版本比Python版本快370倍。速度提升如此夸张的主要原因是Cython版本只是用了C语言代码。所有的变量都是C语言变量，先前需要花费大量时间的Python step函数的调用现在只牵涉对纯C函数的调用，因此极大的减少了循环中Python的时间消耗。

更高级的Python代码加速方法

Cython也可以用来将已存在的C语言代码和Python连接起来，但是我们不会在这里涉及这一部分。

除了Cython，还有其它的用来进行Python代码加速的包。SciPy.weave(http://www.scipy.org/Weave)是SciPy的子包可以用来将C/C++代码嵌入到Python代码中.Numba(http://numba.pydata.org/) 使用即时LLVM编译通过通过动态显示地编译Python代码来达到速度提升的效果。它能很好的和NumPy集成在一起。他的安装需要llvmpy和meta。

相关的工程包括:Theano(http://deeplearning.net/software/theano/),通过在CPU或显卡上显式编译Python代码来高效地定义、优化和评估数学数组表达式。与此相似，Numexpr(https://code.google.com/p/numexpr/)可以编译数组表达式和充分利用矢量CPU指令和多核心处理器。

Blaze (http://blaze.pydata.org/)在本书写作的时候是仍是一个处于早期发展状态的项目。该项目的目标是将所有的这些动态编译技术整合到一个统一的框架中。而且该项目通过允许数组类型和形状不一致、数值缺失、标识尺寸（就像Pandas中那样）等特性扩展了多维数组的概念。该项目现在由NumPy的作者进行开发，在不久的将来，也似乎会成为Python科学计算社区的中心项目。

最后，PyOpenCL(http://mathema.tician.de/software/pyopencl) 和PyCUDA (http://mathema.tician.de/software/pycuda)是OpenCL和CUDA的Python封装。这些库实现了能在现代显卡上编译的类C低级语言，以充分利用显卡的大规模并行架构。实际上显卡包含成百上千的特殊核心，所以能够在大量处理单元按上进行高效的函数处理（单指令多数据流范式[SIMD]）。与纯C代码相比速度提升甚至超过了一个数量级。OpenCL是一个开源的标准语言，CUDA是一个由Nvidia集团持有的专利语言。CUDA代码只能运行在Nvidia显卡上，而大多数显卡和CPU都支持OpenCL。对于OpenCL，相同的代码可以在CPU上编译以充分利用多核心处理器和矢量指令的优势资源。

总结

在本章中，我们介绍了两种Python代码加速方法：一是通过将Python代码转换为低级的C语言代码绕过Python语言的性能瓶颈；二是将Python代码分发到多个计算单元上以充分利用多核心处理器的资源。两种方法甚至可以同时使用。IPython极大地简化了这些技术。没有IPython也能并行化计算和Cython也能使用，但那将会需要更多的引用代码。

在下一章中，我们将探索一些高级选项来定制IPythonm。

posted @ 2016-04-01 21:42  Tacey Wong  阅读(1718)  评论(0)  编辑  收藏  举报