【Python床头书】python NumPy用法示例权威详解

python NumPy用法示例权威详解

学习目标

阅读本文后,您应该能够:

  • 理解NumPy中一维、二维和n维数组之间的区别;
  • 理解如何在不使用for循环的情况下对n维数组应用一些线性代数运算;
  • 理解n维数组的轴和形状属性。

基础知识

NumPy的主要对象是同质的多维数组。它是一个由元素(通常是数字)组成的表格,所有元素具有相同的类型,并由非负整数的元组索引。在NumPy中,维度被称为轴。

例如,一个三维空间中点的坐标数组 [1, 2, 1] 有一个轴。该轴上有3个元素,因此我们说它的长度为3。在下面的示例中,该数组有2个轴。第一个轴的长度为2,第二个轴的长度为3。

[[1., 0., 0.],
 [0., 1., 2.]]

NumPy的数组类称为ndarray。它也被别名array所知。请注意,numpy.array与标准Python库中的array.array不同,后者只处理一维数组并提供较少的功能。ndarray对象的更重要属性包括:

  • ndarray.ndim:数组的轴(维度)数。
  • ndarray.shape:数组的维度。这是一个由整数构成的元组,表示每个维度上数组的大小。对于具有n行和m列的矩阵,形状将是(n,m)。因此,shape元组的长度就是轴的数量,即ndim
  • ndarray.size:数组的总元素数。这等于shape元素的乘积。
  • ndarray.dtype:描述数组中元素类型的对象。可以使用标准Python类型创建或指定dtype。此外,NumPy还提供了自己的类型,例如numpy.int32numpy.int16numpy.float64等。
  • ndarray.itemsize:数组每个元素的字节大小。例如,float64类型的元素数组的itemsize为8(=64/8),而complex32类型的元素数组的itemsize为4(=32/8)。它等同于ndarray.dtype.itemsize
  • ndarray.data:包含数组实际元素的缓冲区。通常情况下,我们不需要使用此属性,因为我们将使用索引访问数组中的元素。

示例

import numpy as np

a = np.arange(15).reshape(3, 5)
print(a)
print("Shape:", a.shape)
print("Dimensions:", a.ndim)
print("Type:", a.dtype.name)
print("Item size:", a.itemsize)
print("Total elements:", a.size)

输出结果如下:

[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]]
Shape: (3, 5)
Dimensions: 2
Type: int64
Item size: 8
Total elements: 15

数组创建

NumPy提供了多种创建数组的方法。例如,您可以使用array函数从常规Python列表或元组创建数组。结果数组的类型将根据序列中元素的类型推断得出。以下是一个示例:

import numpy as np

a = np.array([2, 3, 4])
print(a)
print("Type:", a.dtype)

b = np.array([1.2, 3.5, 5.1])
print(b)
print("Type:", b.dtype)

输出结果如下:

[2 3 4]
Type: int64
[1.2 3.5 5.1]
Type: float64

除了使用序列创建数组外,还可以使用其他函数来创建具有初始占位内容的数组。以下是一些常用的创建数组的函数:

  • zeros:创建一个全零数组。
  • ones:创建一个全为1的数组。
  • empty:创建一个初始内容随机且依赖于内存状态的数组。

这些函数的使用示例如下:

import numpy as np

print(np.zeros((3, 4)))
print(np.ones((2, 3, 4), dtype=np.int16))
print(np.empty((2, 3)))

输出结果如下:

[[0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]]
 
[[[1 1 1 1]
  [1 1 1 1]
  [1 1 1 1]]

 [[1 1 1 1]
  [1 1 1 1]
  [1 1 1 1]]]
 
[[3.73603959e-262 6.02658058e-154 6.55490914e-260]
 [5.30498948e-313 3.14673309e-307 1.00000000e+000]]

NumPy还提供了一些创建数字序列的函数,如arangelinspacearange函数类似于Python内置的range函数,但返回一个数组。linspace函数则用于在指定的范围内生成一定数量的等间隔数值。

