内容


可爱的 Python

Numerical Python

使用 Numeric 软件包和 Numarray 软件包

Comments

系列内容:

此内容是该系列 # 部分中的第 # 部分: 可爱的 Python

敬请期待该系列的后续内容。

此内容是该系列的一部分:可爱的 Python

敬请期待该系列的后续内容。

要了解 Numerical Python 软件包的第一件事情是,Numerical Python 不会让您去做标准 Python 不能完成的任何工作。它只是让您 以快得多的速度去完成标准 Python 能够完成的相同任务。实际上不仅仅如此;许多数组操作用 Numeric 或者 Numarray 来表达比起用标准 Python 数据类型和语法来表达要优雅得多。不过,惊人的速度才是吸引用户使用 Numerical Python 的主要原因。

其实,Numerical Python 只是实现了一个新的数据类型:数组。与可以包含不同类型元素的列表、元组和词典不同的是,Numarray 数组只能包含同一类型的数据。Numarray 数组的另一个优点是,它可以是多维的 -- 但是数组的维度与列表的简单嵌套稍有不同。Numerical Python 借鉴了程序员的实践经验(尤其是那些有科学计算背景的程序员,他们抽象出了 APL、FORTRAN、MATLAB 和 S 等语言中数组的最佳功能),创建了可以灵活改变形状和维度的数组。我们很快会回来继续这一话题。

在 Numerical Python 中对数组的操作是 按元素进行的。虽然二维数组与线性代数中的矩阵类似,但是对它们的操作 (比如乘) 与线性代数中的操作 (比如矩阵乘) 是完全不同的。

让我们来看一个关于上述问题的的具体例子。在纯 Python 中,您可以这样创建一个“二维列表”:

清单 1. Python 的嵌套数组
>>> pyarr = [[1,2,3],
...         [4,5,6],
...         [7,8,9]]
>>> print pyarr
[[1, 2, 3], [4, 5, 6], [7, 8, 9]]
>>> pyarr[1][1] = 0
>>> print pyarr
[[1, 2, 3], [4, 0, 6], [7, 8, 9]]

很好,但是您对这种结构所能做的只是通过单独的 (或者多维的) 索引来设置和检索元素。与此相比,Numarray 数组要更灵活:

清单 2. Numerical Python 数组
>>> from numarray import *
>>> numarr = array(pyarr)
>>> print numarr
[[1 2 3]
 [4 0 6]
 [7 8 9]]

改变并不大,但是使用 Numarray 进行的操作如何呢? 下面是一个例子:

清单 3. 元素操作
>>> numarr2 = numarr * 2
>>> print numarr2
[[ 2  4  6]
 [ 8  0 12]
 [14 16 18]]
>>> print numarr2 + numarr
[[ 3  6  9]
 [12  0 18]
 [21 24 27]]

改变数组的形状:

清单 4. 改变形状
>>> numarr2.shape = (9,)
>>> print numarr2
[ 2  4  6  8  0 12 14 16 18]

Numeric 与 Numarray 之间的区别

总体来看,新的 Numarray 软件包与早期的 Numeric 是 API 兼容的。不过,开发者基于用户经验进行了一些与 Numric 并不兼容的改进。开发者没有破坏任何依赖于 Numeric 的应用程序,而是开创了一个叫做 Numarray 的新项目。在完成本文时,Numarray 还缺少 Numeric 的一些功能,但是已计划实现这些功能。

Numarray 所做的一些改进:

  • 以分层的类结构来组织元素类型,以支持 isinstance() 检验。Numeric 在指定数据类型时只使用字符类型编码 (但是 Numarray 中的初始化软件仍然接受老的字符编码)。
  • 改变了类型强制规则,以保持数组(更为常见)中的类型 ,而不是转换为 Python 标量的类型。
  • 出现了附加的数组属性 (不再只有 getter 和 setter)。
  • 实现了更灵活的异常处理。

新用户不必担心这些变化,就这一点来说,最好一开始就使用 Numarray 而不是 Numeric。

计时的例子

