In [1]:
import numpy as np
a = np.array([0, 1, 2, 3])
a
Out[1]:
例如，数组包含：
- 实验或模拟在离散时间阶段的值
- 测量设备记录的信号，比如声波
- 图像的像素、灰度或颜色
- 用不同X-Y-Z位置测量的3-D数据，例如MRI扫描
...
为什么有用：提供了高速数值操作的节省内存的容器。
In [2]:
L = range(1000)
%timeit [i**2 for i in L]
In [4]:
a = np.arange(1000)
%timeit a**2
1.3.1.1.2 Numpy参考文档¶
- 线上: http://docs.scipy.org/
- 交互帮助:
np.array?
String Form:<built-in function array>
Docstring:
array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0, ...
查找东西：
In [6]:
np.lookfor('create array')
In [8]:
import numpy as np
In [9]:
a = np.array([0, 1, 2, 3])
a
Out[9]:
In [10]:
a.ndim
Out[10]:
In [11]:
a.shape
Out[11]:
In [12]:
len(a)
Out[12]:
- 2-D，3-D，...：
In [13]:
b = np.array([[0, 1, 2], [3, 4, 5]]) # 2 x 3 数组
b
Out[13]:
In [14]:
b.ndim
Out[14]:
In [15]:
b.shape
Out[15]:
In [16]:
len(b) # 返回一个纬度的大小
Out[16]:
In [17]:
c = np.array([[[1], [2]], [[3], [4]]])
c
Out[17]:
In [18]:
c.shape
Out[18]:
In [19]:
a = np.arange(10) # 0 .. n-1 (!)
a
Out[19]:
In [20]:
b = np.arange(1, 9, 2) # 开始，结束（不包含），步长
b
Out[20]:
- 或者通过一些数据点：
In [1]:
c = np.linspace(0, 1, 6) # 起点、终点、数据点
c
Out[1]:
In [2]:
d = np.linspace(0, 1, 5, endpoint=False)
d
Out[2]:
- 普通数组：
In [3]:
a = np.ones((3, 3)) # 提示: (3, 3) 是元组
a
Out[3]:
In [4]:
b = np.zeros((2, 2))
b
Out[4]:
In [5]:
c = np.eye(3)
c
Out[5]:
In [6]:
d = np.diag(np.array([1, 2, 3, 4]))
d
Out[6]:
np.random: 随机数 (Mersenne Twister PRNG) :
In [7]:
a = np.random.rand(4) # [0, 1] 的均匀分布
a
Out[7]:
In [8]:
b = np.random.randn(4) # 高斯
b
Out[8]:
In [9]:
np.random.seed(1234) # 设置随机种子
In [10]:
np.random.rand?
练习：用函数创建数组
- 实验用
arange、linspace、ones、zeros、eye和diag。 - 用随机数创建不同类型的数组。
- 在创建带有随机数的数组前设定种子。
- 看一下函数
np.empty。它能做什么？什么时候会比较有用？
1.3.1.3基础数据类型¶
你可能已经发现，在一些情况下，数组元素显示带有点（即 2. VS 2）。这是因为所使用的数据类型不同：
In [12]:
a = np.array([1, 2, 3])
a.dtype
Out[12]:
In [13]:
b = np.array([1., 2., 3.])
b
Out[13]:
不同的数据类型可以更紧凑的在内存中存储数据，但是大多数时候我们都只是操作浮点数据。注意，在上面的例子中，Numpy自动从输入中识别了数据类型。
你可以明确的指定想要的类型：
In [1]:
c = np.array([1, 2, 3], dtype=float)
c.dtype
Out[1]:
默认数据类型是浮点:
In [2]:
a = np.ones((3, 3))
a.dtype
Out[2]:
其他类型：
复数：
In [4]:
d = np.array([1+2j, 3+4j, 5+6*1j])
d.dtype
Out[4]:
布尔：
In [5]:
e = np.array([True, False, False, True])
e.dtype
Out[5]:
字符：
In [6]:
f = np.array(['Bonjour', 'Hello', 'Hallo',])
f.dtype # <--- 包含最多7个字母的字符
Out[6]:
更多：
- int32
- int64
- unit32
- unit64
1.3.1.4基本可视化¶
现在我们有了第一个数组，我们将要进行可视化。
从pylab模式启动IPython。
ipython --pylab或notebook：
ipython notebook --pylab=inline或者如果IPython已经启动，那么：
In [119]:
%pylab
或者从Notebook中：
In [121]:
%pylab inline
inline 对notebook来说很重要，以便绘制的图片在notebook中显示而不是在新窗口显示。
Matplotlib是2D制图包。我们可以像下面这样导入它的方法：
In [10]:
import matplotlib.pyplot as plt #整洁形式
然后使用（注你需要显式的使用 show ）:
plt.plot(x, y) # 线图
plt.show() # <-- 显示图表（使用pylab的话不需要）
或者，如果你使用 pylab：
plt.plot(x, y) # 线图
在脚本中推荐使用 import matplotlib.pyplot as plt。 而交互的探索性工作中用 pylab。
- 1D作图：
In [12]:
x = np.linspace(0, 3, 20)
y = np.linspace(0, 9, 20)
plt.plot(x, y) # 线图
Out[12]:
In [13]:
plt.plot(x, y, 'o') # 点图
Out[13]:
- 2D 作图:
In [14]:
image = np.random.rand(30, 30)
plt.imshow(image, cmap=plt.cm.hot)
plt.colorbar()
Out[14]:
更多请见matplotlib部分
练习：简单可视化
画出简单的数组：cosine作为时间的一个函数以及2D矩阵。
在2D矩阵上试试使用 gray colormap。
1.3.1.5索引和切片¶
数组的项目可以用与其他Python序列（比如：列表）一样的方式访问和赋值：
In [15]:
a = np.arange(10)
a
Out[15]:
In [16]:
a[0], a[2], a[-1]
Out[16]:
警告：索引从0开始与其他的Python序列（以及C/C++）一样。相反，在Fortran或者Matlab索引从1开始。
使用常用的Python风格来反转一个序列也是支持的：
In [17]:
a[::-1]
Out[17]:
对于多维数组，索引是整数的元组：
In [18]:
a = np.diag(np.arange(3))
a
Out[18]:
In [19]:
a[1, 1]
Out[19]:
In [21]:
a[2, 1] = 10 # 第三行，第二列
a
Out[21]:
In [22]:
a[1]
Out[22]:
注：
- 在2D数组中，第一个维度对应行，第二个维度对应列。
- 对于多维度数组
a，a[0]被解释为提取在指定维度的所有元素
切片：数组与其他Python序列也可以被切片：
In [23]:
a = np.arange(10)
a
Out[23]:
In [24]:
a[2:9:3] # [开始:结束:步长]
Out[24]:
注意最后一个索引是不包含的！：
In [25]:
a[:4]
Out[25]:
切片的三个元素都不是必选：默认情况下，起点是0，结束是最后一个，步长是1：
In [26]:
a[1:3]
Out[26]:
In [27]:
a[::2]
Out[27]:
In [28]:
a[3:]
Out[28]:
Numpy索引和切片的一个小说明...
赋值和切片可以结合在一起：
In [29]:
a = np.arange(10)
a[5:] = 10
a
Out[29]:
In [30]:
b = np.arange(5)
a[5:] = b[::-1]
a
Out[30]:
练习：索引与切片
- 试试切片的特色，用起点、结束和步长：从linspace开始，试着从后往前获得奇数，从前往后获得偶数。 重现上面示例中的切片。你需要使用下列表达式创建这个数组：
In [31]:
np.arange(6) + np.arange(0, 51, 10)[:, np.newaxis]
Out[31]:
练习：数组创建
创建下列的数组（用正确的数据类型）：
[[1, 1, 1, 1],
[1, 1, 1, 1],
[1, 1, 1, 2],
[1, 6, 1, 1]]
[[0., 0., 0., 0., 0.],
[2., 0., 0., 0., 0.],
[0., 3., 0., 0., 0.],
[0., 0., 4., 0., 0.],
[0., 0., 0., 5., 0.],
[0., 0., 0., 0., 6.]]
参考标准：每个数组
提示：每个数组元素可以像列表一样访问，即a[1] 或 a[1, 2]。
提示：看一下 diag 的文档字符串。
练习：创建平铺数组
看一下 np.tile 的文档，是用这个函数创建这个数组：
[[4, 3, 4, 3, 4, 3],
[2, 1, 2, 1, 2, 1],
[4, 3, 4, 3, 4, 3],
[2, 1, 2, 1, 2, 1]]
1.3.1.6 副本和视图¶
切片操作创建原数组的一个视图，这只是访问数组数据一种方式。因此，原始的数组并不是在内存中复制。你可以用 np.may_share_memory() 来确认两个数组是否共享相同的内存块。但是请注意，这种方式使用启发式，可能产生漏报。
当修改视图时，原始数据也被修改：
In [32]:
a = np.arange(10)
a
Out[32]:
In [33]:
b = a[::2]
b
Out[33]:
In [34]:
np.may_share_memory(a, b)
Out[34]:
In [36]:
b[0] = 12
b
Out[36]:
In [37]:
a # (!)
Out[37]:
In [38]:
a = np.arange(10)
c = a[::2].copy() # 强制复制
c[0] = 12
a
Out[38]:
In [39]:
np.may_share_memory(a, c)
Out[39]:
乍看之下这种行为可能有些奇怪，但是这样做节省了内存和时间。
实例：素数筛选
用筛选法计算0-99之间的素数
- 构建一个名为
_prime形状是 (100,) 的布尔数组，在初始将值都设为True：
In [40]:
is_prime = np.ones((100,), dtype=bool)
- 将不属于素数的0，1去掉
In [41]:
is_prime[:2] = 0
对于从2开始的整数 j ，化掉它的倍数：
In [42]:
N_max = int(np.sqrt(len(is_prime)))
for j in range(2, N_max):
is_prime[2*j::j] = False
- 看一眼
help(np.nonzero)，然后打印素数 - 接下来:
- 将上面的代码放入名为
prime_sieve.py的脚本文件 - 运行检查一下时候有效
- 使用埃拉托斯特尼筛法的优化建议
- 跳过已知不是素数的
j - 第一个应该被划掉的数是$j^2$
- 跳过已知不是素数的
- 将上面的代码放入名为
In [44]:
np.random.seed(3)
a = np.random.random_integers(0, 20, 15)
a
Out[44]:
In [45]:
(a % 3 == 0)
Out[45]:
In [47]:
mask = (a % 3 == 0)
extract_from_a = a[mask] # 或, a[a%3==0]
extract_from_a # 用面具抽取一个子数组
Out[47]:
赋值给子数组时，用面具索引非常有用：
In [48]:
a[a % 3 == 0] = -1
a
Out[48]:
1.3.1.7.2 用整型数组索引¶
In [49]:
a = np.arange(0, 100, 10)
a
Out[49]:
索引可以用整型数组完成，其中相同的索引重复了几次：
In [50]:
a[[2, 3, 2, 4, 2]] # 注：[2, 3, 2, 4, 2] 是Python列表
Out[50]:
用这种类型的索引可以分配新值：
In [51]:
a[[9, 7]] = -100
a
Out[51]:
当一个新数组用整型数组索引创建时，新数组有相同的形状，而不是整数数组：
In [52]:
a = np.arange(10)
idx = np.array([[3, 4], [9, 7]])
idx.shape
Out[52]:
In [53]:
a[idx]
Out[53]:
下图展示了多种象征索引的应用
练习：象征索引
- 同样，重新生成上图中所示的象征索引
- 用左侧的象征索引和右侧的数组创建在为一个数组赋值，例如，设置上图数组的一部分为0。
In [54]:
a = np.array([1, 2, 3, 4])
a + 1
Out[54]:
In [55]:
2**a
Out[55]:
所有运算是在元素级别上操作：
In [56]:
b = np.ones(4) + 1
a - b
Out[56]:
In [57]:
a * b
Out[57]:
In [58]:
j = np.arange(5)
2**(j + 1) - j
Out[58]:
这些操作当然也比你用纯Python实现好快得多：
In [60]:
a = np.arange(10000)
%timeit a + 1
In [61]:
l = range(10000)
%timeit [i+1 for i in l]
注意：数组相乘不是矩阵相乘：
In [62]:
c = np.ones((3, 3))
c * c # 不是矩阵相乘！
Out[62]:
注：矩阵相乘：
In [63]:
c.dot(c)
Out[63]:
练习：元素级别的操作
- 试一下元素级别的简单算术操作
- 用
%timeit比一下他们与纯Python对等物的时间 - 生成：
[2**0, 2**1, 2**2, 2**3, 2**4]a_j = 2^(3*j) - j
1.3.2.1.2其他操作¶
对比：
In [64]:
a = np.array([1, 2, 3, 4])
b = np.array([4, 2, 2, 4])
a == b
Out[64]:
In [65]:
a > b
Out[65]:
数组级别的对比：
In [66]:
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)
Out[66]:
In [67]:
np.array_equal(a, c)
Out[67]:
逻辑操作：
In [68]:
a = np.array([1, 1, 0, 0], dtype=bool)
b = np.array([1, 0, 1, 0], dtype=bool)
np.logical_or(a, b)
Out[68]:
In [69]:
np.logical_and(a, b)
Out[69]:
In [71]:
a = np.arange(5)
np.sin(a)
Out[71]:
In [72]:
np.log(a)
Out[72]:
In [73]:
np.exp(a)
Out[73]:
形状不匹配
In [74]:
a = np.arange(4)
a + np.array([1, 2])
广播？我们将在稍后讨论。
变换
In [76]:
a = np.triu(np.ones((3, 3)), 1) # 看一下 help(np.triu)
a
Out[76]:
In [77]:
a.T
Out[77]:
警告：变换是视图
因此，下列的代码是错误的，将导致矩阵不对称：
In [78]:
a += a.T
In [79]:
x = np.array([1, 2, 3, 4])
np.sum(x)
Out[79]:
In [80]:
x.sum()
Out[80]:
行求和和列求和：
In [81]:
x = np.array([[1, 1], [2, 2]])
x
Out[81]:
In [83]:
x.sum(axis=0) # 列 (第一纬度)
Out[83]:
In [84]:
x[:, 0].sum(), x[:, 1].sum()
Out[84]:
In [85]:
x.sum(axis=1) # 行 (第二纬度)
Out[85]:
In [86]:
x[0, :].sum(), x[1, :].sum()
Out[86]:
高维的处理，思路相同：
In [87]:
x = np.random.rand(2, 2, 2)
x.sum(axis=2)[0, 1]
Out[87]:
In [88]:
x[0, 1, :].sum()
Out[88]:
1.3.2.2.2 其他简化¶
- 以相同方式运作（也可以使用
axis=）
极值
In [89]:
x = np.array([1, 3, 2])
x.min()
Out[89]:
In [90]:
x.max()
Out[90]:
In [91]:
x.argmin() # 最小值的索引
Out[91]:
In [92]:
x.argmax() # 最大值的索引
Out[92]:
逻辑运算：
In [93]:
np.all([True, True, False])
Out[93]:
In [94]:
np.any([True, True, False])
Out[94]:
注：可以被应用数组对比：
In [95]:
a = np.zeros((100, 100))
np.any(a != 0)
Out[95]:
In [96]:
np.all(a == a)
Out[96]:
In [97]:
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()
Out[97]:
统计:
In [98]:
x = np.array([1, 2, 3, 1])
y = np.array([[1, 2, 3], [5, 6, 1]])
x.mean()
Out[98]:
In [99]:
np.median(x)
Out[99]:
In [100]:
np.median(y, axis=-1) # 最后的坐标轴
Out[100]:
In [101]:
x.std() # 全体总体的标准差。
Out[101]:
... 以及其他更多（随着你成长最好学习一下）。
练习：简化
- 假定有
sum，你会期望看到哪些其他的函数？ sum和cumsum有什么区别？
实例: 数据统计
populations.txt中的数据描述了过去20年加拿大北部野兔和猞猁的数量（以及胡萝卜）。
你可以在编辑器或在IPython看一下数据（shell或者notebook都可以）：
In [104]:
cat data/populations.txt
首先，将数据加载到Numpy数组：
In [107]:
data = np.loadtxt('data/populations.txt')
year, hares, lynxes, carrots = data.T # 技巧: 将列分配给变量
接下来作图：
In [108]:
from matplotlib import pyplot as plt
plt.axes([0.2, 0.1, 0.5, 0.8])
plt.plot(year, hares, year, lynxes, year, carrots)
plt.legend(('Hare', 'Lynx', 'Carrot'), loc=(1.05, 0.5))
Out[108]:
随时间变化的数量的平均数：
In [109]:
populations = data[:, 1:]
populations.mean(axis=0)
Out[109]:
样本的标准差：
In [110]:
populations.std(axis=0)
Out[110]:
每一年哪个物种有最高的数量？：
In [111]:
np.argmax(populations, axis=1)
Out[111]:
实例：随机游走算法扩散
让我们考虑一下简单的1维随机游走过程：在每个时间点，行走者以相等的可能性跳到左边或右边。我们感兴趣的是找到随机游走者在 t 次左跳或右跳后距离原点的典型距离？我们将模拟许多”行走者“来找到这个规律，并且我们将采用数组计算技巧来计算：我们将创建一个2D数组记录事实，一个方向是经历（每个行走者有一个经历），一个纬度是时间：
In [113]:
n_stories = 1000 # 行走者的数
t_max = 200 # 我们跟踪行走者的时间
我们随机选择步长1或-1去行走：
In [115]:
t = np.arange(t_max)
steps = 2 * np.random.random_integers(0, 1, (n_stories, t_max)) - 1
np.unique(steps) # 验证: 所有步长是1或-1
Out[115]:
我们通过汇总随着时间的步骤来构建游走
In [116]:
positions = np.cumsum(steps, axis=1) # axis = 1: 纬度是时间
sq_distance = positions**2
获得经历轴的平均数：
In [117]:
mean_sq_distance = np.mean(sq_distance, axis=0)
画出结果：
In [126]:
plt.figure(figsize=(4, 3))
plt.plot(t, np.sqrt(mean_sq_distance), 'g.', t, np.sqrt(t), 'y-')
plt.xlabel(r"$t$")
plt.ylabel(r"$\sqrt{\langle (\delta x)^2 \rangle}$")
Out[126]:
我们找到了物理学上一个著名的结果：均方差记录是时间的平方根！
1.3.2.3 广播¶
- numpy数组的基本操作（相加等）是元素级别的
- 在相同大小的数组上仍然适用。 尽管如此, 也可能在不同大小的数组上进行这个操作，假如Numpy可以将这些数组转化为相同的大小：这种转化称为广播。
下图给出了一个广播的例子：
让我们验证一下：
In [127]:
a = np.tile(np.arange(0, 40, 10), (3, 1)).T
a
Out[127]:
In [128]:
b = np.array([0, 1, 2])
a + b
Out[128]:
在不知道广播的时候已经使用过它！：
In [129]:
a = np.ones((4, 5))
a[0] = 2 # 我们将一个数组的纬度0分配给另一个数组的纬度1
a
Out[129]:
In [130]:
a = np.ones((4, 5))
print a[0]
a[0] = 2 # 我们将一个数组的纬度0分配给另一个数组的纬度
a
Out[130]:
一个有用的技巧：
In [133]:
a = np.arange(0, 40, 10)
a.shape
Out[133]:
In [134]:
a = a[:, np.newaxis] # 添加一个新的轴 -> 2D 数组
a.shape
Out[134]:
In [135]:
a
Out[135]:
In [136]:
a + b
Out[136]:
广播看起来很神奇，但是，当我们要解决的问题是输出数据比输入数据有更多纬度的数组时，使用它是非常自然的。
实例：广播
让我们创建一个66号公路上城市之间距离（用公里计算）的数组：芝加哥、斯普林菲尔德、圣路易斯、塔尔萨、俄克拉何马市、阿马里洛、圣塔菲、阿尔布开克、Flagstaff、洛杉矶。
In [138]:
mileposts = np.array([0, 198, 303, 736, 871, 1175, 1475, 1544, 1913, 2448])
distance_array = np.abs(mileposts - mileposts[:, np.newaxis])
distance_array
Out[138]:
许多基于网格或者基于网络的问题都需要使用广播。例如，如果要计算10X10网格中每个点到原点的数据，可以这样：
In [139]:
x, y = np.arange(5), np.arange(5)[:, np.newaxis]
distance = np.sqrt(x ** 2 + y ** 2)
distance
Out[139]:
或者用颜色：
In [141]:
plt.pcolor(distance)
plt.colorbar()
Out[141]:
评论 : numpy.ogrid 函数允许直接创建上一个例子中两个重要纬度向量X和Y：
In [142]:
x, y = np.ogrid[0:5, 0:5]
x, y
Out[142]:
In [143]:
x.shape, y.shape
Out[143]:
In [144]:
distance = np.sqrt(x ** 2 + y ** 2)
因此， np.ogrid 就非常有用，只要我们是要处理网格计算。另一方面， 在一些无法（或者不想）从广播中收益的情况下，np.mgrid 直接提供了由索引构成的矩阵：
In [145]:
x, y = np.mgrid[0:4, 0:4]
x
Out[145]:
In [146]:
y
Out[146]:
In [147]:
a = np.array([[1, 2, 3], [4, 5, 6]])
a.ravel()
Out[147]:
In [148]:
a.T
Out[148]:
In [149]:
a.T.ravel()
Out[149]:
In [150]:
a.shape
Out[150]:
In [152]:
b = a.ravel()
b = b.reshape((2, 3))
b
Out[152]:
或者：
In [153]:
a.reshape((2, -1)) # 不确定的值（-1）将被推导
Out[153]:
警告： ndarray.reshape 可以返回一个视图（参见 help(np.reshape)）, 也可以可以返回副本
In [155]:
b[0, 0] = 99
a
Out[155]:
当心：重排也可以返回一个副本！：
In [156]:
a = np.zeros((3, 2))
b = a.T.reshape(3*2)
b[0] = 9
a
Out[156]:
In [157]:
z = np.array([1, 2, 3])
z
Out[157]:
In [158]:
z[:, np.newaxis]
Out[158]:
In [159]:
z[np.newaxis, :]
Out[159]:
1.3.2.4.4 纬度重组¶
In [160]:
a = np.arange(4*3*2).reshape(4, 3, 2)
a.shape
Out[160]:
In [161]:
a[0, 2, 1]
Out[161]:
In [163]:
b = a.transpose(1, 2, 0)
b.shape
Out[163]:
In [164]:
b[2, 1, 0]
Out[164]:
也是创建了一个视图：
In [165]:
b[2, 1, 0] = -1
a[0, 2, 1]
Out[165]:
1.3.2.4.5 改变大小¶
可以用 ndarray.resize 改变数组的大小：
In [167]:
a = np.arange(4)
a.resize((8,))
a
Out[167]:
但是，它不能在其他地方引用：
In [168]:
b = a
a.resize((4,))
练习：形状操作
- 看一下
reshape的文档字符串，特别要注意其中关于副本和视图的内容。 - 用
flatten来替换ravel。有什么区别？ (提示: 试一下哪个返回视图哪个返回副本) - 试一下用
transpose来进行纬度变换。
1.3.2.5 数据排序¶
按一个轴排序：
In [169]:
a = np.array([[4, 3, 5], [1, 2, 1]])
b = np.sort(a, axis=1)
b
Out[169]:
注：每行分别排序！
原地排序：
In [170]:
a.sort(axis=1)
a
Out[170]:
象征索引排序：
In [171]:
a = np.array([4, 3, 1, 2])
j = np.argsort(a)
j
Out[171]:
In [172]:
a[j]
Out[172]:
找到最大值和最小值：
In [173]:
a = np.array([4, 3, 1, 2])
j_max = np.argmax(a)
j_min = np.argmin(a)
j_max, j_min
Out[173]:
练习：排序
- 试一下原地和非原地排序
- 试一下用不同的数据类型创建数组并且排序。
- 用
all或者array_equal来检查一下结果。 - 看一下
np.random.shuffle，一种更快创建可排序输入的方式。 - 合并
ravel、sort和reshape。 - 看一下
sort的axis关键字，重写一下这个练习。
1.3.2.6 总结¶
入门你需要了解什么？
- 了解如何创建数组：
array、arange、ones、zeros。 - 了解用
array.shape数组的形状，然后使用切片来获得数组的不同视图：array[::2]等。用reshape改变数组形状或者用ravel扁平化。 - 获得数组元素的一个子集，和/或用面具修改他们的值ay and/or modify their values with masks
In [174]:
a[a < 0] = 0
- 了解数组上各式各样的操作，比如找到平均数或最大值 (
array.max()、array.mean())。不需要记住所有东西，但是应该有条件反射去搜索文档 (线上文档,help(),lookfor())!! - 更高级的用法：掌握用整型数组索引，以及广播。了解更多的Numpy函数以便处理多种数组操作。
快读阅读
如果你想要快速通过科学讲座笔记来学习生态系统，你可以直接跳到下一章：Matplotlib: 作图(暂缺)。
本章剩下的内容对于学习介绍部分不是必须的。但是，记得回来完成本章并且完成更多的练习。
1.3.3 数据的更多内容¶
1.3.3.1 更多的数据类型¶
1.3.3.1.1 投射¶
“更大”的类型在混合类型操作中胜出：
In [175]:
np.array([1, 2, 3]) + 1.5
Out[175]:
赋值不会改变类型！
In [176]:
a = np.array([1, 2, 3])
a.dtype
Out[176]:
In [178]:
a[0] = 1.9 # <-- 浮点被截取为整数
a
Out[178]:
强制投射：
In [179]:
a = np.array([1.7, 1.2, 1.6])
b = a.astype(int) # <-- 截取整数
b
Out[179]:
四舍五入：
In [180]:
a = np.array([1.2, 1.5, 1.6, 2.5, 3.5, 4.5])
b = np.around(a)
b # 仍然是浮点
Out[180]:
In [181]:
c = np.around(a).astype(int)
c
Out[181]:
1.3.3.1.2 不同数据类型的大小¶
整数 (带有符号):
| 类型 | 字节数 |
|---|---|
| int8 | 8 bits |
| int16 | 16 bits |
| int32 | 32 bits (与32位平台的int相同) |
| int64 | 64 bits (与64位平台的int相同) |
In [182]:
np.array([1], dtype=int).dtype
Out[182]:
In [183]:
np.iinfo(np.int32).max, 2**31 - 1
Out[183]:
In [184]:
np.iinfo(np.int64).max, 2**63 - 1
Out[184]:
无符号整数:
| 类型 | 字节数 |
|---|---|
| uint8 | 8 bits |
| uint16 | 16 bits |
| uint32 | 32 bits |
| uint64 | 64 bits |
In [185]:
np.iinfo(np.uint32).max, 2**32 - 1
Out[185]:
In [186]:
np.iinfo(np.uint64).max, 2**64 - 1
Out[186]:
浮点数据：
| 类型 | 字节数 |
|---|---|
| float16 | 16 bits |
| float32 | 32 bits |
| float64 | 64 bits (与浮点相同) |
| float96 | 96 bits, 平台依赖 (与 np.longdouble 相同) |
| float128 | 128 bits, 平台依赖 (与 np.longdouble相同) |
In [187]:
np.finfo(np.float32).eps
Out[187]:
In [188]:
np.finfo(np.float64).eps
Out[188]:
In [189]:
np.float32(1e-8) + np.float32(1) == 1
Out[189]:
In [190]:
np.float64(1e-8) + np.float64(1) == 1
Out[190]:
浮点复数：
| 类型 | 字节数 |
|---|---|
| complex64 | 两个 32-bit 浮点 |
| complex128 | 两个 64-bit 浮点 |
| complex192 | 两个 96-bit 浮点, 平台依赖 |
| complex256 | 两个 128-bit 浮点, 平台依赖 |
更小的数据类型
如果你不知道需要特殊数据类型，那你可能就不需要。
比较使用 float32代替 float64:
- 一半的内存和硬盘大小
- 需要一半的宽带（可能在一些操作中更快）
In [191]:
a = np.zeros((1e6,), dtype=np.float64)
b = np.zeros((1e6,), dtype=np.float32)
%timeit a*a
In [192]:
%timeit b*b
- 但是更大的四舍五入误差 - 有时在一些令人惊喜的地方（即，不要使用它们除非你真的需要）
1.3.3.2 结构化的数据类型¶
| 名称 | 类型 |
|---|---|
| sensor_code | (4个字母的字符) |
| position | (浮点) |
| value | (浮点) |
In [194]:
samples = np.zeros((6,), dtype=[('sensor_code', 'S4'),('position', float), ('value', float)])
samples.ndim
Out[194]:
In [195]:
samples.shape
Out[195]:
In [196]:
samples.dtype.names
Out[196]:
In [198]:
samples[:] = [('ALFA', 1, 0.37), ('BETA', 1, 0.11), ('TAU', 1, 0.13),('ALFA', 1.5, 0.37), ('ALFA', 3, 0.11),
('TAU', 1.2, 0.13)]
samples
Out[198]:
用字段名称索引也可以访问字段：
In [199]:
samples['sensor_code']
Out[199]:
In [200]:
samples['value']
Out[200]:
In [201]:
samples[0]
Out[201]:
In [202]:
samples[0]['sensor_code'] = 'TAU'
samples[0]
Out[202]:
一次多个字段：
In [203]:
samples[['position', 'value']]
Out[203]:
和普通情况一样，象征索引也有效：
In [204]:
samples[samples['sensor_code'] == 'ALFA']
Out[204]:
1.3.3.3 面具数组（maskedarray）: 处理缺失值（的传播）¶
- 对于浮点不能用NaN，但是面具对所有类型都适用：
In [207]:
x = np.ma.array([1, 2, 3, 4], mask=[0, 1, 0, 1])
x
Out[207]:
In [208]:
y = np.ma.array([1, 2, 3, 4], mask=[0, 1, 1, 1])
x + y
Out[208]:
- 通用函数的面具版本：
In [209]:
np.ma.sqrt([1, -1, 2, -2])
Out[209]:
注：有许多其他数组的兄弟姐妹
尽管这脱离了Numpy这章的主题，让我们花点时间回忆一下编写代码的最佳实践，从长远角度这绝对是值得的：
最佳实践
- 明确的变量名（不需要备注去解释变量里是什么）
风格：逗号后及＝周围有空格等。
在Python代码风格指南及文档字符串惯例页面中给出了相当数据量如何书写“漂亮代码”的规则（并且，最重要的是，与其他人使用相同的惯例!）。
- 除非在一些极特殊的情况下，变量名及备注用英文。
1.3.4 高级操作¶
1.3.4.1. 多项式
Numpy也包含不同基的多项式：
例如，$3x^2 + 2x - 1$:
In [211]:
p = np.poly1d([3, 2, -1])
p(0)
Out[211]:
In [212]:
p.roots
Out[212]:
In [213]:
p.order
Out[213]:
In [215]:
x = np.linspace(0, 1, 20)
y = np.cos(x) + 0.3*np.random.rand(20)
p = np.poly1d(np.polyfit(x, y, 3))
t = np.linspace(0, 1, 200)
plt.plot(x, y, 'o', t, p(t), '-')
Out[215]:
更多内容见http://docs.scipy.org/doc/numpy/reference/routines.polynomials.poly1d.html。
1.3.4.1.1 更多多项式（有更多的基）¶
Numpy也有更复杂的多项式接口，支持比如切比雪夫基。
$3x^2 + 2x - 1$:
In [216]:
p = np.polynomial.Polynomial([-1, 2, 3]) # 系数的顺序不同！
p(0)
Out[216]:
In [217]:
p.roots()
Out[217]:
In [218]:
p.degree() # 在普通的多项式中通常不暴露'order'
Out[218]:
在切尔雪夫基中使用多项式的例子，多项式的范围在[-1,1]：
In [221]:
x = np.linspace(-1, 1, 2000)
y = np.cos(x) + 0.3*np.random.rand(2000)
p = np.polynomial.Chebyshev.fit(x, y, 90)
t = np.linspace(-1, 1, 200)
plt.plot(x, y, 'r.')
plt.plot(t, p(t), 'k-', lw=3)
Out[221]:
切尔雪夫多项式在插入方面有很多优势。
1.3.4.2 加载数据文件¶
1.3.4.2.1 文本文件¶
例子: populations.txt:
# year hare lynx carrot
1900 30e3 4e3 48300
1901 47.2e3 6.1e3 48200
1902 70.2e3 9.8e3 41500
1903 77.4e3 35.2e3 38200In [222]:
data = np.loadtxt('data/populations.txt')
data
Out[222]:
In [224]:
np.savetxt('pop2.txt', data)
data2 = np.loadtxt('pop2.txt')
注：如果你有一个复杂的文本文件，应该尝试：
np.genfromtxt- 使用Python的I/O函数和例如正则式来解析（Python特别适合这个工作）
提示：用IPython在文件系统中航行
In [225]:
pwd # 显示当前目录
Out[225]:
In [227]:
cd data
In [228]:
ls
1.3.4.2.2 图像¶
使用Matplotlib：
In [233]:
img = plt.imread('data/elephant.png')
img.shape, img.dtype
Out[233]:
In [234]:
plt.imshow(img)
Out[234]:
In [237]:
plt.savefig('plot.png')
plt.imsave('red_elephant', img[:,:,0], cmap=plt.cm.gray)
这只保存了一个渠道（RGB）：
In [238]:
plt.imshow(plt.imread('red_elephant.png'))
Out[238]:
其他包：
In [239]:
from scipy.misc import imsave
imsave('tiny_elephant.png', img[::6,::6])
plt.imshow(plt.imread('tiny_elephant.png'), interpolation='nearest')
Out[239]:
1.3.4.2.3 Numpy的自有格式¶
Numpy有自有的二进制格式，没有便携性但是I/O高效：
In [240]:
data = np.ones((3, 3))
np.save('pop.npy', data)
data3 = np.load('pop.npy')
1.3.4.2.4 知名的（并且更复杂的）文件格式¶
- HDF5: h5py, PyTables
- NetCDF:
scipy.io.netcdf_file, netcdf4-python, ... - Matlab:
scipy.io.loadmat,scipy.io.savemat - MatrixMarket:
scipy.io.mmread,scipy.io.mmread
... 如果有人使用，那么就可能有一个对应的Python库。
练习：文本数据文件
写一个Python脚本从populations.txt加载数据，删除前五行和后五行。将这个小数据集存入 pop2.txt。
Numpy内部
如果你对Numpy的内部感兴趣, 有一个关于Advanced Numpy的很好的讨论。
1.3.5 一些练习¶
1.3.5.1 数组操作¶
- 从2D数组（不需要显示的输入）:
[[1, 6, 11],
[2, 7, 12],
[3, 8, 13],
[4, 9, 14],
[5, 10, 15]]并且生成一个第二和第四行的新数组。
- 将数组a的每一列以元素的方式除以数组b (提示:
np.newaxis):
In [243]:
a = np.arange(25).reshape(5, 5)
b = np.array([1., 5, 10, 15, 20])
- 难一点的题目：创建10 X 3的随机数数组 （在[0, 1]的范围内）。对于每一行，挑出最接近0.5的数。
- 用
abs和argsort找到每一行中最接近的列j。 - 使用象征索引抽取数字。（提示：a[i,j]-数组
i必须包含j中成分的对应行数）
- 用
1.3.5.2 图片操作：给Lena加边框¶
让我们从著名的Lena图（http://www.cs.cmu.edu/~chuck/lennapg/） 上开始，用Numpy数组做一些操作。Scipy在 scipy.lena函数中提供了这个图的二维数组：
In [244]:
from scipy import misc
lena = misc.lena()
注：在旧版的scipy中，你会在 scipy.lena()找到lena。
这是一些通过我们的操作可以获得图片：使用不同的颜色地图，裁剪图片，改变图片的一部分。
- 让我们用pylab的
imshow函数显示这个图片。
In [245]:
import pylab as plt
lena = misc.lena()
plt.imshow(lena)
Out[245]:
- Lena然后以为色彩显示。要将她展示为灰色需要指定一个颜色地图。
In [246]:
plt.imshow(lena, cmap=plt.cm.gray)
Out[246]: