参考文献:

  1. NumPy Reference

  2. NumPy User Guide

  3. Python Numpy Tutorial (with Jupyter and Colab)

  4. 菜鸟教程

  5. 挑战NumPy100关

  我们知道 Python 用 list 结构保存一组值，可以当作数组用，列表的元素可以是任何对象，因此列表中保存的是对象的指针。对数值运算而言，这种结构比较耗费内存和CPU。此外还有 array 模块，能直接保存数值，但不支持多维数组，也没有各种运算函数，也不适合做数值运算。

  NumPy 的诞生弥补了这些不足，它提供 ndarray(n-dimensional array object) 和 ufunc(universal function object) 两种基本对象。以下为本文的导入方式和 NumPy 版本号:

1
2
3
>>> import numpy as np
>>> np.__version__
'1.16.4'

ndarray 对象

创建

  通过给 array()函数传递 Python 的序列对象来创建数组。

1
2
3
>>> a = np.array([1, 2, 3, 4])
>>> b = np.array((5, 6, 7, 8))
>>> c = np.array([[1, 2, 3, 4], [4, 5, 6, 7], [7, 8, 9, 10]])

  数组的形状通过 shape 属性获得 它是一个描述数组各个轴的长度的元组(tuple):

1
2
a.shape b.shape c.shape
(4,)    (4,)    (3,4)

  下面的例子将数组c 的 shape 属性改为(4,3)，并不是对数组进行转置，而是改变每个轴大小，数组元素在内存中的位置并没有改变。此外，当设置某个轴元素个数为-1时，将自动计算此轴长度。

1
2
3
4
5
6
>>> c.shape = 4, 3
>>> c
array([[ 1,  2,  3],
       [ 4,  4,  5],
       [ 6,  7,  7],
       [ 8,  9, 10]])

  使用数组的 reshape()方法，可以创建指定形状的新数组，而原数组形状保持不变。数组 a 和 d 共享数据存储空间，修改任意一个数组元素会同时修改另一个数组的内容。

1
2
3
4
5
6
>>> d = a.reshape(2, 2)
>>> a[1] = 100
     d             a
------------ --------------
[[  1, 100], [1, 100, 3, 4]
 [  3,   4]]

元素类型

  数组的元素类型可以通过 dtype 属性获得。

1
2
>>> c.dtype
dtype('int32')

  可以通过 dtype 参数在创建数组时指定元素类型，注意 float 是 64 位的双精度浮点类型，complex 是 128 位的双精度复数类型。

1
2
3
4
5
>>> ai32 = np.array([1, 2, 3, 4], dtype=np.int32)
>>> af = np.array([1, 2, 3, 4], dtype=float)
>>> ac = np.array([1, 2, 3, 4], dtype=complex)
ai32.dtype af.dtype ac.dtype
int32      float64  complex128

  在需要指定 dtype 参数时，也可以传递一个字符串来表示元素的数值类型。NumPy 中的每个数值类型都有几种字符串表示方式，字符串和类型之间的对应关系储存在 typeDict 字典中。下面例子获得与 float64 类型对应的所有键值:

1
2
>>> [key for key, value in np.typeDict.items() if value is np.float64]
['double', 'd', 12, 'float64', 'Float64', 'f8', 'float_', 'float']

  完整的类型列表通过下面的语句得到。

1
2
3
4
5
6
7
8
set(np.typeDict.values())
---
{numpy.bool_, numpy.bytes_, numpy.complex128, numpy.complex128,
 numpy.complex64, numpy.datetime64, numpy.float16, numpy.float32,
 numpy.float64, numpy.float64, numpy.int16, numpy.int32,
 numpy.int32, numpy.int64, numpy.int8, numpy.object_,
 numpy.str_, numpy.timedelta64, numpy.uint16, numpy.uint32,
 numpy.uint32, numpy.uint64, numpy.uint8, numpy.void}

  上面显示的数值类型与数组的 dtype 属性是不同的对象。通过 dtype 对象的 type 属性可以获得与其对应的数值类型:

1
2
>>> c.dtype.type
<class 'numpy.int32'>

  下面创建一个 16 位(\(-2^{15}\sim2^{15}-1\))的符号整数对象，计算 \(200*200\) 会溢出，得到一个负数 \(200*200-2^{16}\) 。

