hans

hans

【Python】二、Numpy——《用Python做科学计算》


标准安装的 Python 中用列表 (list) 保存一组值,可以用来当作数组使用,不过由于列表的元素可以是任何对象,因此列表中所保存的是对象的指针。这样为了保存一个简单的 [1,2,3],需要有 3 个指针和三个整数对象。对于数值运算来说这种结构显然比较浪费内存和 CPU 计算时间。
此外 Python 还提供了一个 array 模块,array 对象和列表不同,它直接保存数值,和 C 语言的一维数组比较类似。但是由于它不支持多维,也没有各种运算函数,因此也不适合做数值运算。
NumPy 的诞生弥补了这些不足,NumPy 提供了两种基本的对象:ndarray(N-dimensional array object)和
ufunc(universal function
object)。ndarray (下文统一称之为数组) 是存储单一数据类型的多维数组,而 ufunc 则是能够对数组进行处理的函数。#

这里我想讲一下列表(list),元组 or 组元(tuple)和数组(array)的区别:
列表: 元素可以是任何类型的对象,必须初始化,并且初始化后内容可以修改。当所有元素对象一样的时候可以当数组使用。

>>> a = ['a', "hello", 1, 2.2]
>>> a
['a', 'hello', 1, 2.2]
>>> a[1] = "HELLO"
>>> a
['a', 'HELLO', 1, 2.2]

元组: 元素可以是任何类型的对象,必须初始化,并且初始化后内容不可以修改。

>>> b = ('a', "hello", 1, 2.2)
>>> b
('a', 'hello', 1, 2.2)
>>> b[2] = "HELLO"
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'tuple' object does not support item assignment

数组: 元素只可以是一种类型的对象,非必须初始化,并且初始化后内容可以修改。

>>> c = np.array([1,2,3])
>>> d = np.array(['a','b','c'])
>>> e = np.array([1,2,'c','d'])
>>> e
array(['1', '2', 'c', 'd'],
      dtype='|S21')

看到没,数组 e 中将 1,2, 默认为了 string 类型与后面的‘c',’d' 保持一致。
善于观察的话可以看到元组是用的小括号(),列表和数组用的中括号 []。#

1. 函数库导入#

import numpy as np

导入 numpy 函数库,并在后面用缩略 np 表示。

2. 创建#

a,b 是一维数组,c 是多维数组。

>>> 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]])
>>> b
array([5, 6, 7, 8])
>>> c
array([[1, 2, 3, 4],
       [4, 5, 6, 7],
       [7, 8, 9, 10]])
>>> c.dtype
dtype('int32')

数组的大小可以通过其 shape 属性获得:

>>> a.shape
(4,)
>>> c.shape
(3, 4)

数组 a 的 shape 只有一个元素,因此它是一维数组。而数组 c 的 shape 有两个元素,因此它是二维数组,其中第 0 轴(3 个横轴)的长度为 3,第 1 轴(4 个纵轴)的长度为 4。还可以通过修改数组的 shape 属性,在保持数组元素个数不变的情况下,改变数组每个轴的长度。下面的例子将数组 c 的 shape 改为 (4,3),注意从 (3,4) 改为 (4,3) 并不是对数组进行转置,而只是改变每个轴的大小,数组元素在内存中的位置并没有改变:

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

看上面的数组 a 是一维数组,a.shape 只有一个数,这个数字表示数组里面元素的个数。我们还可以用多维数组的表达方法表示 a。注意观察区别:

>>> a = np.array([1, 2, 3, 4])
>>> a
array([1, 2, 3, 4])
>>> a.shape
(4,)
>>> a.shape = 1,4
>>> a
array([[1, 2, 3, 4]])
>>> a.shape
(1, 4)

当某个轴的元素为 - 1 时,将根据数组元素的个数自动计算此轴的长度,因此下面的程序将数组 c 的 shape 改为了 (2,6):

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

使用数组的 reshape 方法,可以创建一个改变了尺寸的新数组,原数组的 shape 保持不变:

>>> d = a.reshape((2,2))
>>> d
array([[1, 2],
       [3, 4]])
>>> a
array([1, 2, 3, 4])

数组 a 和 d 其实共享数据存储内存区域,因此修改其中任意一个数组的元素都会同时修改另外一个数组的内容:

>>> a[1] = 100 # 将数组a的第一个元素改为100
>>> d # 注意数组d中的2也被改变了
array([[  1, 100],
       [  3,   4]])

数组的元素类型可以通过 dtype 属性获得。上面例子中的参数序列的元素都是整数,因此所创建的数组的元素类型也是整数,并且是 32bit 的长整型。可以通过 dtype 参数在创建时指定元素类型:

>>> np.array([[1, 2, 3, 4],[4, 5, 6, 7], [7, 8, 9, 10]], dtype=np.float)
array([[  1.,   2.,   3.,   4.],
       [  4.,   5.,   6.,   7.],
       [  7.,   8.,   9.,  10.]])
>>> np.array([[1, 2, 3, 4],[4, 5, 6, 7], [7, 8, 9, 10]], dtype=np.complex)
array([[  1.+0.j,   2.+0.j,   3.+0.j,   4.+0.j],
       [  4.+0.j,   5.+0.j,   6.+0.j,   7.+0.j],
       [  7.+0.j,   8.+0.j,   9.+0.j,  10.+0.j]])

上面都是通过 array 函数创建数组。下面是三种专门用来创建数组的函数。
arange 函数类似于 python 的 range 函数,通过指定开始值、终值和步长来创建一维数组, 注意数组不包括终值 :

>>> 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 函数通过指定开始值、终值和元素个数来创建一维数组,可以通过 endpoint 关键字指定是否包括终值,缺省(默认)设置是 包括终值 :

>>> np.linspace(0, 1, 12)
array([ 0.        ,  0.09090909,  0.18181818,  0.27272727,  0.36363636,
        0.45454545,  0.54545455,  0.63636364,  0.72727273,  0.81818182,
        0.90909091,  1.        ])

logspace 函数和 linspace 类似,不过它创建等比数列,下面的例子产生 1 (10^0) 到 100 (10^2)、有 20 个元素的等比数列:

>>> np.logspace(0, 2, 20)
array([   1.        ,    1.27427499,    1.62377674,    2.06913808,
          2.6366509 ,    3.35981829,    4.2813324 ,    5.45559478,
          6.95192796,    8.8586679 ,   11.28837892,   14.38449888,
         18.32980711,   23.35721469,   29.76351442,   37.92690191,
         48.32930239,   61.58482111,   78.47599704,  100.        ])

此外,使用 frombuffer, fromstring, fromfile 等函数可以从字节序列创建数组,下面以 fromstring 为例:

>>> s = "abcdefgh"

Python 的字符串实际上是字节序列,每个字符占一个字节,因此如果从字符串 s 创建一个 8bit 的整数数组的话,所得到的数组正好就是字符串中每个字符的 ASCII 编码:

>>> np.fromstring(s, dtype=np.int8)
array([ 97,  98,  99, 100, 101, 102, 103, 104], dtype=int8)

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

>>> np.fromstring(s, dtype=np.int16)
array([25185, 25699, 26213, 26727], dtype=int16)
>>> 98*256+97
25185

我们可以写一个 Python 的函数,它将数组下标转换为数组中对应的值,然后使用此函数创建数组:

>>> def func(i): #定义函数,函数名func,参数i
...   return i%4+1 #返回i除以4的余数加一
...
>>> np.fromfunction(func, (10,))
array([ 1.,  2.,  3.,  4.,  1.,  2.,  3.,  4.,  1.,  2.])

fromfunction 第二个参数 n 可以理解为从 0 开始一共有 n 个取值,赋值给前面的 function。
这个例子是 i 先取 0, 取余加一后返回数值 1
然后 i 取 1, 取余加一后返回数值 2
依次类推
最后 i 取 9,取余加一后返回数值 2#

fromfunction 函数的第一个参数为计算每个数组元素的函数,第二个参数为数组的大小 (shape),因为它支持多维数组,所以第二个参数必须是一个序列,本例中用 (10,) 创建一个 10 元素的一维数组。

下面的例子创建一个二维数组表示九九乘法表,输出的数组 a 中的每个元素 a [i, j] 都等于 func2 (i, j):

