Numpy学习——科学计算
Elementwise操作
基本操作
标量操作：
>>> a = np.array([1, 2, 3, 4])
>>> a + 1
array([2, 3, 4, 5])
>>> 2**a
array([2, 4, 8, 16])
算术操作
>>> b = np.ones(4) + 1
>>> a - b
array([-1., 0., 1., 2.])
>>> a * b
array([2., 4., 6., 8.])
>>> j = np.arange(5)
>>> 2**(j + 1) - j
array([2, 3, 6, 13, 28])
这些操作要比你使用python操作快的多：
>>> a = np.arange(10000)
>>> %timeit a + 1
10000 loops, best of 3: 24.3 us per loop
>>> l = range(10000)
>>> %timeit [i+1 for i in l]
1000 loops, best of 3: 861 us per loop
注意：
数组乘法不等同于矩阵乘法>>> c = np.ones((3, 3)) >>> c * c array([[1., 1., 1.], [1., 1., 1.], [1., 1., 1.]])矩阵乘法结果：
>>> c.dot(c) array([[3., 3., 3.], [3., 3., 3.], [3., 3., 3.]])
其他操作
比较：
>>> a = np.array([1, 2, 3, 4])
>>> b = np.array([4, 2, 2, 4])
>>> a == b
array([False, True, False, True])
>>> a > b
array([False, False, True, False])
数组比较：
>>> a = np.array([1, 2, 3, 4])
>>> b = np.array([4, 2, 2, 4])
>>> c = np.array([1, 2, 3, 4])
>>> np.array_equal(a, b)
False
>>> np.array_equal(a, c)
True
逻辑操作：
>>> a = np.array([1, 1, 0, 0], dtype=bool)
>>> b = np.array([1, 0, 1, 0], dtype=bool)
>>> np.logical_or(a, b)
array([True, True, True, False])
>>> np.logical_and(a, b)
array([True, False, False, False])
超越函数(Transcendental functions)：
>>> a = np.arange(5)
>>> np.sin(a)
array([ 0. , 0.84147098, 0.90929743, 0.14112001, -0.7568025 ])
>>> np.log(a)
array([ -inf, 0. , 0.69314718, 1.09861229, 1.38629436])
>>> np.exp(a)
array([ 1. , 2.71828183, 7.3890561 , 20.08553692, 54.59815003])
维数不匹配：
>>> a = np.arange(4)
>>> a + np.array([1, 2])
Traceback (most recent call last):
File "<pyshell#7>", line 1, in <module>
a + np.array([1, 2])
ValueError: operands could not be broadcast together with shapes (4,) (2,)
我们一会儿再来讨论这个错误
转置：
>>> a = np.triu(np.ones((3, 3)), 1)
>>> a
array([[0., 1., 1.],
[0., 0., 1.],
[0., 0., 0.]])
>>> a.T
array([[0., 0., 0.],
[1., 0., 0.],
[1., 1., 0.]])
triu是返回矩阵的上三角矩阵，后面的参数1代表从左下角开始有多少斜行被置0。
注意：
转置返回的是原始数组的一个视图(view):>>> a = np.arange(9).reshape(3, 3) >>> a.T[0, 2] = 999 >>> a.T array([[ 0, 3, 999], [ 1, 4, 7], [ 2, 5, 8]]) >>> a array([[ 0, 1, 2], [ 3, 4, 5], [999, 7, 8]])
基本约简(reductions)
计算总和
>>> x = np.array([1, 2, 3, 4])
>>> np.sum(x)
10
>>> x.sum()
10
求行列总和：
>>> x = np.array([[1, 1], [2, 2]])
>>> x
array([[1, 1],
[2, 2]])
>>> x.sum(axis=0)
array([3, 3])
>>> x[:, 0].sum(), x[:, 1].sum()
(3, 3)
>>> x.sum(axis=1)
array([2, 4])
>>> x[0, :].sum(), x[1, :].sum()
(2, 4)
设axis=i，则Numpy沿着第i个下标变化的方向进行操作。
相同的思想用在高维数据中：
>>> x = np.random.rand(2, 2, 2)
>>> x
array([[[0.72027886, 0.14531029],
[0.63990608, 0.5921476 ]],
[[0.94029504, 0.36845589],
[0.20379832, 0.75299135]]])
>>> x.sum(axis=2)
array([[0.86558915, 1.23205368],
[1.30875093, 0.95678967]])
>>> x.sum(axis=1)
array([[1.36018494, 0.73745789],
[1.14409336, 1.12144724]])
>>> x.sum(axis=0)
array([[1.6605739 , 0.51376618],
[0.8437044 , 1.34513895]])
>>>
>>> x[0, 1, :].sum()
1.2320536799644604
其他约简
极值(extrema):
>>> x = np.array([1, 3, 2])
>>> x.min()
1
>>> x.max()
3
>>> x.argmin()
0
>>> x.argmax()
1
逻辑操作：
>>> np.all([True, True, False])
False
>>> np.any([True, True, False])
True
注意：可以用于数组比较：
>>> a = np.zeros((100, 100)) >>> np.any(a != 0) False >>> np.any(a != 0) False >>> np.all(a == a) True >>> a == a array([[ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], ..., [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True], [ True, True, True, ..., True, True, True]]) >>> a = np.array([1, 2, 3, 2]) >>> b = np.array([2, 2, 3, 2]) >>> c = np.array([6, 4, 4, 5]) >>> ((a <= b) & (b <= c)).all() True
统计：
>>> x = np.array([1, 2, 3, 1])
>>> y = np.array([[1, 2, 3], [5, 6, 1]])
>>> x.mean()
1.75
>>> np.median(x)
1.5
>>> np.median(y, axis=-1)
array([2., 5.])
>>> x.std()
0.82915619758885
应用举例：
让我们考虑一个简单的一维随机游走的过程，在每一个时刻，一个指针向左和向右跳跃的可能性是相同的。
我们对从指针跳跃t步以后到起点的距离感兴趣。我们将模拟许多指针来找到这个规律，我们将使用数组计算技巧来做到这一点: 我们将创建一个二维数组，其中“经历”(每个指针都有一个经历)在一个方向，时间在另一个方向:
>>> n_experience = 1000
>>> t_max = 200
我们随机选取1和-1，形成以时间为横轴，经历为纵轴的二维数组
>>> t = np.arange(t_max)
>>> steps = 2 * np.random.randint(0, 1 + 1, (n_experience, t_max)) - 1
>>> np.unique(steps)
array([-1, 1])
我们求出累加和按照行的方向：
>>> positions = np.cumsum(steps, axis=1)
>>> sq_distance = positions**2
我们按照列得到经历的平均值
>>> mean_sq_distance = np.mean(sq_distance, axis=0)
画出结果
>>> plt.figure(figsize=(4, 3))
<Figure size 400x300 with 0 Axes>
>>> plt.plot(t, np.sqrt(mean_sq_distance), 'g.', t, np.sqrt(t), 'y-')
[<matplotlib.lines.Line2D object at 0x0000021F445394E0>, <matplotlib.lines.Line2D object at 0x0000021F44539588>]
>>> plt.xlabel(r"$t$")
Text(0.5, 0, '$t$')
>>> plt.ylabel(r"$\sqrt{\langle (\delta x)^2 \rangle}$")
Text(0, 0.5, '$\\sqrt{\\langle (\\delta x)^2 \\rangle}$')
>>> plt,tight_layout()
广播(Broadcasting)
基本的操作都是针对于array数组逐元素的操作，这要求数组要用相同的大小。
此外，如果numpy可以将数组转换成相同的大小，那么这些操作也可以作用在不同大小的数组上，将这种转换称之为广播(Broadcasting)。
下图展现了广播的例子：
让我们证实一下：
>>> a = np.tile(np.arange(0, 40, 10), (3, 1)).T
>>> a
array([[ 0, 0, 0],
[10, 10, 10],
[20, 20, 20],
[30, 30, 30]])
>>> b = np.array([0, 1, 2])
>>> a + b
array([[ 0, 1, 2],
[10, 11, 12],
[20, 21, 22],
[30, 31, 32]])
我们已经尝试过使用广播只是没有注意：
>>> a = np.ones((4, 5))
>>> a[0] = 2
>>> a
array([[2., 2., 2., 2., 2.],
[1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1.]])
一个非常实用的小技巧：
>>> a = np.arange(0, 40, 10)
>>> a.shape
(4,)
>>> a = a[:, np.newaxis]
>>> a.shape
(4, 1)
>>> a
array([[ 0],
[10],
[20],
[30]])
>>> a + b
array([[ 0, 1, 2],
[10, 11, 12],
[20, 21, 22],
[30, 31, 32]])
广播看起来似乎有点神奇，但是实际上在我们处理输出数据维度大于输入数据维度的时候，它很有用处。
应用举例：
让我们构建66号公路城市之间的距离(以英里为单位):芝加哥、斯普林菲尔德、圣路易、塔尔萨、俄克拉荷马城、阿马里洛、圣达菲、阿尔伯克基、弗拉格斯塔夫和洛杉矶。
构造城市之间距离的二维矩阵：
>>> mileposts = np.array([0, 198, 303, 736, 871, 1175, 1475, 1544, 1913, 2448])
>>> distance_array = np.abs(mileposts - mileposts[:, np.newaxis])
>>> distance_array
array([[ 0, 198, 303, 736, 871, 1175, 1475, 1544, 1913, 2448],
[ 198, 0, 105, 538, 673, 977, 1277, 1346, 1715, 2250],
[ 303, 105, 0, 433, 568, 872, 1172, 1241, 1610, 2145],
[ 736, 538, 433, 0, 135, 439, 739, 808, 1177, 1712],
[ 871, 673, 568, 135, 0, 304, 604, 673, 1042, 1577],
[1175, 977, 872, 439, 304, 0, 300, 369, 738, 1273],
[1475, 1277, 1172, 739, 604, 300, 0, 69, 438, 973],
[1544, 1346, 1241, 808, 673, 369, 69, 0, 369, 904],
[1913, 1715, 1610, 1177, 1042, 738, 438, 369, 0, 535],
[2448, 2250, 2145, 1712, 1577, 1273, 973, 904, 535, 0]])
>>> x, y = np.arange(5), np.arange(5)[:, np.newaxis]
>>> distance = np.sqrt(x ** 2 + y ** 2)
>>> distance
array([[0. , 1. , 2. , 3. , 4. ],
[1. , 1.41421356, 2.23606798, 3.16227766, 4.12310563],
[2. , 2.23606798, 2.82842712, 3.60555128, 4.47213595],
[3. , 3.16227766, 3.60555128, 4.24264069, 5. ],
[4. , 4.12310563, 4.47213595, 5. , 5.65685425]])
许多基于网格或者网络的问题都可以使用广播来解决。例如，如果我们要在一个 5 × 5 5 \times 5 5×5的网格当中计算到原点的距离，我们可以用如下的方法：
>>> x, y = np.arange(5), np.arange(5)[:, np.newaxis]
>>> distance = np.sqrt(x ** 2 + y ** 2)
>>> distance
array([[0. , 1. , 2. , 3. , 4. ],
[1. , 1.41421356, 2.23606798, 3.16227766, 4.12310563],
[2. , 2.23606798, 2.82842712, 3.60555128, 4.47213595],
[3. , 3.16227766, 3.60555128, 4.24264069, 5. ],
[4. , 4.12310563, 4.47213595, 5. , 5.65685425]])
>>> plt.pcolor(distance)
>>> plt.colorbar()
>>> plt.show()
注意：numpy.ogrid()函数可以直接生成上述的x，y数组：
>>> x, y = np.ogrid[0:5, 0:5]
>>> x, y
(array([[0],
[1],
[2],
[3],
[4]]), array([[0, 1, 2, 3, 4]]))
>>> x.shape, y.shape
((5, 1), (1, 5))
>>> distance = np.sqrt(x ** 2, y ** 2)
因此，在我们处理表格类问题的时候np.ogrid非常有用。在另一方面，np.mgrid直接提供矩阵的全部索引，当我们不想借助广播的帮助的时候：
>>> x, y = np.mgrid[0:4, 0:4]
>>> x
array([[0, 0, 0, 0],
[1, 1, 1, 1],
[2, 2, 2, 2],
[3, 3, 3, 3]])
>>> y
array([[0, 1, 2, 3],
[0, 1, 2, 3],
[0, 1, 2, 3],
[0, 1, 2, 3]])
数组形状操作
扁平化(Flattening)
>>> a = np.array([[1, 2, 3], [4, 5, 6]])
>>> a.ravel()
array([1, 2, 3, 4, 5, 6])
>>> a.T
array([[1, 4],
[2, 5],
[3, 6]])
>>> a.T.ravel()
array([1, 4, 2, 5, 3, 6])
对于高维数组，最后一个维度优先处理。
重塑(Reshaping)
扁平化的逆处理：
>>> a = np.array([[1, 2, 3], [4, 5, 6]])
>>> a.shape
(2, 3)
>>> b = a.ravel()
>>> b = b.reshape((2, 3))
>>> b
array([[1, 2, 3],
[4, 5, 6]])
或者：
>>> a.reshape((2, -1))
array([[1, 2, 3],
[4, 5, 6]])
>>> b[0, 0] = 99
>>> a
array([[99, 2, 3],
[ 4, 5, 6]])
注意:重塑也可能会返回一个copy!:
>>> a = np.zeros((3, 2))
>>> b = a.T.reshape(3*2)
>>> b[0] = 9
>>> a
array([[0., 0.],
[0., 0.],
[0., 0.]])
为了理解你需要知道更多numpy的存储机制。
增加一个维度
通过使用np.newaxis我们可以使数组增加一个维度：
>>> z = np.array([1, 2, 3])
>>> z
array([1, 2, 3])
>>> z[:, np.newaxis]
array([[1],
[2],
[3]])
>>> z[np.newaxis, :]
array([[1, 2, 3]])
维度混编(Shuffling)
>>> a = np.arange(4*3*2).reshape(4, 3, 2)
>>> a.shape
(4, 3, 2)
>>> a[0, 2, 1]
5
>>> b = a.transpose(1, 2, 0)
>>> b.shape
(3, 2, 4)
>>> b[2, 1, 0]
5
>>> b[2, 1, 0] == a[0, 2, 1]
True
同样也是创建了一个视图：
>>> b[2, 1, 0] = -1
>>> a[0, 2, 1]
-1
更改大小(resizing)
一个数组的大小可以被更改通过ndarray.resize:
>>> a = np.arange(4)
>>> a.resize((8,))
>>> a
array([0, 1, 2, 3, 0, 0, 0, 0])
然而，这要求该数组不能被其他对象所引用：
>>> b = a
>>> a.resize((4,))
Traceback (most recent call last):
File "<pyshell#88>", line 1, in <module>
a.resize((4,))
ValueError: cannot resize an array that references or is referenced
by another array in this way.
Use the np.resize function or refcheck=False
数据排序(Sorting)
沿着一个轴对数组排序：
>>> a = np.array([[4, 3, 5], [1, 2, 1]])
>>> b = np.sort(a, axis=1)
>>> b
array([[3, 4, 5],
[1, 1, 2]])
注意：这是每一行分别排序！
原地排序(In-place sort):
>>> a.sort(axis=1)
>>> a
array([[3, 4, 5],
[1, 1, 2]])
使用花式索引排序：
>>> a = np.array([4, 3, 1, 2])
>>> j = np.argsort(a)
>>> j
array([2, 3, 1, 0], dtype=int64)
>>> a[j]
array([1, 2, 3, 4])
找到最小值和最大值的索引：
>>> a = np.array([4, 3, 1, 2])
>>> j_max = np.argmax(a)
>>> j_min = np.argmin(a)
>>> j_max, j_min
(0, 2)
总结
你需要了解如下知识：
- 知道如何创建一个数组：array, arange, ones, zeros.
- 知道通过array.shape来获取数组的形状，然后使用切割来获取不同的视图: array[::2]。使用reshape或者ravel来改变数组的形状。
- 使用掩膜来获取一个数组的子数组：
>>> a[a < 0] = 0 - 知道数组的多种操作，例如找到平均值或者最大值（array.max(), array.mean())。不需要记住多有东西，但是知道如何在文献中查找。
- 高级使用:掌握整数数组的索引和广播。了解更多NumPy函数来处理各种数组操作。
438
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