1
2
3
4
>>> a = np.int16(200)
>>> a*a
__main__:1: RuntimeWarning: overflow encountered in short_scalars
-25536

  NumPy 的数值对象的运算速度比 Python 的内置类型的运算速度慢很多，如果程序中需要大量地对单个数值运算，应避免使用 NumPy 的数值对象。

1
2
3
4
5
6
7
v1 = 3.14
v2 = np.float64(v1)
%timeit v1*v1
%timeit v2*v2
---
85.1 ns ± 4.1 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)
159 ns ± 10.3 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)

  使用 astype() 方法可以对数组元素类型进行转换:

1
2
3
4
>>> t1 = np.array([1, 2, 3, 4], dtype=np.float)
>>> t2 = np.array([1, 2, 3, 4], dtype=np.complex)
>>> t3 = t1.astype(np.int32)
>>> t4 = t2.astype(np.complex64)

自动生成数组

  前面的例子都是先创建一个 Python 的序列对象，然后通过array()将其转换为数组，这样做显然效率不高。因此NumPy 提供了很多专门用于创建数组的函数。

  arange()类似于内置函数 range()，通过指定开始值、终值和步长来创建表示等差数列的一维数组。

1
2
>>> np.arange(0, 1, 0.1)
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])

  linspace()通过指定开始值、终值和元素个数来创建表示等差数列的一维数组，可通过默认位 True 的 endpoint 参数指定是否包含终值。

1
2
3
4
5
>>> np.linspace(0, 1, 10)
array([0.        , 0.11111111, 0.22222222, 0.33333333, 0.44444444,
       0.55555556, 0.66666667, 0.77777778, 0.88888889, 1.        ])
>>> np.linspace(0, 1, 10, endpoint=False)
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])

  logspace()和 linspace()类似，不过创建的是等比数列。下面的例子为从 \(10^0\) 到 \(10^2\)、有 5 个元素的等比数列。

1
2
3
>>> np.logspace(0, 2, 5)
array([  1.        ,   3.16227766,  10.        ,  31.6227766 ,
       100.        ])

  基数可以通过 base 参数指定，默认值为 10。下面创建一个比例为 \(2^{1/12}\) 的等比数组，其比值为音乐中相差半音的两个音阶之间的频率比值，可用于计算一个八度中所有半音的频率。

1
2
3
4
>>> np.logspace(0, 1, 12, base=2, endpoint=False)
array([1.        , 1.05946309, 1.12246205, 1.18920712, 1.25992105,
       1.33483985, 1.41421356, 1.49830708, 1.58740105, 1.68179283,
       1.78179744, 1.88774863])

  zeros()、ones()、empty()，eye()可以创建指定形状和类型的数组。其中empty()只分配数组所使用的内存，不对数组元素进行初始化操作，因此它的运行速度是最快的。zeros()将元素初始化为 0，ones()将数组元素初始化为 1。

1
2
3
4
5
>>> np.empty((2,3), np.int)
array([[   4063294,    2097214, 1469448192],
       [     32760, 1469454912,      32760]])
>>> np.zeros(4, np.int)
array([0, 0, 0, 0])

  full()将数组元素初始化为指定的值:

1
2
>>> np.full(4, np.pi)
array([3.14159265, 3.14159265, 3.14159265, 3.14159265])

  此外，zeros_like()、ones_like()、empty_like()等函数可创建与参数数组的形状及类型相同的数组。因此，zeros_like(a)和 zeros(a.shape, a.dtype)的效果相同。

1
2
3
>>> a = np.arange(3, dtype=float)
>>> np.zeros_like(a)
array([0., 0., 0.])

  frombuffer()、fromstring()、fromfile()等函数可以从字节序列或文件创建数组，下面以 fromstring()为例介绍它们的用法。Python 的字符串实际上是一个字节序列，每个字符占一个字节，因此如果从字符串s 创建一个 8 位的整数数组，所得到的数组正好就是字符串中每个字符的 ASCII 编码：

1
2
3
4
>>> s = "abcdefgh"
>>> np.fromstring(s, dtype=np.int8)
__main__:1: DeprecationWarning: The binary mode of fromstring is deprecated, as it behaves surprisingly on unicode inputs. Use frombuffer instead
array([ 97,  98,  99, 100, 101, 102, 103, 104], dtype=int8)

  如果从字符串 s 创建 16 位的整数数组，那么两个相邻的字节就表示一个整数，把字节 98 和字节 97 当作一个16 位的整数，它的值就是\(98*256+97=25185\)。可以看出，16 位的整数是以低位字节在前(little-endian)的方式保存在内存中的。

