ufunc是universal function的缩写，它是一种能对数组的每个元素进行运算的函数。NumPy内置的许多ufunc函数都是用C语言实现的，因此它们的计算速度非常快。让我们先看一个例子：

x = np.linspace(0, 2*np.pi, 10)
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.44929360e-16])

先用linspace()产生一个从0到2π的等差数组，然后将其传递给np.sin()函数计算每个元素的正弦值。由于np.sin()是一个ufunc函数，因此在其内部对数组x的每个元素进行循环，分别计算它们的正弦值，并返回一个保存各个计算结果的数组。运算之后数组x中的值并没有改变，而是新创建了一个数组来保存结果。也可以通过out参数指定保存计算结果的数组。因此如果希望直接在数组x中保存结果，可以将它传递给out参数：

t = np.sin(x, out=x)
t is x
True

ufunc函数的返回值仍然是计算的结果，只不过它就是数组x。下面比较np.sin()和Python标准库的math.sin()的计算速度：

import math

x = [i * 0.001 for i in xrange(1000000)]

def sin_math(x):
for i, t in enumerate(x):
x[i] = math.sin(t)

def sin_numpy(x):
np.sin(x, x)

def sin_numpy_loop(x):
for i, t in enumerate(x):
x[i] = np.sin(t)

xl = x[:]
%time sin_math(x)

xa = np.array(x)
%time sin_numpy(xa)

xl = x[:]
%time sin_numpy_loop(x)
Wall time: 302 ms
Wall time: 30 ms
Wall time: 1.28 s

可以看出，np.sin()比math.sin()快10倍多，这得益于np.sin()在C语言级别的循环计算。

列表推导式比循环更快

事实上，标准Python中有比for循环更快的方案：使用列表推导式x = [math.sin(t) for t in x]。但是列表推导式将产生一个新的列表，而不是直接修改原列表中的元素。

np.sin()同样也支持计算单个数值的正弦值。不过值得注意的是，对单个数值的计算，math.sin()则比np.sin()快很多。在Python级别进行循环时，np.sin()的计算速度只有math.sin()的1/6。这是因为：np.sin()为了同时支持数组和单个数值的计算，其C语言的内部实现要比math.sin()复杂很多。此外，对于单个数值的计算，np.sin()的返回值类型和math.sin()的不同，math.sin()返回的是Python的标准float类型，而np.sin()返回float64类型：

type(math.sin(0.5)) type(np.sin(0.5))
------------------- -----------------
float numpy.float64

通过下标运算获取的数组元素的类型为NumPy中定义的类型，将其转换为Python的标准类型还需要花费额外的时间。为了解决这个问题，数组提供了item()方法，它用来获取数组中的单个元素，并且直接返回标准的Python数值类型：

a = np.arange(6.0).reshape(2, 3)
print a.item(1, 2), type(a.item(1, 2)), type(a[1, 2])
5.0 

通过上面的例子我们了解了如何最有效率地使用math模块和NumPy中的数学函数。由于它们各有优缺点，因此在导入时不建议使用import *全部载入，而是应该使用import numpy as np载入，这样可以根据需要选择合适的函数。

一、四则运算

NumPy提供了许多ufunc函数，例如计算两个数组之和的add()函数：

a = np.arange(0, 4)
b = np.arange(1, 5)
np.add(a, b)
array([1, 3, 5, 7])

add()返回一个数组，它的每个元素都是两个参数数组的对应元素之和。如果没有指定out参数，那么它将创建一个新的数组来保存计算结果。如果指定了第三个参数out，则不产生新的数组，而直接将结果保存进指定的数组。

np.add(a, b, a)
a
array([1, 3, 5, 7])

NumPy为数组定义了各种数学运算操作符，因此计算两个数组相加可以简单地写为a + b，而np.add(a, b, a)则可以用a += b来表示。表1列出了数组的运算符以及与之对应的ufunc函数，注意除号的意义根据是否激活__future__.division有所不同。