以下是这些函数的使用示例:

import numpy as np

print(np.arange(10, 30, 5))
print(np.arange(0, 2, 0.3))

from numpy import pi
print(np.linspace(0, 2, 9))

输出结果如下:

[10 15 20 25]
[0.  0.3 0.6 0.9 1.2 1.5 1.8]
[0.   0.25 0.5  0.75 1.   1.25 1.5  1.75 2.  ]

打印数组

当打印数组时,NumPy会以类似于嵌套列表的方式显示它。一维数组将被打印为行,二维数组将被打印为矩阵,三维数组将被打印为矩阵列表。

以下是打印数组的示例:

import numpy as np

a = np.arange(6)
print(a)

b = np.arange(12).reshape(4, 3)
print(b)

c = np.arange(24).reshape(2, 3, 4)
print(c)

输出结果如下:

[0 1 2 3 4 5]

[[ 0  1  2]
 [ 3  4  5]
 [ 6  7  8]
 [ 9 10 11]]

[[[ 0  1  2  3]
  [ 4  5  6  7]
  [ 8  9 10 11]]

 [[12 13 14 15]
  [16 17 18 19]
  [20 21 22 23]]]

请注意,如果数组太大而无法完全打印,NumPy会自动跳过中心部分,只打印角落部分。您可以使用set_printoptions函数更改打印选项以禁用此行为。

import numpy as np
import sys

np.set_printoptions(threshold=sys.maxsize)

基本操作

数组的算术运算是逐元素进行的。将创建一个新数组,并用结果填充。

a = np.array([20, 30, 40, 50])
b = np.arange(4)
b
array([0, 1, 2, 3])
c = a - b
c
array([20, 29, 38, 47])
b**2
array([0, 1, 4, 9])
10 * np.sin(a)
array([ 9.12945251, -9.88031624,  7.4511316 , -2.62374854])
a < 35
array([ True,  True, False, False])

与许多矩阵语言不同,NumPy数组中的乘法运算符 * 是逐元素进行的。可以使用 @ 运算符(在Python >= 3.5 中)或 dot 函数或方法来执行矩阵乘积:

A = np.array([[1, 1],
              [0, 1]])
B = np.array([[2, 0],
              [3, 4]])
A * B     # 逐元素乘积
array([[2, 0],
       [0, 4]])
A @ B     # 矩阵乘积
array([[5, 4],
       [3, 4]])
A.dot(B)  # 另一种矩阵乘积
array([[5, 4],
       [3, 4]])

某些操作(例如 += 和 *=)是就地修改现有数组,而不是创建新数组。

rg = np.random.default_rng(1)  # 创建默认随机数生成器实例
a = np.ones((2, 3), dtype=int)
b = rg.random((2, 3))
a *= 3
a
array([[3, 3, 3],
       [3, 3, 3]])
b += a
b
array([[3.51182162, 3.9504637 , 3.14415961],
       [3.94864945, 3.31183145, 3.42332645]])
a += b  # b 不会自动转换为整型
Traceback (most recent call last):
    ...
numpy._core._exceptions._UFuncOutputCastingError: Cannot cast ufunc 'add' output from dtype('float64') to dtype('int64') with casting rule 'same_kind'

当对不同类型的数组进行操作时,结果数组的类型与更一般或更精确的类型相对应(这种行为称为向上转型)。

a = np.ones(3, dtype=np.int32)
b = np.linspace(0, pi, 3)
b.dtype.name
'float64'
c = a + b
c
array([1.        , 2.57079633, 4.14159265])
c.dtype.name
'float64'
d = np.exp(c * 1j)
d
array([ 0.54030231+0.84147098j, -0.84147098+0.54030231j,
       -0.54030231-0.84147098j])
d.dtype.name
'complex128'

许多一元操作,例如计算数组中所有元素的和,都作为 ndarray 类的方法实现。