1
2
>>> np.fromstring(s, dtype=np.int16)
array([25185, 25699, 26213, 26727], dtype=int16)

  fromfunction()通过函数创建数组，第一个参数为计算数组元素的函数，第二个参数为 shape。

1
2
>>> np.fromfunction(lambda i: i % 4 + 1, (10,), dtype=int)
array([1, 2, 3, 4, 1, 2, 3, 4, 1, 2], dtype=int32)

存取元素

  NumPy 切片语法类似于 Python 列表的标准切片语法，即数组 x 切片获取方式：x[start:stop:step]，因此将数组颠倒即 x[::-1]。除了使用切片存取元素之外，NumPy 还提供了整数列表和布尔数组等下标存取方式。

1
2
3
4
5
6
>>> x = np.arange(10, 1, -1)
>>> a = x[[3, 3, -3, 8]]
>>> x
array([10,  9,  8,  7,  6,  5,  4,  3,  2])
>>> a
array([7, 7, 4, 2])

  整数序列下标也可以用来修改元素的值：

1
2
3
>>> x[[3, 5, 1]] = -1, -2, -3
>>> x
array([10, -3,  8, -1,  6, -2,  4,  3,  2])

  当下标是多维数组时，得到的也是多维数组：

1
2
3
4
5
6
>>> x[np.array([[3, 3, 1, 8], [3, 3, -3, 8]])]
array([[-1, -1, -3,  2],
       [-1, -1,  4,  2]])
>>> x[[3, 3, 1, 8, 3, 3, -3, 8]].reshape(2, 4)
array([[-1, -1, -3,  2],
       [-1, -1,  4,  2]])

  下面是 bool 数组的例子：

1
2
3
4
5
6
7
8
>>> x = np.arange(5, 0, -1)
>>> x[np.array([True, False, True, True, False])]
array([5, 3, 2])
>>> x = np.random.randint(0, 10, 6)
>>> x
array([0, 8, 5, 7, 7, 0])
>>> x[x > 5]
array([8, 7, 7])

多维数组

  多维数组与一维数组类似，NumPy 采用元组作为数组的下标。

1
2
3
4
5
6
7
8
9
10
11
12
>>> a = np.arange(0, 60, 10).reshape(-1, 1) + np.arange(0, 6)
>>> a
array([[ 0,  1,  2,  3,  4,  5],
       [10, 11, 12, 13, 14, 15],
       [20, 21, 22, 23, 24, 25],
       [30, 31, 32, 33, 34, 35],
       [40, 41, 42, 43, 44, 45],
       [50, 51, 52, 53, 54, 55]])
---
a[0, 3:5]     a[4:, 4:]        a[:, 2]                         a[2::2, ::2]
array([3, 4]) array([[44, 45], array([ 2, 12, 22, 32, 42, 52]) array([[20, 22, 24],
                     [54, 55]])                                       [40, 42, 44]])

  如果下标元组中只包含整数和切片，那么得到的数组和原始数组共享数据，它是原数组的视图。下面的例子中，数组 b 是 a 的视图，它们共享数据，因此修改 b[0]时，数组 a 中对应的元素也被修改：

1
2
3
4
>>> b = a[0, 3:5]
>>> b[0] = -b[0]
>>> a[0, 3:5]
array([-3,  4])

  因为数组的下标是一个元组，所以我们可以将下标元组保存起来，用同一个元组存取多个数组。下面的例子中，a[idx]和 a[::2, 2:]相同，a[idx][idx]和 a[::2, 2:][::2, 2:]相同。

1
2
3
4
5
6
7
8
>>> idx = slice(None, None, 2), slice(2, None)
>>> a[idx]
array([[ 2, -3,  4,  5],
       [22, 23, 24, 25],
       [42, 43, 44, 45]])
>>> a[idx][idx]
array([[ 4,  5],
       [44, 45]])

  用 Python 内置的 slice() 函数创建下标比较麻烦，因此 NumPy 提供了一个 s_对象来帮助我们创建数组下标，s_实际上是 IndexExpression 类的一个对象：