表1　数组的运算符以及对应的ufunc函数

数组对象支持操作符，极大地简化了算式的编写，不过要注意如果算式很复杂，并且要运算的数组很大，将会因为产生大量的中间结果而降低程序的运算速度。例如，假设对a、b、c三个数组采用算式x=a*b+c加以计算，那么它相当于：

t = a * b
x = t + c
del t

也就是说，需要产生一个临时数组t来保存乘法的运算结果，然后再产生最后的结果数组x。可以将算式分解为下面的两行语句，以减少一次内存分配：

x = a*b
x += c

二、比较运算和布尔运算

使用==、>等比较运算符对两个数组进行比较，将返回一个布尔数组，它的每个元素值都是两个数组对应元素的比较结果。例如：

np.array([1, 2, 3]) < np.array([3, 2, 1])
array([ True, False, False], dtype=bool)

每个比较运算符也与一个ufunc函数对应，表2是比较运算符与ufunc函数的对照表。

表2　比较运算符与相应的ufunc函数

由于Python中的布尔运算使用and、or和not等关键字，它们无法被重载，因此数组的布尔运算只能通过相应的ufunc函数进行。这些函数名都以logical_开头，在IPython中使用自动补全可以很容易地找到它们：

>>> np.logical # 按Tab键进行自动补全
np.logical_and np.logical_not np.logical_or np.logical_xor

下面是一个使用logical_or()进行“或运算”的例子：

a = np.arange(5)
b = np.arange(4, -1, -1)
print a == b
print a > b
print np.logical_or(a == b, a > b) # 和 a>=b 相同
[False False True False False]
[False False False True True]
[False False True True True]

对两个布尔数组使用and、or和not等进行布尔运算，将抛出ValueError异常。因为布尔数组中有True也有False，所以NumPy无法确定用户的运算目的：

a == b and a > b
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
 in ()
----> 1 a == b and a > b

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

错误信息告诉我们可以使用数组的any()或all()方法，在NumPy中同时也定义了any()和all()函数，它们的用法和Python内置的any()和all()类似。只要数组中有一个元素值为True，any()就返回True；而只有当数组的全部元素都为True时，all()才返回True。

np.any(a == b) np.any(a == b) and np.any(a > b)
-------------- --------------------------------
True True

以bitwise_开头的函数是位运算函数，包括bitwise_and、bitwise_not、bitwise_or和bitwise_xor等。也可以使用&、~、|和^等操作符进行计算。

对于布尔数组来说，位运算和布尔运算的结果相同。但在使用时要注意，位运算符的优先级比比较运算符高，因此需要使用括号提高比较运算符的运算优先级。例如：

(a == b) | (a > b)
array([False, False, True, True, True], dtype=bool)

整数数组的位运算和C语言的位运算相同，在使用时要注意元素类型的符号，例如下面的arange()所创建的数组的元素类型为32位符号整数，因此对正数按位取反将得到负数。以整数0为例，按位取反的结果是0xFFFFFFFF，在32位符号整数中，这个值表示-1。

~ np.arange(5)
array([-1, -2, -3, -4, -5])

而如果对8位无符号整数数组进行按位取反运算：

~ np.arange(5, dtype=np.uint8)
array([255, 254, 253, 252, 251], dtype=uint8)

同样的整数0，按位取反的结果是0xFF，当它是8位无符号整数时，它的值是255。

三、自定义ufunc函数

通过NumPy提供的标准ufunc函数，可以组合出复杂的表达式，在C语言级别对数组的每个元素进行计算。但有时这种表达式不易编写，而对每个元素进行计算的程序却很容易用Python实现，这时可以用frompyfunc()将计算单个元素的函数转换成ufunc函数，这样就可以方便地用所产生的ufunc函数对数组进行计算了。

例如，我们可以用一个分段函数描述三角波，三角波的形状如图1所示，它分为三段：上升段、下降段和平坦段。