a = rg.random((2, 3))
a
array([[0.82770259, 0.40919914, 0.54959369],
       [0.02755911, 0.75351311, 0.53814331]])
a.sum()
3.1057109529998157
a.min()
0.027559113243068367
a.max()
0.8277025938204418

默认情况下,这些操作将适用于数组,就像它是一个数字列表一样,而不考虑其形状。但是,通过指定 axis 参数,您可以沿着指定的轴应用操作:

b = np.arange(12).reshape(3, 4)
b
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

b.sum(axis=0)     # 每列的和
array([12, 15, 18, 21])

b.min(axis=1)     # 每行的最小值
array([0, 4, 8])

b.cumsum(axis=1)  # 沿每行的累积和
array([[ 0,  1,  3,  6],
       [ 4,  9, 15, 22],
       [ 8, 17, 27, 38]])

通用函数(Universal Functions)

NumPy提供了一些熟悉的数学函数,如sin、cos和exp。在NumPy中,这些被称为“通用函数”(ufunc)。在NumPy中,这些函数对数组逐元素操作,产生一个输出数组。

B = np.arange(3)
B
array([0, 1, 2])
np.exp(B)
array([1.        , 2.71828183, 7.3890561 ])
np.sqrt(B)
array([0.        , 1.        , 1.41421356])
C = np.array([2., -1., 4.])
np.add(B, C)
array([2., 0., 6.])

索引、切片和迭代

一维数组可以像列表和其他Python序列一样进行索引、切片和迭代。

a = np.arange(10)**3
a
array([  0,   1,   8,  27,  64, 125, 216, 343, 512, 729])
a[2]
8
a[2:5]
array([ 8, 27, 64])
# 相当于 a[0:6:2] = 1000;
# 从开始到位置6(不包括),设置每隔2个元素为1000
a[:6:2] = 1000
a
array([1000,    1, 1000,   27, 1000,  125,  216,  343,  512,  729])
a[::-1]  # 倒序排列的a
array([ 729,  512,  343,  216,  125, 1000,   27, 1000,    1, 1000])
for i in a:
    print(i**(1 / 3.))

9.999999999999998  # 可能会有所不同
1.0
9.999999999999998
3.0
9.999999999999998
4.999999999999999
5.999999999999999
6.999999999999999
7.999999999999999
8.999999999999998

多维数组可以每个轴有一个索引。这些索引以逗号分隔的元组给出:

def f(x, y):
    return 10 * x + y

b = np.fromfunction(f, (5, 4), dtype=int)
b
array([[ 0,  1,  2,  3],
       [10, 11, 12, 13],
       [20, 21, 22, 23],
       [30, 31, 32, 33],
       [40, 41, 42, 43]])
b[2, 3]
23
b[0:5, 1]  # b中第二列的每一行
array([ 1, 11, 21, 31, 41])
b[:, 1]    # 与上面的例子等效
array([ 1, 11, 21, 31, 41])
b[1:3, :]  # b中第二和第三行的每一列
array([[10, 11, 12, 13],
       [20, 21, 22, 23]])

当提供的索引少于轴的数量时,缺失的索引被视为完整切片:

b[-1]   # 最后一行。相当于 b[-1, :]
array([40, 41, 42, 43])

在 b[i] 的方括号表达式中,i 被视为一个 i ,后面跟着足够多的 : 来表示剩余的轴。NumPy还允许使用省略号 … 来编写 b[i, …]。

省略号 (…) 表示需要产生一个完整的索引元组所需的冒号数量。例如,如果 x 是一个具有 5 个轴的数组,则

x[1, 2, …] 等效于 x[1, 2, :, :, :]

x[…, 3] 等效于 x[:, :, :, :, 3]

x[4, …, 5, :] 等效于 x[4, :, :, 5, :]

c = np.array([[[  0,  1,  2],  # 一个三维数组(两个堆叠的二维数组)
               [ 10, 12, 13]],
              [[100, 101, 102],
               [110, 112, 113]]])