1
2
>>> np.s_[::2, 2:]
(slice(None, None, 2), slice(2, None, None))

  在多维数组的下标元组中，也可以使用整数元组或列表、整数数组和布尔数组。当在下标中使用这些对象时，所获得的数据是原始数据的副本，因此修改结果数组不会改变原始数组。

  当所有轴都用形状相同的整数数组作为下标时，得到的数组和下标数组的形状相同：

1
2
3
4
5
6
7
>>> x = np.array([[0, 1], [2, 3]])
>>> y = np.array([[-1, -2], [-3, -4]])
>>> a[x, y]
array([[ 5, 14],
       [23, 32]])
# or
>>> a[(0, 1, 2, 3), (-1, -2, -3, -4)].reshape(2, 2)

结构数组

  在C语言中可以通过 struct 关键字定义结构类型，结构中的字段占据连续的内存空间。类型相同的两个结构体所占用的内存大小相同，因此可以很容易定义结构数组。与C语言类似，NumPy 中也很容易对这种结构数组进行操作。

1
2
3
4
5
persontype = np.dtype({ #1
    'names':['name', 'age', 'weight'],
    'formats':['S32','i', 'f']}, align=True)
a = np.array([("Zhang", 32, 75.5), ("Wang", 24, 65.2)], #2
             dtype = persontype)

  #1 先创建一个 dtype 对象 persontype，它的参数是一个描述结构类型的各个字段的字典。字典有两个键：'names'和'formats'。每个键对应的值都是一个列表。'names'定义结构中每个字段的名称，而'formats'则定义每个字段的类型。这里使用类型字符串定义字段类型：

  • 'S32' ：长度为32 字节的字符串类型，由于结构中每个元素的大小必须固定，因此需要指定字符串的长度。
  • 'i' ：32 bit 的整数类型，相当于 np.int32。
  • 'f' ：32 bit 的单精度浮点数类型，相当于 np.float32。

  #2 然后调用array()创建数组，通过 dtype 参数指定所创建数组的元素类型为 persontype。数组a 的元素类型：

1
2
>>> a.dtype
dtype([('name', 'S32'), ('age', '<i4'), ('weight', '<f4')], align=True)

  因此还可以用多个元组的列表描述结构的类型：

1
persontype = np.dtype([('name', '|S32'), ('age', '<i4'), ('weight', '<f4')])

  类型字符串前面的'|'、'<'、'>'等字符表示字段值的字节顺序：

  • | ：忽视字节顺序
  • < ：低位字节在前，即小端模式(little endian)
  • > ：高位字节在前，即大端模式(big endian)

  下面的例子是两种修改方式：

1
2
3
4
5
6
7
8
>>> c = a[1]
>>> c["name"] = "Li"
>>> a[1]["name"]
b'Li'
>>> b = a["age"]
>>> b[0] = 40
>>> a[0]["age"]
40

  通过 a.tostring() 或 a.tofile() 方法，可以将数组 a 以二进制的方式转换成字符串或写入文件：

1
a.tofile("test.bin")

  在 IPyton 中运行以下程序读取 test.bin 中数据。%%file 为 IPython 的 magic command, 它将单元格中文本保存成文件read_struct_array.c

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
%%file read_struct_array.c
#include <stdio.h>

struct person {
    char name[30];
    int age;
    float weight;
}p[3];

int main() {
    FILE *fp = fopen("test.bin", "rb");
    fread(p, sizeof(struct person), 2, fp);
    fclose(fp);
    for (int i = 0; i < 2; ++i) {
        printf("%s %d %f\n", p[i].name, p[i].age, p[i].weight);
    }
    return 0;
}

  在 IPython 中可以通过!执行系统命令

1
2
3
4
5
!gcc read_struct_array.c -o read_struct_array.exe
!read_struct_array.exe
---
Zhang 40 75.500000
Li 24 65.199997

  结构类型中可以包括其他的结构类型，例如：

1
2
>>> np.dtype([('f1', [('f2', '<i2')])])
dtype([('f1', [('f2', '<i2')])])

  当某个字段的类型为数组时，用元组的第三个参数表示其形状：

1
2
>>> np.dtype([('f0', '<i4'), ('f1', '<f8', (2, 3))])
dtype([('f0', '<i4'), ('f1', '<f8', (2, 3))])

  用下面的字典参数也可以定义结构类型，字典的键为结构中的字段名，值为字段的类型描述，但是由于字典的键是没有顺序的，因此字段的顺序需要在类型描述中给出。类型描述是一个元组，它的第二个值给出字段的以字节为单位的偏移量。