>>> def func2(i, j): #定义函数,函数名func,参数i,j
...     return (i+1) * (j+1) #函数返回值
...
>>> a = np.fromfunction(func2, (9,9))
>>> a
array([[  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.],
       [  2.,   4.,   6.,   8.,  10.,  12.,  14.,  16.,  18.],
       [  3.,   6.,   9.,  12.,  15.,  18.,  21.,  24.,  27.],
       [  4.,   8.,  12.,  16.,  20.,  24.,  28.,  32.,  36.],
       [  5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,  45.],
       [  6.,  12.,  18.,  24.,  30.,  36.,  42.,  48.,  54.],
       [  7.,  14.,  21.,  28.,  35.,  42.,  49.,  56.,  63.],
       [  8.,  16.,  24.,  32.,  40.,  48.,  56.,  64.,  72.],
       [  9.,  18.,  27.,  36.,  45.,  54.,  63.,  72.,  81.]]) #可以理解为每一点的下标相乘

3. 存取元素#

数组元素的存取方法和 Python 的标准方法相同:

>>> a = np.arange(10)
>>> a[5]    # 用整数作为下标可以获取数组中的某个元素
5
>>> a[3:5]  # 用范围作为下标获取数组的一个切片,包括a[3]不包括a[5]
array([3, 4])
>>> a[:5]   # 省略开始下标,表示从a[0]开始
array([0, 1, 2, 3, 4])
>>> a[:-1]  # 下标可以使用负数,表示从数组最后一个元素往前数
array([0, 1, 2, 3, 4, 5, 6, 7, 8])
>>> a[2:4] = 100,101    # 下标还可以用来修改元素的值
>>> a
array([  0,   1, 100, 101,   4,   5,   6,   7,   8,   9])
>>> a[1:-1:2]   # 范围中的第三个参数表示步长,2表示隔一个元素取一个元素
array([  1, 101,   5,   7])
>>> a[::-1] # 省略范围的开始下标和结束下标,步长为-1,整个数组头尾颠倒
array([  9,   8,   7,   6,   5,   4, 101, 100,   1,   0])
>>> a[5:1:-2] # 步长为负数时,开始下标必须大于结束下标
array([  5, 101])

和 Python 的列表序列不同,通过下标范围获取的新的数组是原始数组的一个视图。它与原始数组共享同一块数据空间:

>>> b = a[3:7] # 通过下标范围产生一个新的数组b,b和a共享同一块数据空间
>>> b
array([101,   4,   5,   6])
>>> b[2] = -10 # 将b的第2个元素修改为-10
>>> b
array([101,   4, -10,   6])
>>> a # a的第5个元素也被修改为10
array([  0,   1, 100, 101,   4, -10,   6,   7,   8,   9])

除了使用下标范围存取元素之外,NumPy 还提供了两种存取元素的高级方法。
使用整数序列
当使用整数序列对数组元素进行存取时,将使用整数序列中的每个元素作为下标,整数序列可以是列表或者数组。使用整数序列作为下标获得的数组不和原始数组共享数据空间。

>>> x = np.arange(10,1,-1)
>>> x
array([10,  9,  8,  7,  6,  5,  4,  3,  2])
>>> x[[3, 3, 1, 8]] # 获取x中的下标为3, 3, 1, 8的4个元素,组成一个新的数组
array([7, 7, 9, 2])
>>> b = x[np.array([3,3,-3,8])]  #下标可以是负数
>>> b[2] = 100
>>> b
array([7, 7, 100, 2])
>>> x   # 由于b和x不共享数据空间,因此x中的值并没有改变
array([10,  9,  8,  7,  6,  5,  4,  3,  2])
>>> x[[3,5,1]] = -1, -2, -3 # 整数序列下标也可以用来修改元素的值
>>> x
array([10, -3,  8, -1,  6, -2,  4,  3,  2])

使用布尔数组
当使用布尔数组 b 作为下标存取数组 x 中的元素时,将收集数组 x 中所有在数组 b 中对应下标为 True 的元素。使用布尔数组作为下标获得的数组不和原始数组共享数据空间,注意这种方式只对应于布尔数组,不能使用布尔列表。(列表,元组,数组区别还记得?)

