>>> from env_helper import info; info()

页面更新时间： 2023-04-15 21:23:28
运行环境：
    Linux发行版本: Debian GNU/Linux 12 (bookworm)
    操作系统内核: Linux-6.1.0-7-amd64-x86_64-with-glibc2.36
    Python版本: 3.11.2

1.2. ufunc运算¶

ufunc是universal function的缩写,它是一种能对数组的每个元素进行操作的函数。NumPy内置的许多ufunc函数都是在C语言级别实现的,因此它们的计算速度非常快。让我们来看一个例子:

>>> import numpy as np
>>> 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.44929360e-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.44929360e-16])

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

>>> 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 np
>>> x = [i * 0.001 for i in range(1000000)]
>>> start = time.time()
>>> for i, t in enumerate(x):
>>>     x[i] = math.sin(t)
>>>
>>> t1 = time.time() - start
>>> print("math.sin:", t1)
>>>
>>> x = [i * 0.001 for i in range(1000000)]
>>> x = np.array(x)
>>> start = time.time()
>>> np.sin(x,x)
>>> t2 = time.time() - start
>>> print("numpy.sin:", t2)

math.sin: 0.16526269912719727
numpy.sin: 0.01620030403137207

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

>>> t1 / t2

10.201209730827532

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

>>> x = [i * 0.001 for i in range(1000000)]
>>> start = time.time()
>>> for i, t in enumerate(x):
>>>     x[i] = np.sin(t)
>>> print("numpy.sin loop:", time.time() - start)

numpy.sin loop: 0.8658764362335205

请注意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))

float

>>> type(np.sin(0.5))

numpy.float64

通过上面的例子我们了解了如何最有效率地使用math库和numpy库中的数学函数。因为它们各有长短,因此在导入时不建议使用*号全部载入,而是应该使用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

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函数对数组进行计算了。让我们来看一个例子。我们想用一个分段函数描述三角波,三角波的样子如图2.4所示: 图 2.4 - 三角波可以用分段函数进行计算我们很容易根据上图所示写出如下的计算三角波某点y坐标的函数:

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

显然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.0
>>>         elif x < c0: r = x / c0 * hc
>>>         else: r = (c-x) / (c-c0) * hc
>>>         return 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)将其转换为双精度浮点数组。

1.2.1. 广播¶

当我们使用ufunc函数对两个数组进行计算时,ufunc函数会对这两个数组的对应元素进行计算,因此它要求这两个数组有相同的大小(shape相同)。如果两个数组的shape不同的话,会进行如下的广播 (broadcasting)处理: 1. 让所有输入数组都向其中shape最长的数组看齐,shape中不足的部分都通过在前面加1补齐 2. 输出数组的shape是输入数组shape的各个轴上的最大值 3. 如果输入数组的某个轴和输出数组的对应轴的长度相同或者其长度为1时,这个数组能够用来计算,否则出错 4. 当输入数组的某个轴的长度为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. ]])

1.2.2. ogrid为什么不是函数¶

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

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

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创建的三维曲面

1.2.3. ufunc的方法¶

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

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

例如:

>>> np.add.reduce([1,2,3]) # 1 + 2 + 3

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

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

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

可以看出result[::2]和a相等,而result[1::2]和np.add.accumulate(a)相等。 outer 方法,.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

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

1.1. NumPy-快速处理数据

1.3. 矩阵运算

Python 3 教程 文档

1.2. ufunc运算¶

1.2.1. 广播¶

1.2.2. ogrid为什么不是函数¶

1.2.3. ufunc的方法¶

Python 3 教程文档