1
2
>>> np.dtype({'name':('S25', 0), 'age':(np.uint8, 25)})
dtype([('name', 'S25'), ('age', 'u1')])

ufunc 函数

  ufunc 是一种能对数组每个元素进行运算的函数。NumPy 中许多函数都是由 C 语言编写，计算速度非常快。

四则运算

表达式 ufunc 函数
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])

  一些例子：

1
2
3
4
5
6
7
8
9
10
>>> x1 = np.arange(9.0).reshape((3, 3))
>>> x2 = np.arange(3.0)
>>> np.multiply(x1, x2)
array([[ 0.,  1.,  4.],
       [ 0.,  4., 10.],
       [ 0.,  7., 16.]])
>>> np.divmod(np.arange(5), 3)
(array([0, 0, 0, 1, 1], dtype=int32), array([0, 1, 2, 0, 1], dtype=int32))
>>> np.gcd(np.arange(6), 20)
array([20,  1,  2,  1,  4,  5])

比较运算和布尔运算

表达式 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 中的 bool 运算使用 and、or 和 not 等关键字，它们无法被重载，因此数组的 bool 运算只能通过相应的 ufunc 函数进行。这些函数名都以“logical_”开头，在 IPython 中使用自动补全找到它们：

1
2
np.logical # Press Tab
np.logical_and np.logical_not np.logical_or np.logical_xor

  位运算函数包括 bitwise_and、bitwise_or、bitwise_xor、invert、left_shift 和 right_shift。

  对整数 0，在 32 位符号整数中按位取反的结果是 0xFFFFFFFF，这个值表示 -1。而在 8 位无符号整数中结果为 0xFF，为 255。

1
2
3
4
>>> ~ np.arange(5)
array([-1, -2, -3, -4, -5], dtype=int32)
>>> ~ np.arange(5, dtype=np.uint8)
array([255, 254, 253, 252, 251], dtype=uint8)

  left_shift 例子如下：

1
2
3
4
>>> np.binary_repr(5)
'101'
>>> np.left_shift(5, [1,2,3])
array([10, 20, 40], dtype=int32)

自定义 ufunc 函数

  如图，我们用一个分段函数计算三角波上某点 y 坐标：

1
2
3
4
5
6
def triangle_wave(x, c, c0, hc):
    x = x - int(x)
    if x >= c: r = 0.0
    elif x < c0: r = x / c0 * hc
    else: r = (c - x) / (c - c0) * hc
    return r

  若考虑以下方法，每次都要用列表推导式调用函数，对于多维数组很麻烦

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

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

1
frompyfunc(func, nin, nout)

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

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

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

1
2
y2.dtype   y2.astype(np.float).dtype
dtype('O') dtype('float64')

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

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

  验证结果如下：

1
2
np.all(y1 == y2) np.all(y2 == y3)
True             True

广播 (Broadcasting)

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

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

  我们举一个例子，先创建二维数组 a，形状为(6, 1)：

1
2
3
4
5
6
7
8
>>> a = np.arange(0, 60, 10).reshape(-1, 1)：
a      a.shape
[[ 0], (6, 1)
 [10],
 [20],
 [30],
 [40],
 [50]]

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

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

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

1
2
3
4
5
6
7
8
>>> 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)，相当于做了如下计算：

1
2
3
>>> 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，这相当于：

1
2
3
4
5
6
7
8
>>> 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]]

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

1
2
3
4
5
6
7
8
>>> 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 对象。

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

  mgrid 对象的用法和 ogrid 对象类似，但是它所返回的是进行广播之后的数组：

1
2
3
4
5
6
7
>>> 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(开始值, 结束值, 长度)类似。
1
2
3
4
5
6
>>> x, y = np.ogrid[:1:4j, :1:3j]
x              y
[[0.        ], [[0. , 0.5, 1. ]]
 [0.33333333],
 [0.66666667],
 [1.        ]]

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

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

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

1
2
3
4
5
6
7
8
>>> a = np.arange(4)
>>> a[None, :]
array([[0, 1, 2, 3]])
>>> a[:, None]
array([[0],
       [1],
       [2],
       [3]])

  因此可以实现广播运算：