图1　三角波可以用分段函数进行计算

根据图1，我们很容易写出计算三角波上某点的Y坐标的函数。显然triangle_wave()只能计算单个数值，不能对数组直接进行处理。

def triangle_wave(x, c, c0, hc):
x = x - int(x) # 三角波的周期为1，因此只取x坐标的小数部分进行计算
if x >= c: r = 0.0
elif x < c0: r = x / c0 * hc
else: r = (c-x) / (c-c0) * hc
return r

我们可以用下面的程序，先使用列表推导式计算出一个列表，然后用array()将列表转换为数组。这种做法每次都需要使用列表推导式语法调用函数，这对于多维数组很麻烦。

x = np.linspace(0, 2, 1000)
y1 = np.array([triangle_wave(t, 0.6, 0.4, 1.0) for t in x])

通过frompyfunc()可以将计算单个值的函数转换为能对数组的每个元素进行计算的ufunc函数。frompyfunc()的调用格式为：

frompyfunc(func, nin, nout)

其中：func是计算单个元素的函数，nin是func的输入参数的个数，nout是func的返回值的个数。下面的程序使用frompyfunc()将triangle_wave()转换为ufunc函数对象triangle_ufunc1：

triangle_ufunc1 = np.frompyfunc(triangle_wave, 4, 1)
y2 = triangle_ufunc1(x, 0.6, 0.4, 1.0)

值得注意的是，triangle_ufunc1()所返回的数组的元素类型是object，因此还需要调用数组的astype()方法，以将其转换为双精度浮点数组：

y2.dtype y2.astype(np.float).dtype
---------- -------------------------
dtype('O') dtype('float64')

使用vectorize()也可以实现和frompyfunc()类似的功能，但它可以通过otypes参数指定返回的数组的元素类型。otypes参数可以是一个表示元素类型的字符串，也可以是一个类型列表，使用列表可以描述多个返回数组的元素类型。下面的程序使用vectorize()计算三角波：

triangle_ufunc2 = np.vectorize(triangle_wave, otypes=[np.float])
y3 = triangle_ufunc2(x, 0.6, 0.4, 1.0)

最后我们验证一下结果：

np.all(y1 == y2) np.all(y2 == y3)
---------------- ----------------
True True

四、广播

当使用ufunc函数对两个数组进行计算时，ufunc函数会对这两个数组的对应元素进行计算，因此它要求这两个数组的形状相同。如果形状不同，会进行如下广播(broadcasting)处理：

1）让所有输入数组都向其中维数最多的数组看齐，shape属性中不足的部分都通过在前面加1补齐。

2）输出数组的shape属性是输入数组的shape属性的各个轴上的最大值。

3）如果输入数组的某个轴的长度为1或与输出数组的对应轴的长度相同，这个数组能够用来计算，否则出错。

4）当输入数组的某个轴的长度为1时，沿着此轴运算时都用此轴上的第一组值。

上述4条规则理解起来可能比较费劲，下面让我们看一个实际的例子。

先创建一个二维数组a，其形状为(6,1)：

a = np.arange(0, 60, 10).reshape(-1, 1)
a a.shape
------ -------
[[ 0], (6, 1)
[10],
[20],
[30],
[40],
[50]]

再创建一维数组b，其形状为(5,)：

b = np.arange(0, 5)
b b.shape
--------------- -------
[0, 1, 2, 3, 4] (5,)

计算a与b的和，得到一个加法表，它相当于计算两个数组中所有元素对的和，得到一个形状为(6,5)的数组：

c = a + b
c c.shape
---------------------- -------
[[ 0, 1, 2, 3, 4], (6, 5)
[10, 11, 12, 13, 14],
[20, 21, 22, 23, 24],
[30, 31, 32, 33, 34],
[40, 41, 42, 43, 44],
[50, 51, 52, 53, 54]]