>>> x = np.arange(5,0,-1)
>>> x
array([5, 4, 3, 2, 1])
>>> x[np.array([True, False, True, False, False])]
>>> # 布尔数组中下标为0,2的元素为True,因此获取x中下标为0,2的元素
array([5, 3])
>>> x[[True, False, True, False, False]]
>>> # 如果是布尔列表,则把True当作1, False当作0,按照整数序列方式获取x中的元素
array([4, 5, 4, 5, 5])
>>> x[np.array([True, False, True, True])]
>>> # 布尔数组的长度不够时,不够的部分都当作False
array([5, 3, 2])
>>> x[np.array([True, False, True, True])] = -1, -2, -3
>>> # 布尔数组下标也可以用来修改元素
>>> x
array([-1,  4, -2, -3,  1])

布尔数组一般不是手工产生,而是使用布尔运算的 ufunc 函数产生,关于 ufunc 函数请参照 ufunc 运算 一节。

>>> x = np.random.rand(10) # 产生一个长度为10,元素值为0-1的随机数的数组
>>> x
array([ 0.72223939,  0.921226  ,  0.7770805 ,  0.2055047 ,  0.17567449,
        0.95799412,  0.12015178,  0.7627083 ,  0.43260184,  0.91379859])
>>> x>0.5
>>> # 数组x中的每个元素和0.5进行大小比较,得到一个布尔数组,True表示x中对应的值大于0.5
array([ True,  True,  True, False, False,  True, False,  True, False,  True], dtype=bool)
>>> x[x>0.5]
>>> # 使用x>0.5返回的布尔数组收集x中的元素,因此得到的结果是x中所有大于0.5的元素的数组
array([ 0.72223939,  0.921226  ,  0.7770805 ,  0.95799412,  0.7627083 ,
        0.91379859])

4. 多维数组#

多维数组的存取和一维数组类似,因为多维数组有多个轴,因此它的下标需要用多个值来表示,NumPy 采用组元 (tuple) 作为数组的下标。如图所示,a 为一个 6x6 的数组,图中用颜色区分了各个下标以及其对应的选择区域。
1668631158233.jpg


a [2::2,::2] 这一句分解开这么看,首先 “2:” 表示从索引为 2 的横轴(0 轴)开始到末尾,然后 “:2” 表示隔两行取一次。其实就是取 2,4 行。后面第一个 “:” 表示所有列(1 轴),“:2” 表示隔两列取一次。就是取 0,2,4 三列。

比如 24 的下标是 (2,4),2 就是 0 轴上的索引,4 就是 1 轴上的索引。


数组 a 是一个加法表,纵轴的值为 0, 10, 20, 30, 40, 50;横轴的值为 0, 1, 2, 3, 4, 5。纵轴的每个元素都和横轴的每个元素求和。

>>> np.arange(0, 60, 10).reshape(-1, 1) + np.arange(0, 6)
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]])

多维数组同样也可以使用整数序列和布尔数组进行存取。
1668631278744.jpg
a[(0,1,2,3,4),(1,2,3,4,5)] :
用于存取数组的下标和仍然是一个有两个元素的组元,组元中的每个元素都是整数序列,分别对应数组的第 0 轴和第 1 轴。从两个序列的对应位置取出两个整数组成下标:
a[0,1], a[1,2], ..., a[4,5]。
a [3:, [0, 2, 5]] : 下标中的第 0 轴是一个范围,它选取第 3 行之后的所有行;第 1 轴是整数序列,它选取第 0, 2, 5 三列。
a [mask, 2] : 下标的第 0 轴是一个布尔数组,它选取第 0,2,5 行;第 1 轴是一个整数,选取第 2 列。

5. 数组结构#

假设我们需要定义一个结构数组,它的每个元素都有 name, age 和 weight 字段。在 NumPy 中可以如下定义:

import numpy as np
persontype = np.dtype({
    'names':['name', 'age', 'weight'],
    'formats':['S32','i', 'f']})
a = np.array([("Zhang",32,75.5),("Wang",24,65.2)],
    dtype=persontype)

我们先创建一个 dtype 对象 persontype,通过其字典参数描述结构类型的各个字段。字典有两个关键字:names,formats。每个关键字对应的值都是一个列表。names 定义结构中的每个字段名,而 formats 则定义每个字段的类型:
S32 : 32 个字节的字符串类型,由于结构中的每个元素的大小必须固定,因此需要指定字符串的长度
i : 32bit 的整数类型,相当于 np.int32
f : 32bit 的单精度浮点数类型,相当于 np.float32
然后我们调用 array 函数创建数组,通过关键字参数 dtype=persontype,
指定所创建的数组的元素类型为结构 persontype。运行上面程序之后,我们可以在 IPython 中执行如下的语句查看数组 a 的元素类型

