NumPy-快速处理数据
2.2 ufunc运算
ufunc是universal function的缩写，它是一种能对数组的每个元素进行操作的函数。NumPy内置的许 多ufunc函数都是在C语言级别实现的，因此它们的计算速度非常快。让我们来看一个例子：
>>> x = np.linspace(0, 2*np.pi, 10)
# 对数组x中的每个元素进行正弦计算，返回一个同样大小的新数组
>>> y = np.sin(x)
>>> y
array([ 0.00000000e+00, 6.42787610e-01, 9.84807753e-01, 8.66025404e-01, 3.42020143e-01, -3.42020143e-01, -8.66025404e-01, -9.84807753e-01, -6.42787610e-01, -2.44921271e-16])
先用linspace产生一个从0到2*PI的等距离的10个数，然后将其传递给sin函数，由于np.sin是一个 ufunc函数，因此它对x中的每个元素求正弦值，然后将结果返回，并且赋值给y。计算之后x中的值并
2.2. ufunc运算 25
没有改变，而是新创建了一个数组保存结果。如果我们希望将sin函数所计算的结果直接覆盖到数组x上 去的话，可以将要被覆盖的数组作为第二个参数传递给ufunc函数。例如：:
>>> t = np.sin(x,x)
>>> x
array([ 0.00000000e+00, 6.42787610e-01, 9.84807753e-01, 8.66025404e-01, 3.42020143e-01, -3.42020143e-01, -8.66025404e-01, -9.84807753e-01, -6.42787610e-01, -2.44921271e-16])
>>> id(t) == id(x)
True
import numpy as np
x = [i * 0.001 for i in xrange(1000000)]
start = time.clock() for i, t in enumerate(x):
x[i] = math.sin(t)
print "math.sin:", time.clock() - start
x = [i * 0.001 for i in xrange(1000000)]
x = np.array(x) start = time.clock() np.sin(x,x)
print "numpy.sin:", time.clock() - start
# 输出
# math.sin: 1.15426932753
# numpy.sin: 0.0882399858083
在我的电脑上计算100万次正弦值，numpy.sin比math.sin快10倍多。这得利于numpy.sin在C语言级 别的循环计算。numpy.sin同样也支持对单个数值求正弦，例如：numpy.sin(0.5)。不过值得注意的 是，对单个数的计算math.sin则比numpy.sin快得多了，让我们看下面这个测试程序：
x = [i * 0.001 for i in xrange(1000000)]
start = time.clock()
26 第 2 章 NumPy-快速处理数据
用Python做科学计算
for i, t in enumerate(x):
x[i] = np.sin(t)
print "numpy.sin loop:", time.clock() - start
# 输出
# numpy.sin loop: 5.72166965355
请注意numpy.sin的计算速度只有math.sin的1/5。这是因为numpy.sin为了同时支持数组和单个值的 短，因此在导入时不建议使用*号全部载入，而是应该使用import numpy as np的方式载入，这样我 们可以根据需要选择合适的函数调用。
NumPy中有众多的ufunc函数为我们提供各式各样的计算。除了sin这种单输入函数之外，还有许多多 个输入的函数，add函数就是一个最常用的例子。先来看一个例子:
>>> a = np.arange(0,4)
>>> a
array([0, 1, 2, 3])
>>> b = np.arange(1,5)
>>> b
2.2. ufunc运算 27
y = x1 + x2: add(x1, x2 [, y]) y = x1 - x2: subtract(x1, x2 [, y]) y = x1 * x2: multiply (x1, x2 [, y])
y = x1 / x2: divide (x1, x2 [, y]), 如果两个数组的元素为整数，那么用整数除法 y = x1 / x2: true divide (x1, x2 [, y]), 总是返回精确的商
y = x1 // x2: floor divide (x1, x2 [, y]), 总是对返回值取整 y = -x: negative(x [,y])
y = x1**x2: power(x1, x2 [, y])
y = x1 % x2: remainder(x1, x2 [, y]), mod(x1, x2, [, y])
数组对象支持这些操作符，极大地简化了算式的编写，不过要注意如果你的算式很复杂，并且要运算 的数组很大的话，会因为产生大量的中间结果而降低程序的运算效率。例如：假设a b c三个数组采用 算式x=a*b+c计算，那么它相当于:
t = a * b x = t + c del t
也就是说需要产生一个数组t保存乘法的计算结果，然后再产生最后的结果数组x。我们可以通过手工将 一个算式分解为x = a*b; x += c，以减少一次内存分配。
通过组合标准的ufunc函数的调用，可以实现各种算式的数组计算。不过有些时候这种算式不易编写，
而针对每个元素的计算函数却很容易用Python实现，这时可以用frompyfunc函数将一个计算单个元 素的函数转换成ufunc函数。这样就可以方便地用所产生的ufunc函数对数组进行计算了。让我们来看 一个例子。
我们想用一个分段函数描述三角波，三角波的样子如下：
图 2.4 - 三角波可以用分段函数进行计算
我们很容易根据上图所示写出如下的计算三角波某点y坐标的函数:
28 第 2 章 NumPy-快速处理数据
用Python做科学计算
def triangle_wave(x, c, c0, hc):
x = x - int(x) # 三角波的周期为1，因此只取x坐标的小数部分进行计算
x = np.linspace(0, 2, 1000)
y = np.array([triangle_wave(t, 0.6, 0.4, 1.0) for t in x])
这种做法每次都需要使用列表包容语法调用函数，对于多维数组是很麻烦的。让我们来看看如何用 frompyfunc函数来解决这个问题：
triangle_ufunc = np.frompyfunc( lambda x: triangle_wave(x, 0.6, 0.4, 1.0), 1, 1) y2 = triangle_ufunc(x)
frompyfunc的调用格式为frompyfunc(func, nin, nout)，其中func是计算单个元素的函数，nin是此 函数的输入参数的个数，nout是此函数的返回值的个数。虽然triangle_wave函数有4个参数，但是由 于后三个c, c0, hc在整个计算中值都是固定的，因此所产生的ufunc函数其实只有一个参数。为了满足 这个条件，我们用一个lambda函数对triangle_wave的参数进行一次包装。这样传入frompyfunc的函 数就只有一个参数了。这样子做，效率并不是太高，另外还有一种方法：
def triangle_func(c, c0, hc):
def trifunc(x):
# 用trifunc函数创建一个ufunc函数，可以直接对数组进行计算, 不过通过此函数
# 计算得到的是一个Object数组，需要进行类型转换 return np.frompyfunc(trifunc, 1, 1)
y2 = triangle_func(0.6, 0.4, 1.0)(x)
我们通过函数triangle_func包装三角波的三个参数，在其内部定义一个计算三角波的函数trifunc，
trifunc函数在调用时会采用triangle_func的参数进行计算。最后triangle_func返回用frompyfunc转 换结果。
2.2. ufunc运算 29
值得注意的是用frompyfunc得到的函数计算出的数组元素的类型为object，因为frompyfunc函数无
>>> a = np.arange(0, 60, 10).reshape(-1, 1)
>>> a
array([[ 0], [10], [20], [30], [40], [50]])
>>> a.shape
(6, 1)
再创建一维数组b，其shape为(5,)：
>>> b = np.arange(0, 5)
>>> b
用Python做科学计算
>>> b = b.repeat(6,axis=0)
>>> b
>>> a = a.repeat(5, axis=1)
>>> a
2.2. ufunc运算 31
>>> x,y = np.ogrid[0:5,0:5]
>>> x, y = np.ogrid[0:1:4j, 0:1:3j]
>>> x
import numpy as np
from enthought.mayavi import mlab
x, y = np.ogrid[-2:2:20j, -2:2:20j]
z = x * np.exp( - x**2 - y**2)
pl = mlab.surf(x, y, z, warp_scale="auto")
32 第 2 章 NumPy-快速处理数据
用Python做科学计算
mlab.axes(xlabel='x', ylabel='y', zlabel='z') mlab.outline(pl)
此程序使用mayavi的mlab库快速绘制3D曲面，关于mlab的相关内容将在今后的章节进行介绍。
图 2.5 -使用Mayavi将二维数组绘制成3D曲面
2.2.2 ufunc的方法
ufunc函数本身还有些方法，这些方法只对两个输入一个输出的ufunc函数有效，其它的ufunc对象调 用这些方法时会抛出ValueError异常。
reduce 方法和Python的reduce函数类似，它沿着axis轴对array进行操作，相当于将<op>运算符插 入到沿axis轴的所有子数组或者元素当中。
<op>.reduce (array=, axis=0, dtype=None)
例如：
>>> np.add.reduce([1,2,3]) # 1 + 2 + 3
6
>>> np.add.reduce([[1,2,3],[4,5,6]], axis=1) # 1,4 + 2,5 + 3,6
array([ 6, 15])
accumulate 方法和reduce方法类似，只是它返回的数组和输入的数组的shape相同，保存所有的中间 计算结果：
>>> np.add.accumulate([1,2,3])
array([1, 3, 6])
2.2. ufunc运算 33
>>> np.add.accumulate([[1,2,3],[4,5,6]], axis=1) array([[ 1, 3, 6],
[ 4, 9, 15]])
reduceat 方法计算多组reduce的结果，通过indices参数指定一系列reduce的起始和终了位置。
reduceat的计算有些特别，让我们通过一个例子来解释一下：
>>> a = np.array([1,2,3,4])
>>> result = np.add.reduceat(a,indices=[0,1,0,2,0,3,0])
>>> result
array([ 1, 2, 3, 3, 6, 4, 10])
对于indices中的每个元素都会调用reduce函数计算出一个值来，因此最终计算结果的长度和indices 的长度相同。结果result数组中除最后一个元素之外，都按照如下计算得出：
if indices[i] < indices[i+1]:
result[i] = np.reduce(a[indices[i]:indices[i+1]]) else:
result[i] = a[indices[i]
而最后一个元素如下计算:
>>> a.shape += (1,)*b.ndim
>>> <op>(a,b)
>>> a = a.squeeze()
其中squeeze的功能是剔除数组a中长度为1的轴。如果你看不太明白这个等同程序的话，让我们来看 一个例子：
34 第 2 章 NumPy-快速处理数据
用Python做科学计算
>>> a = np.matrix([[1,2,3],[5,5,6],[7,9,9]])
>>> a*a**-1
matrix([[ 1.00000000e+00, 1.66533454e-16, -8.32667268e-17], [ -2.77555756e-16, 1.00000000e+00, -2.77555756e-17], [ 1.66533454e-16, 5.55111512e-17, 1.00000000e+00]])
因为a是用matrix创建的矩阵对象，因此乘法和幂运算符都变成了矩阵运算，于是上面计算的是矩 阵a和其逆矩阵的乘积，结果是一个单位矩阵。
矩阵的乘积可以使用dot函数进行计算。对于二维数组，它计算的是矩阵乘积，对于一维数组，它计算 的是其点积。当需要将一维数组当作列矢量或者行矢量进行矩阵运算时，推荐先使用reshape函数将一