可爱的 Python: Numerical Python
linxh
|
1#
linxh 发表于 2007-01-25 01:43
可爱的 Python: Numerical Python使用 Numeric 软件包和 Numarray 软件包 作者： David Mertz, Ph.D. , 开发人员, Gnosis Software, Inc. Numerical Python (通常称为 NumPy) 是一个广为应用的 Python 扩展库，用于快速处理任意维数的固定类型数组。由于底层代码是充分优化的 C 语言代码，因而对数组的主要操作在 NumPy 调用中执行时，速度不再受到 Python 解释器的限制。因为 NumPy 已经取得了这样的成功，所以 NumPy 的开发者将用一个叫做 Numarray 的新模块来取代 NumPy，新模块基本上 (但并不是完全) 与 NumPy 兼容。在本文中，David 介绍了 NumPy 的一般功能，以及 Numarray 将要带来的一些特殊改进。 要了解 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 所做的一些改进：
新用户不必担心这些变化，就这一点来说，最好一开始就使用 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 调用得到的是长度而不是最大的索引值。 zmin 、 ymin 和 xmin 确保了我们的最小值索引值为 0，不会取到负数； zmax 、 ymax 和 xmax 实际上并不需要，因为数组每一维的大小减去 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。它不难上手，并且它提供的对数组的快速操作可以应用于极广泛的领域 -- 往往是您开始时意想不到的。 参考资料
|