Python中的Numeric包和Numarray包使用教程
要了解NumericalPython软件包的第一件事情是，NumericalPython不会让您去做标准Python不能完成的任何工作。它只是让您以快得多的速度去完成标准Python能够完成的相同任务。实际上不仅仅如此；许多数组操作用Numeric或者Numarray来表达比起用标准Python数据类型和语法来表达要优雅得多。不过，惊人的速度才是吸引用户使用NumericalPython的主要原因。
其实，NumericalPython只是实现了一个新的数据类型：数组。与可以包含不同类型元素的列表、元组和词典不同的是，Numarray数组只能包含同一类型的数据。Numarray数组的另一个优点是，它可以是多维的--但是数组的维度与列表的简单嵌套稍有不同。NumericalPython借鉴了程序员的实践经验（尤其是那些有科学计算背景的程序员，他们抽象出了APL、FORTRAN、MATLAB和S等语言中数组的最佳功能），创建了可以灵活改变形状和维度的数组。我们很快会回来继续这一话题。
在NumericalPython中对数组的操作是按元素进行的。虽然二维数组与线性代数中的矩阵类似，但是对它们的操作(比如乘)与线性代数中的操作(比如矩阵乘)是完全不同的。
让我们来看一个关于上述问题的的具体例子。在纯Python中，您可以这样创建一个“二维列表”：
清单1.Python的嵌套数组
>>>pyarr=[[1,2,3], ...[4,5,6], ...[7,8,9]] >>>printpyarr [[1,2,3],[4,5,6],[7,8,9]] >>>pyarr[1][1]=0 >>>printpyarr [[1,2,3],[4,0,6],[7,8,9]]
很好，但是您对这种结构所能做的只是通过单独的(或者多维的)索引来设置和检索元素。与此相比，Numarray数组要更灵活：
清单2.NumericalPython数组
>>>fromnumarrayimport* >>>numarr=array(pyarr) >>>printnumarr [[123] [406] [789]]
改变并不大，但是使用Numarray进行的操作如何呢?下面是一个例子：
清单3.元素操作
>>>numarr2=numarr*2 >>>printnumarr2 [[246] [8012] [141618]] >>>printnumarr2+numarr [[369] [12018] [212427]]
改变数组的形状：
清单4.改变形状
>>>numarr2.shape=(9,) >>>printnumarr2 [2468012141618]
Numeric与Numarray之间的区别
总体来看，新的Numarray软件包与早期的Numeric是API兼容的。不过，开发者基于用户经验进行了一些与Numric并不兼容的改进。开发者没有破坏任何依赖于Numeric的应用程序，而是开创了一个叫做Numarray的新项目。在完成本文时，Numarray还缺少Numeric的一些功能，但是已计划实现这些功能。
Numarray所做的一些改进：
- 以分层的类结构来组织元素类型，以支持isinstance()检验。Numeric在指定数据类型时只使用字符类型编码(但是Numarray中的初始化软件仍然接受老的字符编码)。
- 改变了类型强制规则，以保持数组(更为常见)中的类型，而不是转换为Python标量的类型。
- 出现了附加的数组属性(不再只有getter和setter)。
- 实现了更灵活的异常处理。
新用户不必担心这些变化，就这一点来说，最好一开始就使用Numarray而不是Numeric。
计时的例子
让我们来感受一下在NumericalPython中的操作相对于标准Python的速度优势。作为一个“演示任务”，我们将创建一个数字序列，然后使它们加倍。首先是标准Python方法的一些变体：
清单5.对纯Python操作的计时
deftimer(fun,n,comment=""): fromtimeimportclock start=clock() printcomment,len(fun(n)),"elements", print"in%.2fseconds"%(clock()-start) defdouble1(n):returnmap(lambdan:2*n,xrange(n)) timer(double1,5000000,"Runningmap()onxrangeiterator:") defdouble2(n):return[2*nforninxrange(n)] timer(double2,5000000,"Runninglistcomponxrangeiter:") defdouble3(n): double=[] forninxrange(n): double.append(2*n) returndouble timer(double3,5000000,"Buildingnewlistfromiterator:")
我们可以看出map()方法、listcomprehension和传统循环方法之间的速度差别。那么，需要同类元素类型的标准array模块呢？它可能会更快一些：
清单6.对标准array模块的计时
importarray
defdouble4(n):return[2*nforninarray.array('i',range(n))]
timer(double4,5000000,"Runninglistcomponarray.array:")
最后我们来看Numarray的速度如何。作为额外对照，我们来看如果必须要将数组还原为一个标准的列表时，它是否同样具有优势：
清单7.对Numarray操作的计时
fromnumarrayimport* defdouble5(n):return2*arange(n) timer(double5,5000000,"Numarrayscalarmultiplication:") defdouble6(n):return(2*arange(n)).tolist() timer(double6,5000000,"Numarraymult,returninglist:")
现在运行它：
清单8.比较结果
$python2.3timing.py Runningmap()onxrangeiterator:5000000elementsin13.61seconds Runninglistcomponxrangeiter:5000000elementsin16.46seconds Buildingnewlistfromiterator:5000000elementsin20.13seconds Runninglistcomponarray.array:5000000elementsin25.58seconds Numarrayscalarmultiplication:5000000elementsin0.61seconds Numarraymult,returninglist:5000000elementsin3.70seconds
处理列表的不同技术之间的速度差异不大，也许还是值得注意，因为这是尝试标准的array模块时的方法问题。但是Numarray一般用不到1/20的时间内就可以完成操作。将数组还原为标准列表损失了很大的速度优势。
不应通过这样一个简单的比较就得出结论，但是这种加速可能是典型的。对大规模科学计算来说，将计算的时间由几个月下降到几天或者从几天下降到几个小时，是非常有价值的。
系统建模
NumericalPython的典型用例是科学建模，或者可能是相关领域，比如图形处理和旋转，或者信号处理。我将通过一个比较实际的问题来说明Numarray的许多功能。假设您有一个参量可变的三维物理空间。抽象地说，任何参数化空间，不论有多少维，Numarray都适用。实际上很容易想像，比如一个房间，它的各个点的温度是不同的。我在NewEngland的家已经到了冬天，因而这个问题似乎更有现实意义。
为简单起见，下面我给出的例子中使用的是较小的数组(虽然这可能是显然的，但是还是有必要明确地指出来)。不过，即使是处理有上百万个元素而不仅仅是几十个元素的数组，Numarray也还是很快；前者可能在真正的科学模型中更为常见。
首先，我们来创建一个“房间”。有很多方法可以完成这项任务，但是最常用的还是使用可调用的array()方法。使用这个方法，我们可以生成具有多种初始化参数(包括来自任何Python序列的初始数据)的Numerical数组。不过对于我们的房间来说，用zeros()函数就可以生成一个温度均匀的寒冷房间：
清单9.初始化房间的温度
>>>fromnumarrayimport* >>>room=zeros((4,3,5),Float) >>>printroom [[[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 >>>printroom [[[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 >>>printroom [[[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#thefireplacecellitself >>>printroom [[[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.查看北面的墙
>>>printnorth [[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.微气流模拟
>>>defequalize(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) ...returnroom
这个模型当然有一些不现实：单元不会只根据它周围的温度进行调整而不去影响它相邻的单元。尽管如此，还是让我们来看一下它执行的情况。首先我们选择一个随机的单元--或者实际上我们选取的是单元本身在每一维度上的索引值加上1，因为我们通过.shape调用得到的是长度而不是最大的索引值。zmin、ymin和xmin确保了我们的最小值索引值为0，不会取到负数；zmax、ymax和xmax实际上并不需要，因为数组每一维的大小减去1之后的索引值就被当作最大值来使用(如同Python中的列表)。
然后，我们需要定义邻近单元的区域。由于我们的房间很小，所以经常会选择到房间的表面、边沿或者一角--单元的region可能会比最大的27元素(3x3x3)子集要小。这没关系；我们只需要使用正确的分母来计算平均值。这个新的平均温度值被赋给前面随机选择的单元。
您可以在您的模型中执行任意多次的平均化过程。每一次调用只调整一个单元。多次调用会使用房间的某些部分的温度逐渐趋于平均。即使数组是动态改变的，equalize()函数照样可以返回它的数组。当您只想平均化模型的一个拷贝时这将非常有用：
清单17.执行equalize()
>>>printequalize(room.copy()) [[[70.70.70.70.70.] [70.70.70.70.70.] [70.70.70.70.70.]] [[70.70.71.33333370.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。它不难上手，并且它提供的对数组的快速操作可以应用于极广泛的领域--往往是您开始时意想不到的。