由于a和b的维数不同，根据规则1)，需要让b的shape属性向a对齐，于是在b的shape属性前加1，补齐为(1,5)。相当于做了如下计算：

b.shape = 1, 5
b b.shape
----------------- -------
[[0, 1, 2, 3, 4]] (1, 5)

这样，加法运算的两个输入数组的shape属性分别为(6,1)和(1,5)，根据规则2)，输出数组的各个轴的长度为输入数组各个轴的长度的最大值，可知输出数组的shape属性为(6,5)。

由于b的第0轴的长度为1，而a的第0轴的长度为6，为了让它们在第0轴上能够相加，需要将b的第0轴的长度扩展为6，这相当于：

b = b.repeat(6, axis=0)
b b.shape
----------------- -------
[[0, 1, 2, 3, 4], (6, 5)
[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4]]

这里的repeat()方法沿着axis参数指定的轴复制数组中各个元素的值。由于a的第1轴的长度为1，而b的第1轴的长度为5，为了让它们在第1轴上能够相加，需要将a的第1轴的长度扩展为5，这相当于：

a = a.repeat(5, axis=1)
a a.shape
---------------------- -------
[[ 0, 0, 0, 0, 0], (6, 5)
[10, 10, 10, 10, 10],
[20, 20, 20, 20, 20],
[30, 30, 30, 30, 30],
[40, 40, 40, 40, 40],
[50, 50, 50, 50, 50]]

经过上述处理之后，a和b就可以按对应元素进行相加运算了。当然，在执行a + b运算时，NumPy内部并不会真正将长度为1的轴用repeat()进行扩展，这样太浪费内存空间了。由于这种广播计算很常用，因此NumPy提供了ogrid对象，用于创建广播运算用的数组。

x, y = np.ogrid[:5, :5]
x y
----- -----------------
[[0], [[0, 1, 2, 3, 4]]
[1],
[2],
[3],
[4]]

此外，NumPy还提供了mgrid对象，它的用法和ogrid对象类似，但是它所返回的是进行广播之后的数组：

x, y = np.mgrid[:5, :5]
x y
----------------- -----------------
[[0, 0, 0, 0, 0], [[0, 1, 2, 3, 4],
[1, 1, 1, 1, 1], [0, 1, 2, 3, 4],
[2, 2, 2, 2, 2], [0, 1, 2, 3, 4],
[3, 3, 3, 3, 3], [0, 1, 2, 3, 4],
[4, 4, 4, 4, 4]] [0, 1, 2, 3, 4]]

ogrid是一个很有趣的对象，它像多维数组一样，用切片元组作为下标，返回的是一组可以用来广播计算的数组。其切片下标有两种形式：

· 开始值:结束值:步长，和np.arange(开始值,结束值,步长)类似。

· 开始值:结束值:长度j，当第三个参数为虚数时，它表示所返回的数组的长度，和np.linspace(开始值,结束值,长度)类似。

x, y = np.ogrid[:1:4j, :1:3j]
x y
--------------- --------------------
[[ 0. ], [[ 0. , 0.5, 1. ]]
[ 0.33333333],
[ 0.66666667],
[ 1. ]]

利用ogrid的返回值，我们很容易计算二元函数在等间距网格上的值。下面是绘制三维曲面（x，y）=xex2-y2的程序：

x, y = np.ogrid[-2:2:20j, -2:2:20j]
z = x * np.exp( - x**2 - y**2)

图2为使用ogrid计算的三维曲面。

图2　使用ogrid计算二元函数的曲面

为了充分利用ufunc函数的广播功能，我们经常需要调整数组的形状，因此数组支持特殊的下标对象None，它表示在None对应的位置创建一个长度为1的新轴，例如对于一维数组a，a[None, :]和a.reshape(1, -1)等效，而a[:, None]和a.reshape(-1, 1)等效：

a = np.arange(4)
a[None, :] a[:, None]
-------------- ----------
[[0, 1, 2, 3]][[0],
[1],
[2],
[3]]

