ufunc 是 universal function （ universal ： 通用的，一般的）的缩写，它是一种能对数组的每个元素进行操作的函数。 NumPy内置的许多 ufunc 函数底层都是通过C语言实现的，因此它们的计算速度非常快。

ufunc运算

让我们来看一个例子，先准备一个数组：

对数组x中的每个元素进行正弦计算，返回一个同样大小的新数组：

先用 linspace 产生一个从 0 到 2*PI 的等距离的10个数,然后将其传递给 sin 函数。

由于 np.sin 是一个 ufunc 函数,因此它对 x 中的每个元素求正弦值,然后将结果返回,并且赋值给 y 。 计算之后 x 中的值并没有改变，而是新创建了一个数组保存结果。

NumPy 提供了一种节省内存的方式。 如果我们希望将 sin 函数所计算的结果直接覆盖到数组 x 上去的话，可以将要被覆盖的数组作为第二个参数传递给 ufunc函数。

例如:

查看一下内存中 t 与 x 的地址ID：

sin 函数的第二个参数也是x，那么它所做的事情就是对x中的每给值求正弦值，并且把结果保存到x中的对应的位置中。 此时函数的返回值仍然是整个计算的结果，只不过它就是x，因此两个变量的id是相同的(变量t和变量x指向同一块内存区域)。

为了对运算速度进行说明，使用下面的程序比较了 numpy.math 和 Python 标准库中 math.sin 的计算速度:

在我的电脑上计算100万次正弦值，numpy.sin 比 math.sin 快10倍多：

这得利于 numpy.sin 在C语言级别的循环计算。 numpy.sin 同样也支持对单个数值求正弦，例如: numpy.sin(0.5) 。 不过值得注意的是，对单个数的计算 math.sin 则比 numpy.sin 快得多了，让我们看下面这个测试程序:

请注意 numpy.sin 的计算速度只有 math.sin 的1/5。

这是因为 numpy.sin 为了同时支持数组和单个值的计算，其C语言的内部实现要比 math.sin 复杂很多, 如果我们同样在Python级别进行循环的话，就会看出其中的差别了。

此外，numpy.sin 返回的数的类型和 math.sin 返回的类型有所不同， math.sin 返回的是Python的标准 float 类型，而 numpy.sin 则返回一个 numpy.float64 类型:

通过上面的例子我们了解了如何最有效率地使用 math 库和 numpy 库中的数学函数。 因为它们各有长短，因此在导入时不建议使用 * 号全部载入， 而是应该使用 import numpy as np 的方式载入，这样我们可以根据需要选择合适的函数调用。

运算符和其对应的 ufunc 函数

NumPy中有众多的 ufunc 函数为我们提供各式各样的计算。 除了 sin 这种单输入函数之外，还有许多多个输入的函数， add函数就是一个最常用的例子。先来看一个例子:

add 函数返回一个新的数组，此数组的每个元素都为两个参数数组的对应元素之和。 它接受第3个参数指定计算结果所要写入的数组，如果指定的话， add 函数就不再产生新的数组。 由于Python的操作符重载功能，计算两个数组相加可以简单地写为 a+b ，而 np.add(a,b,a) 则可以用 a+=b 来表示。 下面是数组的运算符和其对应的 ufunc 函数的一个列表， 注意除号 ''/'' 的意义根据是否激活 __future__.division 有所不同。

数组对象支持这些操作符，极大地简化了算式的编写，不过要注意如果算式很复杂， 并且要运算的数组很大，会因为产生大量的中间结果而降低程序的运算效率。 例如:假设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所示:

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

我们很容易根据上图所示写出如下的计算三角波某点y坐标的函数:

三角波的周期为1，因此只取x坐标的小数部分进行计算

显然 triangle_wave 函数只能计算单个数值，不能对数组直接进行处理。 我们可以用下面的方法先使用列表包容(List comprehension)，计算出一个 list ， 然后用 array 函数将列表转换为数组：

这种做法每次都需要使用列表包容语法调用函数，对于多维数组是很麻烦的。 让我们来看看如何用 frompyfunc 函数来解决这个问题:

frompyfunc 的调用格式为 frompyfunc(func, nin, nout) ，其中 func 是计算单个元素的函数， nin 是此函数的输入参数的个数， nout 是此函数的返回值的个数。虽然 triangle_wave 函数有4个参数， 但是由于后三个c, c0, hc在整个计算中值都是固定的，因此所产生的 ufunc 函数其实只有一个参数。 为了满足这个条件，我们用一个 lambda 函数对 triangle_wave 的参数进行一次包装。 这样传入 frompyfunc 的函数就只有一个参数了。这样做，效率并不是太高，另外还有一种方法:

我们通过函数 triangle_func 包装三角波的三个参数，在其内部定义一个计算三角波的函数 trifunc ， trifunc 函数在调用时会采用 triangle_func 的参数进行计算。 最后 triangle_func 返回用 frompyfunc 转换结果。