1
2
3
4
5
6
>>> 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_()将两个一维数组转换成可广播的二维数组：

1
2
3
4
5
>>> 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]]

  注意数组 y 对应广播运算结果中的第 0 轴，而数组 x 与第 1 轴对应。ix_()的参数可以是 N 个一维数。

ufunc 的方法

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

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

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

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

1
2
3
4
5
>>> 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 参数指定一系列的起始和终止位置。它的计算有些特别，参考下面例子：

1
2
3
4
>>> a = np.array([1, 2, 3, 4])
>>> res = np.add.reduceat(a, indices=[0, 1, 0, 2, 0, 3, 0])
>>> res
array([ 1,  2,  3,  3,  6,  4, 10], dtype=int32)

  它按照如下计算得出：

1
2
3
4
if indices[i] < indices[i+1]:
    res[i] = <op>.reduce(a[indices[i]:indices[i+1]])
else:
    res[i] = a[indices[i]]

  而最后一个元素则按照<op>.reduce(a[indices[-1]:])计算得出。

  在例子中，可以看出 res[::2]和 a 相等，而 res[1::2]和 np.add.accumulate(a) 相等。

  ufunc 函数的 outer()方法等同于：

1
2
3
a.shape += (1,)*b.ndim
<op>(a, b)
a = a.squeee()

  其中 squeeze()方法剔除数组 a 中长度为 1 的轴，例如：

1
2
3
4
5
6
>>> 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()计算得到的结果是乘法表。

多维数组的下标存取

下标对象

  多维数组的下标应该是一个长度上与数组的维数相同的元组。如果下标元组的长度比数组的维数大，就会出错。如果小，就需要在下标元组的后面补“:”，使得它的长度与数组维数相同。

  如果下标对象不是元组，NumPy 会首先把它转换为元组。这种转换可能会和用户所希望的不一致。数组 a 是一个三维数组，下面的例子用二维列表 lidx 和一个二维数组 aidx 作为下标，得到的结果是不一样的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
>>> a = np.arange(3 * 4 * 5).reshape(3, 4, 5)
>>> lidx = [[0], [1]]
>>> aidx = np.array(lidx)
>>> a[lidx]
__main__:1: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
array([[5, 6, 7, 8, 9]])
>>> a[aidx]
array([[[[ 0,  1,  2,  3,  4],
         [ 5,  6,  7,  8,  9],
         [10, 11, 12, 13, 14],
         [15, 16, 17, 18, 19]]],


       [[[20, 21, 22, 23, 24],
         [25, 26, 27, 28, 29],
         [30, 31, 32, 33, 34],
         [35, 36, 37, 38, 39]]]])

  下标元组各个元素有如下几种类型：切片、整数、整数数组和布尔数组。如果元素不是这些类型，如列表或元组，就将其转换成整数数组。如果下标元组的所有元素都是切片和整数，那么用它作为下标得到的是原始数组的一个视图，即它和原始数组共享数据存储空间。

整数数组作为下标

  下面讨论下标元组中的元素由切片和整数数组构成的情况。假设整数数组有 \(N_c\) 个，切片有 \(N_s\) 个。\(N_c+N_s\) 为数组的维数 \(D\)。首先这 \(N_c\) 个整数数组必须满足广播条件，假设它们进行广播之后的维数为 \(M\)，形状为\((d_0,d_1,\cdots,d_{M-1})\)。如果 \(N_s\) 为0，即没有切片元素，那么下标得到的结果数组 \(res\) 的形状和整数数组广播之后的形状相同。它的每个元素值可按照下面的公式得出： \[res[i_0,i_1,\cdots,i_{M-1}]=X[ind_0[i_0,i_1,\cdots,i_{M-1}],\cdots,ind_{N_t-1}[i_0,i_1,\cdots,i_{M-1}]]\]

  其中 \(ind0\) 到 \(ind_{N_t-1}\) 为进行广播之后的整数数组。

若只需沿着指定轴通过整数数组获取元素，可以使用numpy.take()函数，运算速度比整数数组下标略快，支持下标越界处理。

1
2
3
4
5
6
7
8
9
10
11
12
>>> i0 = np.array([[1, 2, 1], [0, 1, 0]])
>>> i1 = np.array([[[0]], [[1]]])
>>> i2 = np.array([[[2, 3, 2]]])
>>> b = a[i0, i1, i2]
>>> b
array([[[22, 43, 22],
        [ 2, 23,  2]],

       [[27, 48, 27],
        [ 7, 28,  7]]])