下面的例子利用None作为下标，实现广播运算：

x = np.array([0, 1, 4, 10])
y = np.array([2, 3, 8])
x[None, :] + y[:, None]
array([[ 2, 3, 6, 12],
[ 3, 4, 7, 13],
[ 8, 9, 12, 18]])

还可以使用ix_()将两个一维数组转换成可广播的二维数组：

gy, gx = np.ix_(y, x)
gx gy gx + gy
------------------ ----- ------------------
[[ 0, 1, 4, 10]] [[2], [[ 2, 3, 6, 12],
[3], [ 3, 4, 7, 13],
[8]] [ 8, 9, 12, 18]]

在上面的例子中，通过ix_()将数组x和y转换成能进行广播运算的二维数组。注意数组y对应广播运算结果中的第0轴，而数组x与第1轴对应。ix_()的参数可以是N个一维数组，它将这些数组转换成N维空间中可广播的N维数组。

五、ufunc的方法

ufunc函数对象本身还有一些方法函数，这些方法只对两个输入、一个输出的ufunc函数有效，其他的ufunc对象调用这些方法时会抛出ValueError异常。

reduce()方法和Python的reduce()函数类似，它沿着axis参数指定的轴对数组进行操作，相当于将运算符插入到沿axis轴的所有元素之间：.reduce(array, axis=0, dtype=None)。例如：

r1 = np.add.reduce([1, 2, 3]) # 1 + 2 + 3
r2 = np.add.reduce([[1, 2, 3], [4, 5, 6]], axis=1) # (1+2+3),(4+5+6)
r1 r2
-- --------
6 [ 6, 15]

accumulate()方法和reduce()类似，只是它返回的数组和输入数组的形状相同，保存所有的中间计算结果：

a1 = np.add.accumulate([1, 2, 3])
a2 = np.add.accumulate([[1, 2, 3], [4, 5, 6]], axis=1)
a1 a2
--------- --------------
[1, 3, 6] [[ 1, 3, 6],
[ 4, 9, 15]]

reduceat()方法计算多组reduce()的结果，通过indices参数指定一系列的起始和终止位置。它的计算有些特别，让我们通过例子详细解释一下：

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参数中的每个元素都会计算出一个值，因此最终的计算结果和indices参数的长度相同。结果数组result中除最后一个元素之外，都按照如下计算得出：

if indices[i] < indices[i+1]:
result[i] = .reduce(a[indices[i]:indices[i+1]])
else:
result[i] = a[indices[i]]

而最后一个元素如下计算:

.reduce(a[indices[-1]:])

因此在上面的例子中，数组result的每个元素按照如下计算得出：

1 : a[0] -> 1
2 : a[1] -> 2
3 : a[0] + a[1] -> 1 + 2
3 : a[2] -> 3
6 : a[0] + a[1] + a[2] -> 1 + 2 + 3 = 6
4 : a[3] -> 4
10: a[0] + a[1] + a[2] + a[4] -> 1 + 2 + 3 + 4 = 10

可以看出result[::2]和a相等，而result[1::2]和np.add.accumulate(a)相等。

ufunc函数对象的outer()方法等同于如下程序：

a.shape += (1,)*b.ndim
(a,b)
a = a.squeeze()

其中squeeze()方法剔除数组a中长度为1的轴。让我们看一个例子：

np.multiply.outer([1, 2, 3, 4, 5], [2, 3, 4])
array([[ 2, 3, 4],
[ 4, 6, 8],
[ 6, 9, 12],
[ 8, 12, 16],
[10, 15, 20]])

可以看出通过outer()计算的结果是如下乘法表：

*| 2 3 4
------------
1| 2 3 4
2| 4 6 8
3| 6 9 12
4| 8 12 16
5|10 15 20

如果将这两个数组按照等同程序一步一步地进行计算，就会发现乘法表最终是通过广播的方式计算出来的。