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

来源：互联网发布：国企股东是谁知乎编辑：程序博客网时间：2024/06/05 08:12

标准安装的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'])>>> earray(['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]])>>> barray([5, 6, 7, 8])>>> carray([[1, 2, 3, 4],       [4, 5, 6, 7],       [7, 8, 9, 10]])>>> c.dtypedtype('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>>> carray([[ 1,  2,  3],       [ 4,  4,  5],       [ 6,  7,  7],       [ 8,  9, 10]])

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

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

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

>>> d = a.reshape((2,2))>>> darray([[1, 2],       [3, 4]])>>> aarray([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+9725185

我们可以写一个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))>>> aarray([[  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    # 下标还可以用来修改元素的值>>> aarray([  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共享同一块数据空间>>> barray([101,   4,   5,   6])>>> b[2] = -10 # 将b的第2个元素修改为-10>>> barray([101,   4, -10,   6])>>> a # a的第5个元素也被修改为10array([  0,   1, 100, 101,   4, -10,   6,   7,   8,   9])

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

>>> x = np.arange(10,1,-1)>>> xarray([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>>> barray([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 # 整数序列下标也可以用来修改元素的值>>> xarray([10, -3,  8, -1,  6, -2,  4,  3,  2])

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

>>> x = np.arange(5,0,-1)>>> xarray([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])]>>> # 布尔数组的长度不够时，不够的部分都当作Falsearray([5, 3, 2])>>> x[np.array([True, False, True, True])] = -1, -2, -3>>> # 布尔数组下标也可以用来修改元素>>> xarray([-1,  4, -2, -3,  1])

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

>>> x = np.random.rand(10) # 产生一个长度为10，元素值为0-1的随机数的数组>>> xarray([ 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.5array([ 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的数组，图中用颜色区分了各个下标以及其对应的选择区域。

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------

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

多维数组同样也可以使用整数序列和布尔数组进行存取。

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 nppersontype = 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.dtypedtype([('name', '|S32'), ('age', '<i4'), ('weight', '<f4')])

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

>>> a[0]('Zhang', 32, 75.5)>>> a[0].dtypedtype([('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"]>>> barray([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)>>> yarray([  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)>>> xarray([  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 timeimport mathimport numpy as npx = [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() - startx = [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)>>> aarray([0, 1, 2, 3])>>> b = np.arange(1,5)>>> barray([1, 2, 3, 4])>>> np.add(a,b)array([1, 3, 5, 7])>>> np.add(a,b,a)array([1, 3, 5, 7])>>> aarray([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 * bx = t + cdel t

也就是说需要产生一个数组t保存乘法的计算结果，然后再产生最后的结果数组x。我们可以通过手工将一个算式分解为x = a*b; x += c，以减少一次内存分配。

通过组合标准的ufunc函数的调用，可以实现各种算式的数组计算。不过有些时候这种算式不易编写，而针对每个元素的计算函数却很容易用Python实现，这时可以用frompyfunc函数将一个计算单个元素的函数转换成ufunc函数。这样就可以方便地用所产生的ufunc函数对数组进行计算了。让我们来看一个例子。

我们想用一个分段函数描述三角波，三角波的样子如图所示：

我们很容易根据上图所示写出如下的计算三角波某点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)>>> aarray([[ 0], [10], [20], [30], [40], [50]])>>> a.shape(6, 1)再创建一维数组b，其shape为(5,)：>>> b = np.arange(0, 5)>>> barray([0, 1, 2, 3, 4])>>> b.shape(5,)

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

>>> c = a + b>>> carray([[ 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>>> barray([[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)>>> barray([[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)>>> aarray([[ 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]>>> xarray([[0],       [1],       [2],       [3],       [4]])>>> yarray([[0, 1, 2, 3, 4]])

ogrid是一个很有趣的对象，它像一个多维数组一样，用切片组元作为下标进行存取，返回的是一组可以用来广播计算的数组。其切片下标有两种形式：
1.开始值:结束值:步长，和np.arange(开始值, 结束值, 步长)类似
2.开始值:结束值:长度j，当第三个参数为虚数时，它表示返回的数组的长度，和np.linspace(开始值, 结束值, 长度)类似：

>>> x, y = np.ogrid[0:1:4j, 0:1:3j]>>> xarray([[ 0.        ],       [ 0.33333333],       [ 0.66666667],       [ 1.        ]])>>> yarray([[ 0. ,  0.5,  1. ]])

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

import numpy as npfrom enthought.mayavi import mlabx, 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曲面

8.ufunc的方法

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

>>> np.add.reduce([1,2,3]) # 1 + 2 + 36>>> np.add.reduce([[1,2,3],[4,5,6]], axis=1) # 1,4 + 2,5 + 3,6array([ 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 方法，<op>.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>>> aarray([[ 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的dtypedtype('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" )>>> carray([[ 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"] # 数组aarray([[1, 2, 3],       [4, 5, 6]])>>> r["arr_1"] # 数组barray([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9])>>> r["sin_array"] # 数组carray([ 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

阅读全文

0 0