让我们来感受一下在 Numerical Python 中的操作相对于标准 Python 的速度优势。作为一个“演示任务”,我们将创建一个数字序列,然后使它们加倍。首先是标准 Python 方法的一些变体:

清单 5. 对纯 Python 操作的计时
def timer(fun, n, comment=""):
    from time import clock
    start = clock()
    print comment, len(fun(n)), "elements",
    print "in %.2f seconds" % (clock()-start)
def double1(n): return map(lambda n: 2*n, xrange(n))
timer(double1, 5000000, "Running map() on xrange iterator:")
def double2(n): return [2*n for n in xrange(n)]
timer(double2, 5000000, "Running listcomp on xrange iter: ")
def double3(n):
    double = []
    for n in xrange(n):
        double.append(2*n)
    return double
timer(double3, 5000000, "Building new list from iterator: ")

我们可以看出 map() 方法、list comprehension 和传统循环方法之间的速度差别。那么,需要同类元素类型的标准 array 模块呢?它可能会更快一些:

清单 6. 对标准 array 模块的计时
import array
def double4(n): return [2*n for n in array.array('i',range(n))]
timer(double4, 5000000, "Running listcomp on array.array: ")

最后我们来看 Numarray 的速度如何。作为额外对照,我们来看如果必须要将数组还原为一个标准的列表时,它是否同样具有优势:

清单 7. 对 Numarray 操作的计时
from numarray import *
def double5(n): return 2*arange(n)
timer(double5, 5000000, "Numarray scalar multiplication:  ")
def double6(n): return (2*arange(n)).tolist()
timer(double6, 5000000, "Numarray mult, returning list:   ")

现在运行它:

清单 8. 比较结果
$ python2.3 timing.py
Running map() on xrange iterator: 5000000 elements in 13.61 seconds
Running listcomp on xrange iter:  5000000 elements in 16.46 seconds
Building new list from iterator:  5000000 elements in 20.13 seconds
Running listcomp on array.array:  5000000 elements in 25.58 seconds
Numarray scalar multiplication:   5000000 elements in 0.61 seconds
Numarray mult, returning list:    5000000 elements in 3.70 seconds

处理列表的不同技术之间的速度差异不大,也许还是值得注意,因为这是尝试标准的 array 模块时的方法问题。但是 Numarray 一般用不到 1/20 的时间内就可以完成操作。将数组还原为标准列表损失了很大的速度优势。

不应通过这样一个简单的比较就得出结论,但是这种加速可能是典型的。对大规模科学计算来说,将计算的时间由几个月下降到几天或者从几天下降到几个小时,是非常有价值的。

系统建模

Numerical Python 的典型用例是科学建模,或者可能是相关领域,比如图形处理和旋转,或者信号处理。我将通过一个比较实际的问题来说明 Numarray 的许多功能。假设您有一个参量可变的三维物理空间。抽象地说,任何参数化空间,不论有多少维,Numarray 都适用。实际上很容易想像,比如一个房间,它的各个点的温度是不同的。我在 New England 的家已经到了冬天,因而这个问题似乎更有现实意义。

为简单起见,下面我给出的例子中使用的是较小的数组(虽然这可能是显然的,但是还是有必要明确地指出来)。不过,即使是处理有上百万个元素而不仅仅是几十个元素的数组,Numarray 也还是很快;前者可能在真正的科学模型中更为常见。

首先,我们来创建一个“房间”。有很多方法可以完成这项任务,但是最常用的还是使用可调用的 array() 方法。使用这个方法,我们可以生成具有多种初始化参数 (包括来自任何 Python 序列的初始数据) 的 Numerical 数组。不过对于我们的房间来说,用 zeros() 函数就可以生成一个温度均匀的寒冷房间:

清单 9. 初始化房间的温度
>>> from numarray import *
>>> room = zeros((4,3,5),Float)
>>> print room
[[[ 0. 0. 0. 0. 0.]
  [ 0. 0. 0. 0. 0.]
  [ 0. 0. 0. 0. 0.]]
 [[ 0. 0. 0. 0. 0.]
  [ 0. 0. 0. 0. 0.]
  [ 0. 0. 0. 0. 0.]]
 [[ 0. 0. 0. 0. 0.]
  [ 0. 0. 0. 0. 0.]
  [ 0. 0. 0. 0. 0.]]
 [[ 0. 0. 0. 0. 0.]
  [ 0. 0. 0. 0. 0.]
  [ 0. 0. 0. 0. 0.]]]