值得注意的是用 frompyfunc 得到的函数计算出的数组元素的类型为 object ， 因为 frompyfunc 函数无法保证Python函数返回的数据类型都完全一致。 因此还需要再次 y2.astype(np.float64) 将其转换为双精度浮点数组。

广播 （broadcasting）

broadcasting 的翻译，很难理解其意义。

对于英文来讲，其实是 broad + cast 两个单词的拼写。 这两个单词放在一起恰巧是英文单词 broadcast ，但据此将其翻译为“广播”，也应理解其意义并非“广泛传播”。

cast 在编程中非常重要，中文一般称为“造型”，一般是指程序中数据类型的转换。 而在 NumPy 中， 此处的操作并非是数据类型的转换，而是数组长度的对齐，以满足数组计算。 而这个对齐都是由短往长变的，即为 broad 的意义。

当我们使用 ufunc 函数对两个数组进行计算时， ufunc 函数会对这两个数组的对应元素进行计算， 因此它要求这两个数组有相同的大小(shape相同)。如果两个数组的 shape 不同的话， 会进行如下的广播(broadcasting)处理:

1. 让所有输入数组都向其中 shape 最长的数组看齐， shape 中不足的部分都通过在前面加1补齐
2. 输出数组的 shape 是输入数组 shape 的各个轴上的最大值
3. 如果输入数组的某个轴和输出数组的对应轴的长度相同或者其长度为1时，这个数组能够用来计算，否则出错
4. 当输入数组的某个轴的长度为1时,沿着此轴运算时都用此轴上的第一组值

上述4条规则理解起来可能比较费劲,让我们来看一个实际的例子。 先创建一个二维数组 a ,其shape为 (6,1) :

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

计算a和b的和,得到一个加法表，它相当于计算a,b中所有元素组的和， 得到一个 shape 为 (6,5) 的数组:

由于a和b的 shape 长度(也就是 ndim 属性)不同，根据规则1，需要让b的 shape 向a对齐， 于是将b的 shape 前面加1，补齐为(1,5)。相当于做了如下计算:

这样加法运算的两个输入数组的 shape 分别为 (6,1) 和 (1,5) ，根据规则2， 输出数组的各个轴的长度为输入数组各个轴上的长度的最大值，可知输出数组的 shape 为 (6,5) 。 由于b的第0轴上的长度为1，而a的第0轴上的长度为6，因此为了让它们在第0轴上能够相加， 需要将b在第0轴上的长度扩展为6，这相当于:

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

经过上述处理之后，a和b就可以按对应元素进行相加运算了。当然，numpy 在执行 a+b 运算时， 其内部并不会真正将长度为1的轴用 repeat 函数进行扩展，如果这样做的话就太浪费空间了。

由于这种广播计算很常用，因此 numpy 提供了一个快速产生如上面a,b数组的方法: ogrid 对象:

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

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

ogrid 为什么不是函数

根据Python的语法，只有在中括号中才能使用用冒号隔开的切片语法，如果 ogrid 是函数的话, 那么这些切片必须使用 slice 函数创建，这显然会增加代码的长度。

利用 ogrid 的返回值，我能很容易计算x, y网格面上各点的值，或者x, y, z网格体上各点的值。 下面是绘制三维曲面 x * exp(x**2 - y**2) 的程序:

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")
mlab.axes(xlabel='x', ylabel='y', zlabel='z')
mlab.outline(pl)

此程序使用 mayavi 的 mlab 库快速绘制如图2.5所示的3D曲面， 关于 mlab 的相关内容将在今后的章节进行介绍。

图 2.5 - 使用 ogrid 创建的三维曲面

ufunc 的方法

ufunc 函数本身还有些方法，这些方法只对两个输入一个输出的 ufunc 函数有效， 其它的 ufunc 对象调用这些方法时会抛出 ValueError 异常。 reduce 方法和Python的 reduce 函数类似，它沿着 axis 轴对 array 进行操作， 相当于将运算符插入到沿 axis 轴的所有子数组或者元素当中。

<op>.reduce (array=, axis=0, dtype=None)

例如:

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

reduceat 方法计算多组 reduce 的结果，通过 indices 参数指定一系列 reduce 的起始和终了位置。 reduceat 的计算有些特别，让我们通过一个例子来解释一下:

对于 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]

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

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

因此上面例子中，结果的每个元素如下计算而得：

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) 相等。 outer 方法， <op>.outer(a,b) 方法的计算等同于如下程序:

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

其中 squeeze 的功能是剔除数组a中长度为1的轴。如果你看不太明白这个等同程序的话， 让我们来看一个例子:

可以看出通过 outer 方法计算的结果是如下的乘法表:

#
2, 3, 4
# 1
# 2
# 3
# 4
# 5

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