c.shape
(2, 2, 3)
c[1, ...]  # 等同于 c[1, :, :] 或 c[1]
array([[100, 101, 102],
       [110, 112, 113]])
c[..., 2]  # 等同于 c[:, :, 2]
array([[  2,  13],
       [102, 113]])

对多维数组进行迭代是根据第一个轴进行的:

for row in b:
    print(row)

[0 1 2 3]
[10 11 12 13]
[20 21 22 23]
[30 31 32 33]
[40 41 42 43]

但是,如果要对数组中的每个元素执行操作,可以使用 flat 属性,它是数组的所有元素的迭代器:

for element in b.flat:
    print(element)

0
1
2
3
10
11
12
13
20
21
22
23
30
31
32
33
40
41
42
43

形状操作

改变数组的形状

可以使用各种命令来改变数组的形状。注意,以下三个命令都会返回一个修改后的数组,但不会改变原始数组:

a = np.floor(10 * rg.random((3, 4)))
a
array([[3., 7., 3., 4.],
       [1., 4., 2., 2.],
       [7., 2., 4., 9.]])
a.shape
(3, 4)
# 返回扁平化的数组
a.ravel()
array([3., 7., 3., 4., 1., 4., 2., 2., 7., 2., 4., 9.])
# 返回形状修改后的数组
a.reshape(6, 2)
array([[3., 7.],
       [3., 4.],
       [1., 4.],
       [2., 2.],
       [7., 2.],
       [4., 9.]])
# 返回转置后的数组
a.T
array([[3., 1., 7.],
       [7., 4., 2.],
       [3., 2., 4.],
       [4., 2., 9.]])
a.T.shape
(4, 3)
a.shape
(3, 4)

ravel 方法返回的数组通常以“C风格”存储,也就是说,最右边的索引“变化最快”,所以 a[0, 0] 后面的元素是 a[0, 1]。如果将数组重新形状为其他形状,同样会按照“C风格”处理数组。NumPy通常以这种顺序创建数组,因此ravel通常不需要复制其参数,但如果数组是通过取另一个数组的切片或使用其他选项创建的,则可能需要进行复制。ravelreshape函数还可以使用可选参数指示要使用FORTRAN风格的数组,其中最左边的索引变化得最快。

reshape函数返回带有修改形状的参数,而ndarray.resize方法修改数组本身:

a
array([[3., 7., 3., 4.],
       [1., 4., 2., 2.],
       [7., 2., 4., 9.]])
# 修改数组的形状
a.resize((2, 6))
a
array([[3., 7., 3., 4., 1., 4.],
       [2., 2., 7., 2., 4., 9.]])

如果在重塑操作中给定了 -1 的维度,则会自动计算其他维度:

a.reshape(3, -1)
array([[3., 7., 3., 4.],
       [1., 4., 2., 2.],
       [7., 2., 4., 9.]])

堆叠不同的数组

可以沿不同的轴将多个数组堆叠在一起:

a = np.floor(10 * rg.random((2, 2)))
a
array([[9., 7.],
       [5., 2.]])
b = np.floor(10 * rg.random((2, 2)))
b
array([[1., 9.],
       [5., 1.]])
np.vstack((a, b))  # 沿垂直轴堆叠
array([[9., 7.],
       [5., 2.],
       [1., 9.],
       [5., 1.]])
np.hstack((a, b))  # 沿水平轴堆叠
array([[9., 7., 1., 9.],
       [5., 2., 5., 1.]])

column_stack函数将一维数组作为列堆叠成二维数组。对于二维数组,它等同于hstack

np.column_stack((a, b))  # 对于二维数组
array([[9., 7., 1., 9.],
       [5., 2., 5., 1.]])
a = np.array([4., 2.])
b = np.array([3., 8.])
np.column_stack((a, b))  # 返回一个二维数组
array([[4., 3.],
       [2., 8.]])
