ufunc运算

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.44921271e-16])

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

>>> t = np.sin(x,out=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

ufunc函数的返回值仍然是计算的结果,只不过它就是数组x,因此两个数组的id是相同的。

下面比较np.sin和Python标准库的math.sin的计算速度:

02-numpy/numpy_speed_test.py

NumPy计算和Python标准库的计算速度比较

import time
import math
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

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

程序的输出为:

math.sin: 0.779217749742
numpy.sin: 0.0772958574344
numpy.sin loop: 5.25540843878

可以看出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 'float'>
>>> type(np.sin(0.5))
<type 'numpy.float64'>

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

>>> a = np.arange(6.0).reshape(2,3)
>>> a.item(1,2) # 和a[1,2]类似
5.0
>>> type(a.item(1,2)) # item()所返回的是Python的标准float类型
<type 'float'>
>>> type(a[1,2]) # 下标方式获得的是NumPy的float64类型
<type 'numpy.float64'>

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

四则运算

NumPy提供许多ufunc函数,例如计算两个数组之和的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()返回一个数组,它的每个元素都是两个参数数组的对应元素之和。如果没有指定out参数,那么它将创建一个新的数组保存计算结果。如果指定了第三个参数out,则不产生新的数组,而直接将结果保存进指定的数组。

NumPy为数组定义了各种数学运算操作符,因此计算两个数组相加可以简单地写为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

比较和布尔运算

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

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

每个比较运算符也与一个ufunc函数对应,下面是比较运算符和其ufunc函数的对照表:

y = x1 == x2 equal(x1, x2 [, y])
y = x1 != x2 not_equal(x1, x2 [, y])
y = x1 < x2 less(x1, x2, [, y])
y = x1 <= x2 less_equal(x1, x2, [, y])
y = x1 > x2 greater(x1, x2, [, y])
y = x1 >= x2 greater_equal(x1, x2, [, y])

由于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)
>>> a == b
array([False, False,  True, False, False], dtype=bool)
>>> a > b
array([False, False, False,  True,  True], dtype=bool)
>>> np.logical_or(a==b, a>b) # 和 a>=b 相同
array([False, False,  True,  True,  True], dtype=bool)

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

>>> 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()方法[1]。只要数组中有一个值为True,则any()返回True;而只有数组的全部元素都为True,all()才返回True。

>>> np.any(a==b)
True
>>> np.any(a==b) and np.any(a>b)
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。

Footnotes

[1]在NumPy中同时也定义了any()和all()函数。

自定义ufunc函数

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

例如,我们可以用一个分段函数描述三角波,三角波的样子如【图:三角波可以用分段函数进行计算】所示,它分为三段:上升、下降和平坦段。

/tech/static/books/scipy/_images//numpy_frompyfunc.png

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

02-numpy/numpy_frompyfunc.py

用frompyfunc()计算三角波的波形数组

根据【图:三角波可以用分段函数进行计算】,我们很容易写出计算三角波上某点的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()方法将其转换为双精度浮点数组:

>>> run numpy_frompyfunc.py # 在IPython中运行计算三角波的程序
>>> y2.dtype
dtype('object')
>>> y2 = y2.astype(np.float)
>>> y2.dtype
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)
True
>>> np.all(y2==y3)
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
array([[ 0], [10], [20], [30], [40], [50]])
>>> a.shape
(6, 1)

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

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

计算a和b的和,得到一个加法表,它相当于计算两个数组中所有元素组的和,得到一个形状为(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的维数不同,根据规则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的第1轴长度为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就可以按对应元素进行相加运算了。

当然,在执行“a+b”运算时,NumPy内部并不会真正将长度为1的轴用repeat()进行扩展,这样太浪费空间了。

由于这种广播计算很常用,因此NumPy提供了快速产生能进行广播运算的数组的ogrid对象。

>>> x,y = np.ogrid[:5,:5]
>>> x
array([[0], [1], [2], [3], [4]])
>>> y
array([[0, 1, 2, 3, 4]])
mgrid对象的用法和ogrid对象类似,但是它所返回的是进行广播之后的数组。请读者运行“np.mgrid[:5,:5]”试试看。

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

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

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

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

利用ogrid的返回值,我们能很容易计算二元函数在等间距网格上的值。下面是绘制三维曲面f(x,y) = x e^{x^2-y^2}的程序:

02-numpy/numpy_ogrid_mlab.py

用ogird产生二维坐标网格,计算三维空间的曲面

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)
mlab.show()

此程序使用Mayavi的mlab模块快速绘制如【图:使用ogrid创建的三维曲面】所示的3D曲面,关于Mayavi的相关内容将在今后的章节进行介绍。

/tech/static/books/scipy/_images//numpy_ogrid_mlab.png

使用ogrid创建的三维曲面

如果已经有了多个表示不同轴上的取值的一维数组。想计算出它们所构成的网格上的每点的函数值,可以使用ix_():

>>> x = np.array([0, 1, 4, 10])
>>> y = np.array([2, 3, 8])
>>> gy, gx = np.ix_(y, x)
>>> gx
array([[ 0,  1,  4, 10]])
>>> gy
array([[2], [3], [8]])
>>> gx+gy # 广播运算
array([[ 2,  3,  6, 12],
       [ 3,  4,  7, 13],
       [ 8,  9, 12, 18]])

在上面的例子中,通过ix_()将数组x和y转换成能进行广播运算的二维数组。注意数组y对应广播运算结果中的第0轴,而数组x与第1轴对应。

ufunc的方法

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

reduce()方法和Python的reduce()函数类似,它沿着axis参数指定的轴对数组进行操作,相当于将<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+2+3),(4+5+6)
array([ 6, 15])

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

>>> 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参数指定一系列的起始和终了位置。它的计算有些特别,让我们通过例子详细解释一下:

>>> 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] = <op>.reduce(a[indices[i]:indices[i+1]])
else:
    result[i] = a[indices[i]]

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

<op>.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)相等。使用reduceat()可以对数组中的多个片段同时进行reduce运算。

ufunc函数对象的outer()方法对其两个参数数组的每两对元素的组合进行运算。若数组a的维数为M,数组b的维数为N,则ufunc函数op的outer()方法对a、b数组计算所生成的数组c的维数为M+N。c的形状是a、b的形状的结合。例如a的形状为(2,3),b的形状为(4,5),则c的形状为(2,3,4,5)。让我们看一个例子:

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

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

內容目录

上一个主题

ndarray对象

下一个主题

多维数组的下标存取

本页

loading...