自上而下每一个二维的“矩阵”代表三维房间的一个水平层面。

首先,我们将整个房间的温度提高到比较舒适的 70 华氏度 (大约是 20 摄氏度):

清单 10. 打开加热器
>>> room += 70
>>> print room
[[[ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]]
 [[ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]]
 [[ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]]
 [[ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]]]

请注意,在我们接下来对 Numarray 数组和 Python 列表进行操作时有很重要的区别。当您选取数组的层面时 -- 我们将会看到,多维数组中的分层方法非常灵活且强大 -- 您得到的不是一个拷贝而是一个“视图”。指向相同的数据可以有多种途径。

让我们具体来看。假设我们房间有一个通风装置,会将地面的温度降低四度:

清单 11. 温度的变化
>>> floor = room[3]
>>> floor -= 4
>>> print room
[[[ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]]
 [[ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]]
 [[ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]]
 [[ 66. 66. 66. 66. 66.]
  [ 66. 66. 66. 66. 66.]
  [ 66. 66. 66. 66. 66.]]]

与此相对,北面墙上的壁炉将每个邻近位置的温度升高了 8 度,而它所在位置的温度为 90 度。

清单 12. 使用壁炉取暖
>>> north = room[:,0]
>>> near_fireplace = north[2:4,2:5]
>>> near_fireplace += 8
>>> north[3,2] = 90  # the fireplace cell itself
>>> print room
[[[ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]]
 [[ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]]
 [[ 70. 78. 78. 78. 70.]
  [ 70. 70. 70. 70. 70.]
  [ 70. 70. 70. 70. 70.]]
 [[ 66. 74. 90. 74. 66.]
  [ 66. 66. 66. 66. 66.]
  [ 66. 66. 66. 66. 66.]]]

这里我们使用了一些比较巧妙的索引方法,可以沿多维的方向来指定层面。这些视图应该保留,以后还会用到。例如,您可能希望知道整个北面墙上的当前温度:

清单 13. 查看北面的墙
>>> print north
[[ 70. 70. 70. 70. 70.]
 [ 70. 70. 70. 70. 70.]
 [ 70. 78. 78. 78. 70.]
 [ 66. 74. 90. 74. 66.]]

更多操作

以上介绍的仅仅是 Numarray 中便捷的函数和数组方法/属性中的一小部分。我希望能给您一些初步的认识;Numarray 文档是深入学习的极好参考资料。

既然我们的房间现在各处的温度不再相同,我们可能需要判断全局的状态。例如,当前房间内的平均温度:

清单 14. 查看平均化后的数组
>>> add.reduce(room.flat)/len(room.flat)
70.066666666666663

这里需要解释一下。您可以对数组进行的所有操作都有相对应的 通用函数 (ufunc)。所以,我们在前面的代码中使用的 floor -= 4 ,可以替换为 subtract(floor,4,floor) 。指定 subtract() 的三个参数,操作就可以正确完成。您还可以用 floor=subtract(floor,4) 来创建 floor 的一个拷贝,但这可能不是您所期望的,因为变化将发生在一个新的数组中,而不是 room 的一个子集中。

然而,unfunc 不仅仅是函数。它们还可以是可调用的对象,具有自己的方法:其中 .reduce() 可能是最为有用的一个。 reduce() 的工作方式如同 Python 中的内置函数 reduce() ,每个操作都是基本的 ufunc (不过这些方法在应用于 Numerical 数组时会快得多)。换句话说, add.reduce() 表示的是 sum()multiply.reduce() 表示的是 product() (这些快捷名称也是定义好了的)。

在求房间各单元温度的和之前,您需要先得到数据的一个一维视图。不然,您得到的是第一维的和,并生成一个降低了维数的新数组。例如:

清单 15. 非平面数组的错误结果
>>> add.reduce(room)
array([[ 276.,  292.,  308.,  292.,  276.],
       [ 276.,  276.,  276.,  276.,  276.],
       [ 276.,  276.,  276.,  276.,  276.]])

这样一个空间和可能会有用,但它并不是我们这里想要得到的。

既然我们是在对一个物理系统建模,我们来让它更真实一些。房间内有微小的气流,使得温度发生变化。在建模时我们可以假设每一个小的时间段内,每个单元会根据它周围的温度进行调整:

清单 16. 微气流模拟
>>> def equalize(room):
...    z,y,x = map(randint, (1,1,1), room.shape)
...    zmin,ymin,xmin = maximum([z-2,y-2,x-2],[0,0,0]).tolist()
...    zmax,ymax,xmax = [z+1,y+1,x+1]
...    region = room[zmin:zmax,ymin:ymax,xmin:xmax].copy()
...    room[z-1,y-1,x-1] = sum(region.flat)/len(region.flat)
...    return room

这个模型当然有一些不现实:单元不会只根据它周围的温度进行调整而不去影响它相邻的单元。尽管如此,还是让我们来看一下它执行的情况。首先我们选择一个随机的单元 -- 或者实际上我们选取的是单元本身在每一维度上的索引值加上 1,因为我们通过 .shape 调用得到的是长度而不是最大的索引值。 zminyminxmin 确保了我们的最小值索引值为 0,不会取到负数; zmaxymaxxmax 实际上并不需要,因为数组每一维的大小减去 1 之后的索引值就被当作最大值来使用(如同 Python 中的列表)。

然后,我们需要定义邻近单元的区域。由于我们的房间很小,所以经常会选择到房间的表面、边沿或者一角 -- 单元的 region 可能会比最大的 27 元素 (3x3x3) 子集要小。这没关系;我们只需要使用正确的分母来计算平均值。这个新的平均温度值被赋给前面随机选择的单元。

您可以在您的模型中执行任意多次的平均化过程。每一次调用只调整一个单元。多次调用会使用房间的某些部分的温度逐渐趋于平均。即使数组是动态改变的, equalize() 函数照样可以返回它的数组。当您只想平均化模型的一个 拷贝时这将非常有用:

清单 17. 执行 equalize()
>>> print equalize(room.copy())
[[[ 70.       70.       70.       70.       70.     ]
  [ 70.       70.       70.       70.       70.     ]
  [ 70.       70.       70.       70.       70.     ]]
 [[ 70.       70.       71.333333  70.       70.     ]
  [ 70.       70.       70.       70.       70.     ]
  [ 70.       70.       70.       70.       70.     ]]
 [[ 70.       78.       78.       78.       70.     ]
  [ 70.       70.       70.       70.       70.     ]
  [ 70.       70.       70.       70.       70.     ]]
 [[ 66.       74.       90.       74.       66.     ]
  [ 66.       66.       66.       66.       66.     ]
  [ 66.       66.       66.       68.       66.     ]]]

结束语

本文仅介绍了 Numarray 的部分功能。它的功能远不止这些。例如,您可以使用填充函数来填充数组,这对于物理模型来说非常有用。您不但可以通过层面而且可以通过索引数组来指定数组的子集 -- 这使您不但可以对数组中不连续的片断进行操作,而且可以 -- 通过 take() 函数 -- 以各种方式重新定义数组的维数和形状。

前面我所描述的大部分操作都是针对于标量和数组的;您还可以执行数组之间的操作,包括那些不同维度的数组之间。这涉及到的内容很多,但通过 API 可以直观地完成所有这些操作。

我鼓励您在自己的系统上安装 Numarray 和 / 或 Numeric。它不难上手,并且它提供的对数组的快速操作可以应用于极广泛的领域 -- 往往是您开始时意想不到的。


相关主题


评论

添加或订阅评论,请先登录注册

static.content.url=http://www.ibm.com/developerworks/js/artrating/
SITE_ID=10
Zone=Linux
ArticleID=22202
ArticleTitle=可爱的 Python: Numerical Python
publish-date=12012003