np.hstack((a, b))        # 结果不同
array([4., 2., 3., 8.])
a[:, np.newaxis]  # 将 `a` 视为二维列向量
array([[4.],
       [2.]])
np.column_stack((a[:, np.newaxis], b[:, np.newaxis]))
array([[4., 3.],
       [2., 8.]])
np.hstack((a[:, np.newaxis], b[:, np.newaxis]))  # 结果相同
array([[4., 3.],
       [2., 8.]])

对于多维数组,堆叠函数的行为取决于沿着哪个轴进行堆叠。hstack沿第二个轴堆叠,vstack沿第一个轴堆叠,concatenate允许指定要进行连接的轴的可选参数。

将一个数组分割成多个较小的数组

使用hsplit可以沿水平轴将一个数组分割成多个子数组,可以通过指定要返回的等形状的数组数量或指定分割点之后的列来完成:

a = np.floor(10 * rg.random((2, 12)))
a
array([[6., 7., 6., 9., 0., 5., 4., 0., 6., 8., 5., 2.],
       [8., 5., 5., 7., 1., 8., 6., 7., 1., 8., 1., 0.]])
# 将 `a` 分割成 3 个数组
np.hsplit(a, 3)
[array([[6., 7., 6., 9.],
       [8., 5., 5., 7.]]), array([[0., 5., 4., 0.],
       [1., 8., 6., 7.]]), array([[6., 8., 5., 2.],
       [1., 8., 1., 0.]])]
# 在第三列和第四列之后将 `a` 分割
np.hsplit(a, (3, 4))
[array([[6., 7., 6.],
       [8., 5., 5.]]), array([[9.],
       [7.]]), array([[0., 5., 4., 0., 6., 8., 5., 2.],
       [1., 8., 6., 7., 1., 8., 1., 0.]])]

vsplit沿垂直轴分割,array_split允许指定要进行分割的轴。

复制和视图

在操作和处理数组时,有时会将数据复制到新的数组中,有时则不会。这常常会使初学者感到困惑。有三种情况:

完全没有复制

简单的赋值不会复制对象或其数据。

a = np.array([[ 0,  1,  2,  3],
              [ 4,  5,  6,  7],
              [ 8,  9, 10, 11]])
b = a            # 没有创建新的对象
b is a           # a 和 b 是同一个 ndarray 对象的两个名称
True

Python将可变对象作为引用传递,因此函数调用不会复制对象。

def f(x):
    print(id(x))

id(a)  # id 是对象的唯一标识符
148293216  # 可能不同
f(a)
148293216  # 可能不同

视图或浅复制

不同的数组对象可以共享相同的数据。view方法创建一个新的数组对象,它查看相同的数据。

c = a.view()
c is a
False
c.base is a            # c 是 a 拥有的数据的视图
True
c.flags.owndata
False

c = c.reshape((2, 6))  # a 的形状不会改变
a.shape
(3, 4)
c[0, 4] = 1234         # a 的数据发生了改变
a
array([[   0,    1,    2,    3],
       [1234,    5,    6,    7],
       [   8,    9,   10,   11]])

对数组进行切片会返回它的视图:

s = a[:, 1:3]
s[:] = 10  # s[:] 是 s 的视图。注意 s = 10 和 s[:] = 10 的区别
a
array([[   0,   10,   10,    3],
       [1234,   10,   10,    7],
       [   8,   10,   10,   11]])

深复制

copy方法创建一个数组及其数据的完整副本。

d = a.copy()  # 创建一个具有新数据的新数组对象
d is a
False
d.base is a  # d 不与 a 共享任何内容
False
d[0, 0] = 9999
a
array([[   0,   10,   10,    3],
       [1234,   10,   10,    7],
       [   8,   10,   10,   11]])

如果切片后需要复制原始数组,则应在切片之后使用copy。例如,假设 a 是一个巨大的中间结果,最终结果 b 只包含 a 的一小部分,当使用切片构建 b 时,应进行深复制:

a = np.arange(int(1e8))
b = a[:100].copy()
del a  # 可以释放 `a` 的内存。

如果使用b = a[:100],则 a 被 b 引用,并且即使执行del a,它仍然会保留在内存中。

函数和方法概述

下面是一些有用的NumPy函数和方法名称的列表,按类别排序。有关完整列表,请参阅相关文档。

数组创建
arange, array, copy, empty, empty_like, eye, fromfile, fromfunction, identity, linspace, logspace, mgrid, ogrid, ones, ones_like, r_, zeros, zeros_like

转换
ndarray.astype, atleast_1d, atleast_2d, atleast_3d, mat

操作
array_split, column_stack, concatenate, diagonal, dsplit, dstack, hsplit, hstack, ndarray.item, newaxis, ravel, repeat, reshape, resize, squeeze, swapaxes, take, transpose, vsplit, vstack

问题
all, any, nonzero, where

排序
argmax, argmin, argsort, max, min, ptp, searchsorted, sort

运算
choose, compress, cumprod, cumsum, inner, ndarray.fill, imag, prod, put, putmask, real, sum

基本统计
cov, mean, std, var

基本线性代数
cross, dot, outer, linalg.svd, vdot

更高级的功能

广播规则

广播允许通用函数以有意义的方式处理形状不同的输入。广播的第一条规则是,如果所有输入数组的维度数目不相同,则“1”将重复添加到较小数组的形状的最前面,直到所有数组具有相同的维度数目。

广播的第二条规则确保在特定维度上尺寸为 1 的数组的行为就像它们具有该尺寸的最大形状的数组一样。假设数组 a 具有大小为 1 的特定维度,那么该数组的元素值在该维度上被认为是相同的,对于“广播”数组来说。

应用广播规则后,所有数组的大小必须匹配。更多详细信息请参见广播。

高级索引和索引技巧

NumPy提供了比常规Python序列更多的索引功能。除了整数和切片索引之外,还可以使用整数数组和布尔数组进行索引。

使用索引数组进行索引

a = np.arange(12)**2  # 前12个平方数
i = np.array([1, 1, 3, 8, 5])  # 索引数组
a[i]  # 在位置 `i` 处的元素
array([ 1,  1,  9, 64, 25])

j = np.array([[3, 4], [9, 7]])  # 二维索引数组
a[j]  # 形状与 `j` 相同
array([[ 9, 16],
       [81, 49]])

当被索引的数组 a 是多维的时候,单个索引数组指的是 a 的第一个维度。以下示例通过使用调色板将标签图像转换为彩色图像来展示此行为。

palette = np.array([[0, 0, 0],         # 黑色
                    [255, 0, 0],       # 红色
                    [0, 255, 0],       # 绿色
                    [0, 0, 255],       # 蓝色
                    [255, 255, 255]])  # 白色
image = np.array([[0, 1, 2, 0],  # 每个值对应调色板中的颜色
                  [0, 3, 4, 0]])
palette[image]  # 形状为 (2, 4, 3) 的彩色图像
array([[[  0,   0,   0],
        [255,   0,   0],
        [  0, 255,   0],
        [  0,   0,   0]],

       [[  0,   0,   0],
        [  0,   0, 255],
        [255, 255, 255],
        [  0,   0,   0]]])

可以给出每个维度的多个索引数组。各个维度的索引数组必须具有相同的形状。

a = np.arange(12).reshape(3, 4)
a
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
i = np.array([[0, 1],  # `a` 的第一维的索引
              [1, 2]])
j = np.array([[2, 1],  # 第二个维度的索引
              [3, 3]])

a[i, j]  # `i` 和 `j` 必须具有相同的形状
array([[ 2,  5],
       [ 7, 11]])

a[i, 2]
array([[ 2,  6],
       [ 6, 10]])

a[:, j]
array([[[ 2,  1],
        [ 3,  3]],

       [[ 6,  5],
        [ 7,  7]],

       [[10,  9],
        [11, 11]]])