>>> a.dtype
dtype([('name', '|S32'), ('age', '<i4'), ('weight', '<f4')])

结构数组的存取方式和一般数组相同,通过下标能够取得其中的元素,注意元素的值看上去像是组元,实际上它是一个结构:

>>> a[0]
('Zhang', 32, 75.5)
>>> a[0].dtype
dtype([('name', '|S32'), ('age', '<i4'), ('weight', '<f4')])

a [0] 是一个结构元素,它和数组 a 共享内存数据,因此可以通过修改它的字段,改变原始数组中的对应字段:

>>> c = a[1]
>>> c["name"] = "Li"
>>> a[1]["name"]
"Li"

结构像字典一样可以通过字符串下标获取其对应的字段值:

>>> a[0]["name"]
'Zhang'

我们不但可以获得结构元素的某个字段,还可以直接获得结构数组的字段,它返回的是原始数组的视图,因此可以通过修改 b [0] 改变 a [0]["age"]:

>>> b=a[:]["age"] # 或者a["age"]
>>> b
array([32, 24])
>>> b[0] = 40
>>> a[0]["age"]
40

通过调用 a.tostring 或者 a.tofile 方法,可以直接输出数组 a 的二进制形式:

>>> a.tofile("test.bin")

6.ufunc 运算#

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

>>> 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.44921271e-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.44921271e-16])
>>> 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 xrange(1000000)]
start = time.clock()
for i, t in enumerate(x):
    x[i] = math.sin(t)
print "math.sin:", time.clock() - start

x = [i * 0.001 for i in xrange(1000000)]
x = np.array(x)
start = time.clock()
np.sin(x,x)
print "numpy.sin:", time.clock() - start

# 输出
# math.sin: 1.15426932753
# numpy.sin: 0.0882399858083

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

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

# 输出
# numpy.sin loop: 5.72166965355

请注意 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))
<type 'float'>
>>> type(np.sin(0.5))
<type 'numpy.float64'>

通过上面的例子我们了解了如何最有效率地使用 math 库和 numpy 库中的数学函数。
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 函数对数组进行计算了。让我们来看一个例子。

我们想用一个分段函数描述三角波,三角波的样子如图所示:
1668631406858.jpg
我们很容易根据上图所示写出如下的计算三角波某点 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 函数只能计算单个数值,不能对数组直接进行处理。
另外还有一种方法:

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) 将其转换为双精度浮点数组。

7. 广播#

当我们使用 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 是一个很有趣的对象,它像一个多维数组一样,用切片组元作为下标进行存取,返回的是一组可以用来广播计算的数组。其切片下标有两种形式:
1. 开始值:结束值:步长,和 np.arange (开始值,结束值,步长) 类似
2. 开始值:结束值:长度 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. ]])

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

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 库快速绘制如图所示的 3D 曲面
1668631471880.jpg

8.ufunc 的方法#

reduce
方法和 Python 的 reduce 函数类似,它沿着 axis 轴对 array 进行操作,相当于将运算符插入到沿 axis 轴的所有子数组或者元素当中。
例如:

>>> np.add.reduce([1,2,3]) # 1 + 2 + 3
6
>>> 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]])

outer 方法,.outer (a,b) 方法:

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

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

9. 矩阵运算#

矩阵的乘积可以使用 dot 函数进行计算。对于二维数组,它计算的是矩阵乘积,对于一维数组,它计算的是其点积。当需要将一维数组当作列矢量或者行矢量进行矩阵运算时,推荐先使用 reshape 函数将一维数组转换为二维数组:

>>> a = array([1, 2, 3])
>>> a.reshape((-1,1))
array([[1],
       [2],
       [3]])
>>> a.reshape((1,-1))
array([[1, 2, 3]])

dot : 对于两个一维的数组,计算的是这两个数组对应下标元素的乘积和 (数学上称之为内积);其实就是点乘!

10. 文件存储#

NumPy 提供了多种文件操作函数方便我们存取数组内容。文件存取的格式分为两类:二进制和文本。而二进制格式的文件又分为 NumPy 专用的格式化二进制类型和无格式类型。

使用数组的方法函数 tofile 可以方便地将数组中数据以二进制的格式写进文件。tofile 输出的数据没有格式,因此用 numpy.fromfile 读回来的时候需要自己格式化数据:

>>> a = np.arange(0,12)
>>> a.shape = 3,4
>>> a
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
>>> a.tofile("a.bin")
>>> b = np.fromfile("a.bin", dtype=np.float) # 按照float类型读入数据
>>> b # 读入的数据是错误的
array([  2.12199579e-314,   6.36598737e-314,   1.06099790e-313,
         1.48539705e-313,   1.90979621e-313,   2.33419537e-313])
>>> a.dtype # 查看a的dtype
dtype('int32')
>>> b = np.fromfile("a.bin", dtype=np.int32) # 按照int32类型读入数据
>>> b # 数据是一维的
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])
>>> b.shape = 3, 4 # 按照a的shape修改b的shape
>>> b # 这次终于正确了
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

从上面的例子可以看出,需要在读入的时候设置正确的 dtype 和 shape 才能保证数据一致。并且 tofile 函数不管数组的排列顺序是 C 语言格式的还是 Fortran 语言格式的,统一使用 C 语言格式输出。

此外如果 fromfile 和 tofile 函数调用时指定了 sep 关键字参数的话,数组将以文本格式输入输出。

numpy.load 和 numpy.save 函数以 NumPy 专用的二进制类型保存数据,这两个函数会自动处理元素类型和 shape 等信息,使用它们读写数组就方便多了,但是 numpy.save 输出的文件很难和其它语言编写的程序读入:

>>> np.save("a.npy", a)
>>> c = np.load( "a.npy" )
>>> c
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

如果你想将多个数组保存到一个文件中的话,可以使用 numpy.savez 函数。savez 函数的第一个参数是文件名,其后的参数都是需要保存的数组,也可以使用关键字参数为数组起一个名字,非关键字参数传递的数组会自动起名为 arr_0,
arr_1,
...。savez 函数输出的是一个压缩文件 (扩展名为 npz),其中每个文件都是一个 save 函数保存的 npy 文件,文件名对应于数组名。load 函数自动识别 npz 文件,并且返回一个类似于字典的对象,可以通过数组名作为关键字获取数组的内容:

>>> a = np.array([[1,2,3],[4,5,6]])
>>> b = np.arange(0, 1.0, 0.1)
>>> c = np.sin(b)
>>> np.savez("result.npz", a, b, sin_array = c)
>>> r = np.load("result.npz")
>>> r["arr_0"] # 数组a
array([[1, 2, 3],
       [4, 5, 6]])
>>> r["arr_1"] # 数组b
array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9])
>>> r["sin_array"] # 数组c
array([ 0.        ,  0.09983342,  0.19866933,  0.29552021,  0.38941834,
        0.47942554,  0.56464247,  0.64421769,  0.71735609,  0.78332691])

如果你用解压软件打开 result.npz 文件的话,会发现其中有三个文件:arr_0.npy, arr_1.npy,
sin_array.npy,其中分别保存着数组 a, b, c 的内容。

使用 numpy.savetxt 和 numpy.loadtxt 可以读写 1 维和 2 维的数组:

>>> a = np.arange(0,12,0.5).reshape(4,-1)
>>> np.savetxt("a.txt", a) # 缺省按照'%.18e'格式保存数据,以空格分隔
>>> np.loadtxt("a.txt")
array([[  0. ,   0.5,   1. ,   1.5,   2. ,   2.5],
       [  3. ,   3.5,   4. ,   4.5,   5. ,   5.5],
       [  6. ,   6.5,   7. ,   7.5,   8. ,   8.5],
       [  9. ,   9.5,  10. ,  10.5,  11. ,  11.5]])
>>> np.savetxt("a.txt", a, fmt="%d", delimiter=",") #改为保存为整数,以逗号分隔
>>> np.loadtxt("a.txt",delimiter=",") # 读入的时候也需要指定逗号分隔
array([[  0.,   0.,   1.,   1.,   2.,   2.],
       [  3.,   3.,   4.,   4.,   5.,   5.],
       [  6.,   6.,   7.,   7.,   8.,   8.],
       [  9.,   9.,  10.,  10.,  11.,  11.]])

本文参考自: http://old.sebug.net/paper/books/scipydoc/numpy_intro.html

Loading...
Ownership of this post data is guaranteed by blockchain and smart contracts to the creator alone.