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 中的值并没有改变，而是新创建了一个数组保存结果。如果希望将 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

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

比较了一下 numpy.math 和 Python标准库的 math.sin 的计算速度：:

import time
import math
import numpy as npx = [i * 0.001 for i in range(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 range(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()
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 为了同时支持数组和单个值的计算，其 C 语言的内部实现要比math.sin复杂很多，同样在 Python 级别进行循环的话，就会看出其中的差别了。此外，numpy.sin返回的数的类型和math.sin返回的类型有所不同，math.sin返回的是Python的标准float类型，而numpy.sin则返回一个numpy.float64类型：

>>> type(math.sin(0.5))
<type 'float'>
>>> type(np.sin(0.5))
<type 'numpy.float64'>

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

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

add 函数返回一个新的数组，此数组的每个元素都为两个参数数组的对应元素之和。它接受第3个参数指定计算结果所要写入的数组，如果指定的话，add 函数就不再产生新的数组。

由于 Python 的操作符重载功能，计算两个数组相加可以简单地写为 a+b，而np.add(a,b,a) 则可以用 a+=b 来表示。下面是数组的运算符和其对应的 ufunc 函数的一个列表，注意除号"/"的意义根据是否激活__future__.division有所不同。

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函数对数组进行计算了。让我们来看一个例子。

三角波的样子如下图所示：

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

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

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

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):x = x - int(x) # 三角波的周期为1，因此只取x坐标的小数部分进行计算if x >= c: r = 0.0elif x < c0: r = x / c0 * hcelse: r = (c-x) / (c-c0) * hcreturn r# 用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转换结果。

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

广播

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

让所有输入数组都向其中shape最长的数组看齐，shape 中不足的部分都通过在前面加1补齐
输出数组的shape是输入数组shape的各个轴上的最大值
如果输入数组的某个轴和输出数组的对应轴的长度相同或者其长度为1时，这个数组能够用来计算，否则出错
当输入数组的某个轴的长度为1时，沿着此轴运算时都用此轴上的第一组值
上述4条规则理解起来可能比较费劲，让我们来看一个实际的例子。

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

>>> 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
array([0, 1, 2, 3, 4])
>>> b.shape
(5,)

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

>>> c = a + b
>>> c
array([[ 0,  1,  2,  3,  4],[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]])
>>> c.shape
(6, 5)

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

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

这样加法运算的两个输入数组的 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
array([[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],[0, 1, 2, 3, 4]])

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

>>> a = a.repeat(5, axis=1)
>>> a
array([[ 0,  0,  0,  0,  0],[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 就可以按对应元素进行相加运算了。

当然，numpy 在执行 a+b 运算时，其内部并不会真正将长度为1的轴用repeat函数进行扩展，如果这样做的话就太浪费空间了。

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

>>> x,y = np.ogrid[0:5,0:5]
>>> x
array([[0],[1],[2],[3],[4]])
>>> y
array([[0, 1, 2, 3, 4]])

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

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

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

>>> x, y = np.ogrid[0:1:4j, 0:1:3j]
>>> x
array([[ 0.        ],[ 0.33333333],[ 0.66666667],[ 1.        ]])
>>> y
array([[ 0. ,  0.5,  1. ]])

ogrid 为什么不是函数

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

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

import numpy as np
from enthought.mayavi import mlabx, 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库快速绘制如下图所示的3D曲面。

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])
>>> 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]]

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

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的轴。如果你看不太明白这个等同程序的话，栗子：

>>> 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
# 5