在 Python 中,arr[i, j]arr[(i, j)] 是完全相同的,所以我们可以将 ij 放入一个元组中,然后使用该元组进行索引。

l = (i, j)
# 与 a[i, j] 等效
a[l]
array([[ 2,  5],
       [ 7, 11]])

但是,如果将 i 和 j 放入一个数组中,则无法进行此操作,因为该数组将被解释为对 a 的第一个维度的索引。

s = np.array([i, j])
# 不是我们想要的结果
a[s]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
IndexError: index 3 is out of bounds for axis 0 with size 3
# 与 a[i, j] 相同
a[tuple(s)]
array([[ 2,  5],
       [ 7, 11]])

使用布尔数组进行索引

当我们使用整数索引数组来索引数组时,我们提供了要选择的索引列表。而使用布尔索引的方法不同;我们明确选择数组中我们想要的项目和不想要的项目。

最自然的布尔索引方式是使用与原始数组具有相同形状的布尔数组:

a = np.arange(12).reshape(3, 4)
b = a > 4
b  # `b` 是一个布尔数组,与 `a` 相同的形状
array([[False, False, False, False],
       [False,  True,  True,  True],
       [ True,  True,  True,  True]])
a[b]  # 选取的元素为一维数组
array([ 5,  6,  7,  8,  9, 10, 11])

这个特性在赋值时非常有用:

a[b] = 0  # 所有大于4的元素变为0
a
array([[0, 1, 2, 3],
       [4, 0, 0, 0],
       [0, 0, 0, 0]])

你可以查看下面的示例,了解如何使用布尔索引生成Mandelbrot集合的图像:

import numpy as np
import matplotlib.pyplot as plt
def mandelbrot(h, w, maxit=20, r=2):
    """返回大小为(h,w)的Mandelbrot分形图像。"""
    x = np.linspace(-2.5, 1.5, 4*h+1)
    y = np.linspace(-1.5, 1.5, 3*w+1)
    A, B = np.meshgrid(x, y)
    C = A + B*1j
    z = np.zeros_like(C)
    divtime = maxit + np.zeros(z.shape, dtype=int)

    for i in range(maxit):
        z = z**2 + C
        diverge = abs(z) > r                    # 谁发散了
        div_now = diverge & (divtime == maxit)  # 现在谁在发散
        divtime[div_now] = i                    # 记录时间
        z[diverge] = r                          # 避免过度发散

    return divtime
plt.clf()
plt.imshow(mandelbrot(400, 400))
../_images/quickstart-1.png

使用布尔索引的第二种方式更类似于整数索引;对于数组的每个维度,我们给出一个1D布尔数组来选择我们想要的切片:

a = np.arange(12).reshape(3, 4)
b1 = np.array([False, True, True])         # 第一维度的选择
b2 = np.array([True, False, True, False])  # 第二维度的选择