>>> b.shape
(2, 2, 3)

  i0、i1、i2 三个整数数组的 shape 属性分别为(2,3)、(2,1,1)、(1,1,3)，根据广播规则，先在长度不足 3 的 shape 属性前面补 1，使它们的维数相同，广播之后的shape 属性为各个轴的最大值，即三个整数数组广播之后的 shape 属性为(2,2,3)。

  可以使用 broadcast_arrays()查看广播之后的数组：

1
2
3
4
5
6
7
>>> ind0, ind1, ind2 = np.broadcast_arrays(i0, i1, i2)
ind0                 ind1                 ind2
array([[[1, 2, 1],   array([[[0, 0, 0],   array([[[2, 3, 2],
        [0, 1, 0]],          [0, 0, 0]],          [2, 3, 2]],

       [[1, 2, 1],          [[1, 1, 1],          [[2, 3, 2],
        [0, 1, 0]]])         [1, 1, 1]]])         [2, 3, 2]]])

  对于数组 b 中的任意一个元素 b[i,j,k]，它是数组 a 中经过 ind0、ind1 和 ind2 进行下标转换之后的值：

1
2
3
>>> i, j, k = 1, 1, 1
>>> print(b[i, j, k], a[ind0[i, j, k], ind1[i, j, k], ind2[i, j, k]])
28 28

  下面考虑 \(N_s\) 不为 0 的情况。当存在切片下标时，可以细分为两种情况：下标元组中的整数数组之间没有切片，即整数数组只有一个或者是连续的。这时结果数组的 shape 属性为：将原始数组的 shape 属性中整数数组所占据的部分替换为它们广播之后的 shape 属性。例如，假设原始数组 a 的 shape 属性为(3,4,5)，i0 和 i1 广播之后的形状为(2, 2, 3)，则 a[1:3, i0, i1]的形状为(2, 2, 2, 3)：

1
2
3
>>> c = a[1:3, i0, i1]
>>> c.shape
(2, 2, 2, 3)

  其中，数组 c 的 shape 属性中的第一个 2 是切片 “1:3” 的长度，后面的(2,2,3)则是 i0 和 i1 广播之后数组的形状。

  当下标元组中的整数数组不是连续的，结果数组的 shape 属性为整数数组广播之后的形状后面再加上切片元素对应的形状。例如，a[i0, :, i1]的 shape 属性为(2, 2, 3, 4)。其中(2, 2, 3)是 i0 和 i1 广播之后的形状，而 4 是数组 a 第 1 轴的长度：

1
2
3
>>> d = a[i0, :, i1]
>>> d.shape
(2, 2, 3, 4)

布尔数组作为下标

  当使用布尔数组直接作为下标对象或者元组下标对象中有布尔数组时，都相当于用 nonzero()将布尔数组转换成一组整数数组，然后使用整数数组进行下标运算。

  nonzeros(a)返回数组 a 中值不为零的元素的下标，它的返回值是一个长度为 a.ndim(数组 a 的轴数)的元组，元组的每个元素都是一个整数数组，其值为非零元素的下标在对应轴上的值。对于一维布尔数组b1，nonzero(b1)得到的是一个长度为 1 的元组:

若只需沿着指定轴通过整数数组获取元素，可以使用numpy.compress()函数。

1
2
3
>>> b1 = np.array([True, False, True, False])
>>> np.nonzero(b1)
(array([0, 2], dtype=int64),)

  对于二维数组 b2，nonzero(b2)得到的是一个长度为 2 的元组。它的第 0 个元素是数组 a 中值不为 0 的元素的第 0 轴的下标，第 1 个元素则是第 1 轴的下标，

1
2
3
>>> b2 = np.array([[True, False, True], [True, False, False]])
>>> np.nonzero(b2)
(array([0, 0, 1], dtype=int64), array([0, 2, 0], dtype=int64))

庞大的函数库

随机数

函数名 功能 函数名 功能
rand 0到1之间的随机数 randn 标准正态分布的随机数
randint 指定范围内随机整数 normal 正态分布
uniform 均匀分布 poisson 泊松分布
permutation 随机排列 shuffle 随机打乱顺序
choice 随机抽取样本 seed 设置随机数种子