a[b1, :]                                   # 选择行
array([[ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

a[b1]                                      # 同样的结果
array([[ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

a[:, b2]                                   # 选择列
array([[ 0,  2],
       [ 4,  6],
       [ 8, 10]])

a[b1, b2]                                  # 一个奇怪的操作
array([ 4, 10])

请注意,1D布尔数组的长度必须与您要切片的维度(或轴)的长度相同。在上面的示例中,b1的长度为3(a的行数),b2的长度为4适用于索引a的第二个轴(列)。

ix_() 函数

ix_()函数可用于组合不同的向量,以获得每个n元组的结果。例如,如果要计算从向量a、b和c中取出的所有三元组的 a+b*c:

a = np.array([2, 3, 4, 5])
b = np.array([8, 5, 4])
c = np.array([5, 4, 6, 8, 3])
ax, bx, cx = np.ix_(a, b, c)
ax
array([[[2]],

       [[3]],

       [[4]],

       [[5]]])
bx
array([[[8],
        [5],
        [4]]])
cx
array([[[5, 4, 6, 8, 3]]])
ax.shape, bx.shape, cx.shape
((4, 1, 1), (1, 3, 1), (1, 1, 5))
result = ax + bx * cx
result
array([[[42, 34, 50, 66, 26],
        [27, 22, 32, 42, 17],
        [22, 18, 26, 34, 14]],

       [[43, 35, 51, 67, 27],
        [28, 23, 33, 43, 18],
        [23, 19, 27, 35, 15]],

       [[44, 36, 52, 68, 28],
        [29, 24, 34, 44, 19],
        [24, 20, 28, 36, 16]],

       [[45, 37, 53, 69, 29],
        [30, 25, 35, 45, 20],
        [25, 21, 29, 37, 17]]])
result[3, 2, 4]
17
a[3] + b[2] * c[4]
17

你还可以按以下方式实现reduce:

def ufunc_reduce(ufct, *vectors):
   vs = np.ix_(*vectors)
   r = ufct.identity
   for v in vs:
       r = ufct(r, v)
   return r
and then use it as:

ufunc_reduce(np.add, a, b, c)
array([[[15, 14, 16, 18, 13],
        [12, 11, 13, 15, 10],
        [11, 10, 12, 14,  9]],

       [[16, 15, 17, 19, 14],
        [13, 12, 14, 16, 11],
        [12, 11, 13, 15, 10]],

       [[17, 16, 18, 20, 15],
        [14, 13, 15, 17, 12],
        [13, 12, 14, 16, 11]],

       [[18, 17, 19, 21, 16],
        [15, 14, 16, 18, 13],
        [14, 13, 15, 17, 12]]])

与普通的ufunc.reduce相比,这个版本的reduce的优势在于它利用了广播规则,以避免创建一个大小为输出大小乘以向量数量的参数数组。

使用字符串进行索引

参见Structured arrays.

技巧和提示

以下是一些简短而有用的技巧。

自动重塑

要更改数组的维度,您可以省略其中一个尺寸,它将自动推断出来:

a = np.arange(30)
b = a.reshape((2, -1, 3))  # -1 表示“任意大小”
b.shape
(2, 5, 3)
b
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]]])

向量堆叠

如何从一组相同长度的行向量构建二维数组?在MATLAB中,这很容易:如果x和y是两个长度相同的向量,则只需执行m=[x;y]。在NumPy中,可以使用column_stack、dstack、hstack和vstack函数,具体取决于要进行堆叠的维度。例如:

x = np.arange(0, 10, 2)
y = np.arange(5)
m = np.vstack([x, y])
m
array([[0, 2, 4, 6, 8],
       [0, 1, 2, 3, 4]])
xy = np.hstack([x, y])
xy
array([0, 2, 4, 6, 8, 0, 1, 2, 3, 4])

在多于两个维度的情况下,这些函数的逻辑可能会有些奇怪。

直方图

将NumPy的histogram函数应用于数组会返回一对向量:数组的直方图和一个边界的向量。注意:matplotlib还有一个用于构建直方图的函数(称为hist,类似于Matlab),它与NumPy中的函数不同。主要区别在于pylab.hist会自动绘制直方图,而numpy.histogram只生成数据。

import numpy as np
rg = np.random.default_rng(1)
import matplotlib.pyplot as plt
# 使用方差为0.5^2和均值为2的正态分布生成10000个随机数
mu, sigma = 2, 0.5
v = rg.normal(mu, sigma, 10000)
# 使用50个柱状图绘制归一化直方图
plt.hist(v, bins=50, density=True)       # matplotlib 版本(绘图)
(array...)
# 使用NumPy计算直方图,然后绘制它
(n, bins) = np.histogram(v, bins=50, density=True)  # NumPy 版本(不绘图)
plt.plot(.5 * (bins[1:] + bins[:-1]), n) 
../_images/quickstart-2.png

在Matplotlib >=3.4版本中,还可以使用plt.stairs(n, bins)绘制直方图。

参考链接

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

BigDataMLApplication

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值