NumPy学习笔记

http://blog.chinaunix.net/uid-7596337-id-125813.html 

理解多维矩阵的"求和"、"平均"操作确实太恶心了，numpy提供的函数里还有一堆参数，搞得晕头转向的，这里做个笔记，提醒一下自己， 下面是例程

import numpy as np
X = np.array([[1, 2], [4, 5], [7, 8]])
print np.mean(X, axis=0, keepdims=True)
print np.mean(X, axis=1, keepdims=True)

结果是分别是

[[ 1.5]
[[ 4.  5.]]      [ 4.5]
[ 7.5]]

我个人比较raw的认识就是，axis=0，那么输出矩阵是1行，求每一列的平均（按照每一行去求平均）；axis=1，输出矩阵是1列，求每一行的平均（按照每一列去求平均）。还可以这么理解，axis是几，那就表明哪一维度被压缩成1。

再举个更复杂点的例子，比如我们输入为batch = [128, 28, 28]，可以理解为batch=128，图片大小为28×28像素，我们相求这128个图片的均值，应该这么写

m = np.mean(batch, axis=0)

输出结果m的shape为(28,28)，就是这128个图片在每一个像素点平均值。

NumPy学习笔记（一） 2011-01-07
11:16:26

分类： LINUX

读的是NumPy User Guide Release 1.3，官方另有一本Reference，600多页，还是先读这个省事儿的。

凡例：a. [float]表示对象类型是float，用时不用加'['和']'。b. 代码中的跳格表示换行。

Chap 2 NumPy Basics

1. import numpy as np，沿用之，很喜欢np这个名字，呵呵

2. np支持的数据类型：int8, int16, int32, int64, unit8（无符号的整数，unit8是0～255）, unit16~64, float32, 64, complex64, 128, boolean。字母表示数字类型（dtype），后面的数字表示位数，即数值在内存中所占的bit的数目（不明白～）。还有int, float和complex，其位数取决于平台类型，即16或者32位，complex也是，不过因为是两个float表示的，所以位数是64或者128。还有更高级的类型。

3. 转换数据类型：y=np.int_([1,2,3])，返回一个np的array。

4. 生成array时指定数据类型：z=np.arange(3,dtype='f'/float)，推荐后一种写法。此处虽然写的float，但np会自动把dtype转换成np.flaot。z的转换有两种做法，或者通过类方法：z.astype(int8)，或者通过函数：np.int8(z)。

5. 查看数据类型。也有两种做法，或通过类方法：z.dtype()，或通过函数：np.dtype(z)，或者np.issubtype(d,int)，返回bool值。

6. 创建array的五种方法：a. 从list转换（tuple特殊对待，不转换，见26）；b. 内建函数创建（arange，ones，zeros（这个构建structured array时常用，见30）等等）；c. 从文件读取；d. 直接从raw bytes构建（啥意思？）；e. 借助某些库里的函数生成。

7. 通过转换生成array：x=np.array([1,2,3],[4,5,6])

8. 通过内建方法生成array：np.zeros((2,3))，返回array([0,0,0],[0,0,0])，即维度为2*3的矩阵。

9. arange，类似range，参数可以是小数或者负数，很实用！

10. np.linspaces，规定起点、终点（包含）、返回array的长度，返回一个两端点间数值平均分布的array。

11. np.indices没看懂，似乎某些情况下可以很有用。

12. 直接从文件读取需借用其他库里的方法，不熟，过。

13. x.shape = (2,5)，指定维度，相当于R里的c(2,5) <- dim(x)。也可以这样用：y = x.reshape(5,7)。

14. array indexing和slicing跟Python里的一样也是从0开始的。不同的是如果用index array去做reference的话返回的是一个复本，而不是目标array里的原数。这个原则适用于下面的数值和mask index array。这个特性使得赋值（assignment）运算有点奇特，见24。

15. 说到indexing，记一笔前几天新学到的逆序排列的一招：x[::-1]，其中x不限于list，只要是能slice的都行。其原理是x[[int1]:[int2]:[int3]]，其中int1是起点，int2终点（不包含），int3是步长，负数则表示倒着来。如果int3<0 and int1<int2，那么返回一个len为0的对象，type和x一样。< span="" style="word-wrap: break-word;">

16. array indexing的输入对象也要求是一个array，比如：x[np.array([3, 3, 1, 8])]。也可以是：x[np.array([[1,1],[2,3]])]，返回一个同样维度的array。

17. 多维array的indexing比较麻烦一些，个人认为基本方法就不大直观：y[np.array([0,2,4]), np.array([0,1,2])]返回的是y[0,0]，y[2,1]和y[4,2]。

18. 当然可能也有好处，比如如果取一行/列就为其中每个数指定行/列数，因为有自动broadcasting（不知道啥意思，猜测不止扩展或补齐的意思）功能：y[np.array([0,2,4]), 1]。

19. 不指定某一个维度的index则取该维度所有值，如：y[np.array([0,2,4])]，返回array的shape是index array的shape串联上所有没指定的维度的部分的shape。如这个例子中，y = np.arange(35).reshape(5,7)，shape是5*7，没有指定第二个维度则返回array的shape是3（index array的shape）*7（所有没指定的维度的部分，即第二个维度，的shape）。

20. array数据筛选。我管这个叫filtering，User Guide上叫mask index arrays。一个指功用，一个指手段。用法：y[y>=20]，跟R里的用法一样。也可以绕一点，b=y>20，这样生成一个dtype=bool的array，然后又套在y里：y[b]。还可以对b进行slicing，作为index array指定位置，如：y[b[:,5]]，其中b[:,5]返回第一维度不限，第二维度是5的array，由于y是5*7的，返回array长度是5，shape是5*1。把这个array作为index
array对y筛值，由于是个一维array，对于y来说第二维没有指定，所以取全部的。b是array([False, False, False, True, True], dtype=bool)，所以仅取第一维度的最后两个，对于二维array来说就是最后两行，每行7个数，返回的是一个2*7的array。好绕啊，要是多于2维估计我就晕了~

21. index array和slicing的混搭：y[np.array([0,2,4]),1:3]，工作原理还是跟index array的差不多，只不过1:3这个部分对第一个array里的每个数所代表的行都起作用（broadcast了）。也可以跟mask index array一起混搭。

22. 插入维度：y[:,np.newaxis,:].shape成了(5,1,7)。新增的维度里面没有元素，但运算时行为遵照新的array的适用原则，合并array时很有用：（x原为1*5的array）x[:,np.newaxis] + x[np.newaxis,:] ，返回的是一个5*5的array。不然得赋值z，再指定z的shape等等，较为麻烦。还有一种比较无语的写法：（z的shape是(3,3,3,3)）z[1,...,2]，相当于z[1,:,:,2]。（懒惰是向上的车轮！）

23. array slice赋值，跟Python差不多，需要注意的是每个array都有自己的dtype，赋值类型不符的话会“降级”或者干脆“出乱子”（raise exception啦）。

24. 赋值（assignment）的时候和reference时不一样，总是array里的原数发生变化（不然还赋个啥劲儿的呀），不过规则有点奇特。比如：x[np.array([1, 1, 3, 1]) += 1，array并非直接赋值，而是先被提取出来（extracted），然后逐个赋值，然后assign回去，在1这个位置的那个值被赋值3次，每次都是+=1之后的那个值，所以结果是1和3位置的值+=1。纠结呀~

25. 也可以将index array赋值给变量，然后用变量做indexing，就像20里的那样。可以混搭，可以懒人：indices = (1, Ellipsis, 1)，然后z[indices]

26. list和tuple做indexing结果是不一样的，比如z[(1,1,1,1)]相当于z[1,1,1,1]，仅返回一个数，z[[1,1,1,1]]则返回一个4*3*3*3的array（为什么？z才3*3*3*3，怎么会比z的还多？？）。

27. broadcast的作用是使得不同shape的array在算术运算的时候shape相符，机制上broadcast使得循环在C中发生，这样就比在Python里发生要快。缺点是额外占用内存，所以不要滥用啦。

28. 想要能够利用broadcast必须满足两个条件：从后向前匹配，同一维度上的数值相等，或者至少其中为1。维度数目不一定要求一样。如果是一维的则会自动往另外一个array数值相等或者为1的那个维度上靠。

29. 预设结构模板的array（structured array或称record array）：元素是一组对象组成的tuple，类型已经指定。

30. structured array的dtype参数类型有四种：a. 字符串，包括b1, i1, i2, i4, i8, u1, u2, u4, u8, f4, f8, c8, c16, a，分别代表bytes，整数，无符号整数，浮点数和byte长度一定的字符串，还有一些数：int8,...,uint8,...,float32, float64, complex64, complex128。定义的时候很灵活，比如：x = np.zeros(3, dtype=’3int8, float32, (2,3)float64’)，返回的array长度为3，元素是array([([0,
0, 0], 0.0, [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])。

先告一段落吧，ms User Guide看完第二章就没我啥事儿了，后面都是扩展之类的，离那境界还挺远呢。

31. 第二个参数类型是tuple，仅在将structure映射到一个已有的数据类型上时使用，即在tuple里已有的数据类型以及数据类型定义，比如：x = zeros(3, dtype=(’i4’,[(’r’,’u1’), (’g’,’u1’), (’b’,’u1’), (’a’,’u1’)]))，这样就覆盖了默认名称（f0，f1之类），不过返回没搞明白，为啥是一个1*3的array？行为有点像字典，比如可以这样引用：x[’r’]（默认的话是x[’f1’]），返回的是原数据，而非拷贝。

32. 第三个类型是列表，元素是tuple，每个tuple里，可以有名称、类型、shape（可选）的定义。用法：x = np.zeros(3, dtype=[(’x’,’f4’),(’y’,np.float32),(’value’,’f4’,(2,2))])。

33. 更加复杂的是字典，有两种定义方式。一是给键名names和formats指定长度相等的列表为键值，其中names必须为字符串。还有两个可选的键名offsets和titles。例子：x = np.zeros(3, dtype={’names’:[’col1’, ’col2’], ’formats’:[’i4’,’f4’]})。另一种方式是键名是name的值， 键值其他键名的值是组成的tuple。这里offset是必需的。例子：x = np.zeros(3, dtype={’col1’:(’i1’,0,’title
1’), ’col2’:(’f4’,1,’title 2’)})。

34. 查看name：x.dtype.names（dtype是个对象）。也可以通过赋值进行调整：x.dtype.names = (’x’, ’y’)。还可以像字典一样查看field：x.dtype.fields[’x’][2]。

35. 接下来的内容跟“类”这个概念有关，我到现在也没搞明白，当然也没怎么真正努力去试图理解它。这次就硬着头皮逼一逼自己，看能不能得到一鳞半爪的。首先，子类化ndarray返回的是一个拷贝，正因为此子类化的方法与标准的做法不太一样。

36. 比如初始化时最开始用__new__，然后通常是调用__init__，二者都可以初始化类实例，不过ndarray没有__init__方法，所有初始化都是由__new__完成的。使用__new__的好处是可以返回其他类的实例，而__init__不行。

37. 子类对象的属性需要继承原对象的，但是__new__不能完成这个功能，所以需要__array_finalize__方法把原对象的属性转移过来。

38. 如果用ufunc（通用函数，universal function）处理自定的和ndarray类实例的话，会调用到__array_wrap__方法，把结果转换成子类的实例。缺省状态下__array_wrap__调用__array_finalize__。这个方法需要至少一个参数，是ufunc处理的对象，还有一个可选的参数parameter。有的ufunc返回一个三元tuple，包括ufunc名称、参数和domain，作为parameter。

39. np通过base属性来追踪对象在内存中的位置，这样可以使得在删除view的赋值是不删除原数据本事。原来的array的base一般是None，如果是view，那么其base是所view的原array的名称。

User Guide目前需要看的到此为止，感觉净是array的介绍，没讲什么ufunc，可能是相关的部分还没写吧。网上找了一些资料，改日再看。另外，现在开始看那些Python教程里关于类和OOP的部分应该能加深理解了。

NumPy教程，中文，提供pdf

ndarray类内容一览表

NumPy数学函数一览表

《用Python做科学计算》做得很漂亮，字体瞧着顺眼儿，代码部分不仅照例是等宽字符，还换了灰背景，作者是个有心人。今天答辩委员会的专家还提到了写论文要细致严谨，除了文字之外的其他细节包括图表的刻度、图例和图注还有论文字体大小对齐等等都要力争做到至少规范。因此在投入程度方面儿我要向作者致敬并学习。当然，还有他/她的开源精神~

从目录看只准备看NumPy，SciPy和matplotlib的部分。此外比较感兴趣的是3D演示方面的，因为之前没怎么接触过，不过那部分等先学biopython和wxPython、然后继续修炼一阵Python再说，现在还远照顾不到。吾生也有涯！

凡例：a. [float]表示对象类型是float，用时不用加'['和']'。

Chap 2 NumPy-快速处理数据

1. 两个科学计算Python合集：Python(x,y)和Enthought Python Distribution (EPD)。交互解释器：iPython（装了一个，偶尔用用，还没体会到比Python自带的shell的强大所在及程度）、spyder（Python(x,y)自带的IDE，模仿MATLAB的workspace功能）。比较有用的函数库：Numpy，SciPy，SymPy（符号运算系虾米？），Traits（界面设计库，说得我都想弃wxPython开始学这个了。在计算机领域学什么是个很重要的问题。），TVTK（用Traits库封装的标准VTK库），Visual（制作3D动画演示），OpenCV（这个前几天还碰到过，感觉跟图像生成渲染之类的有关，对于Python，也提供了API）。

2. 调整shape时数组（就是官方文档里的array）元素在内存中的位置并没有改变，也就是提供了一种view而非copy。如果某一维数值是-1，则按数组的长度和其他维度的数值自动计算。

3. reshape方法也是，返回一个shape修改过的数组，但是还是view而非copy。

4. 和linspace相似，logspace创建在log尺度均匀分布的数值组成的数列，即等比数列。

5. 此外，还可以使用frombuffer, fromstring, fromfile等函数可以从字节序列创建数组，也就是上次没有深入的领域。以fromstring为例，在内存中以字节方式存储。一个字节（byte）占8个bit（终于搞明白byte和bit的关系啦，呼瑞！），如果要转化成int16的话有个两个字节串联顺序的问题。Python里是以little endian（又见little endian！当年写tiqho的时候没整明白的问题终于搞明白了，happy）方式在内存中保存数据的，即低位字节在前，故int16值等于之前的字节*256（应该是因为是一个byte占8个bit，故用2^8，不过还是没彻底明白）加上后一个字节。

6. 还可以使用fromfunction利用自己写的函数进行转换。比如：np.fromfunction(func, (10,))，其中func是自写函数，(10,)是返回数组的shape。func的参数默认从0开始，直到填满指定的shape。

7. 数组的reference还可以用boolean array（一般不是手动一个个敲出来的，见9），此时必须是array，不能再用列表蒙混了，broadcast时都按False算。返回其中为True的元素的相应位置上的目标array中的元素。

8. np.random.rand([int]) ：产生一个长度为[int]，元素值为0-1的随机数的数组。

9. 对于数组，如果进行大小比较，和Python里类似，返回的也是bool，不过是个数组。这个数组就可以用来做7里的reference。见7。

10. seriously？——“组元不需要圆括号虽然我们经常在Python中用圆括号将组元括起来，但是其实组元的语法定义只需要用逗号隔开即可。

11. 结构数组（structured array）有点像类（或者我理解的类），定义一个dtype对象，有一定模式，然后指定其他对象的dtype是它。

12. 查看dtype：obj.dtype，可能返回如下：dtype([('name', '|S32'), ('age', '<i4'), ('weight',="" '

13. 作者认为“在C语言中我们可以通过struct关键字定义结构类型，结构中的字段占据连续的内存空间，每个结构体占用的内存大小都相同”。受教了。

14. 通过调用数组的tostring或者tofile方法，可以直接输出数组的二进制形式。比如：a.tofile("test.bin")。

15. 和C语言协作使用时需要注意的是C语言的结构体为了内存寻址方便，会自动的添加一些填充用的字节，即内存对齐，以保持最终的结构体大小变。因此如果NumPy中的所配置的内存大小不符合C语言的对齐规范的话，将会出现数据错位。在创建dtype对象时，可以传递参数align=True以确保不发生数据错位。

16. 用字典参数也可以定义dtype对象时，因为字典的关键字是没有顺序的，所以需要在类型描述中给出字段的顺序。用法是在字段（field）后给出偏移量（即offset）的值，单位是byte，如：dtype([('surname', '|S25'), ('age', '|u1')])。

17. ndarray在内存中的数据结构太深，不过为了将来万一和其他语言的串接，硬着头皮看一看吧。ndarray数据结构引用两个对象：数据存储区和dtype对象存储区，具体包括dtype，dim count，dimensions，strides和data。dim count指维度的数目；dimesion指各维度的数字；strides指每个轴的下标增加1时数据存储区中的指针所增加的字节数，比如有个3*3的数组，元素类型是float32，那么每个数占4（32/8）字节。在C语言格式中，第二个维度的数字比第一个变得快，所以第二个维度数字增加1指针增加4字节，第一个的话就是12（4*3）个字节，所以strides（本身就有步伐”的意思）分别为12和4。不过这要求数据在内存中连续存储。另一种数字跳的方式是Fortran语言格式的，需要重设参数值，改变默认设定：order=''F''。

18. ufunc原来是说能对数组里每个元素都操作的函数，ft~

19. ufunc创建新数组并返回，如果要覆盖的话可以在第二个参数里指定被覆盖的对象，比如：t = np.sin(x,x)，x表现为经过了更新。

20. 算数计算运算符和函数，其中[, y]是可选参数，用于指定覆盖对象，除号运算符的处理根据是否激活__future__.division有所不同：x1+x2，相当于add(x1, x2 [, y])；x1-x2，相当于subtract(x1, x2 [, y])；x1*x2，multiply(x1, x2 [, y])；x1/x2: divide(x1, x2 [, y])，整数除法；x1/x2: true divide (x1, x2 [, y]), 返回精确的商；x1//x2，floor divide
(x1, x2 [, y])，地板除；-x，negative(x [,y])；x1**x2，power(x1, x2 [, y])；x1%x2，remainder(x1, x2 [, y]), mod(x1, x2, [, y])。

21. 复杂的算式处理大数组会产生大量的中间结果而导致速度减缓，一个小技巧是手工拆分算式分解，以尽量减少内存分配，提高速度。

22. 下一个是应该算是神器级别的函数：frompyfunc，把Python里的函数（可以是自写的）转化成ufunc，用法是frompyfunc(func,
nin, nout)，其中func是需要转换的函数，nin是函数的输入参数的个数，nout是此函数的返回值的个数。注意frompyfunc函数无法保证返回的数据类型都完全一致,因此返回一个中间类型object，需要再次obj.astype(np.float64)之类将其元素类型强制调齐。

23. 作者对broadcast的处理是直接翻成“广播”。补齐的规则除了以前记下的，还有：维度差异的部分通过在已有的维度前面加1进行补齐；输出数组的shape是输入数组shape的各个轴上的最大值；当输入数组的某个轴的长度为1时，沿着此轴运算时都用此轴上的第一组值。

24. repeat方法：b.repeat(6,axis=0)，在第0轴（shape里的第一个维度）的方向上长度扩展为6，也就是第2个维度重复5次（有点绕，呵呵）。

25. ogrid对象：用切片组元作为下标进行存取，返回一到多个可以用来广播计算的数组。有两种用法：开始值:结束值:步长，和np.arange类似（不包含结束值喽）； 开始值:结束值:长度j，当第三个参数为虚数时，它表示返回的数组的长度，和np.linspace类似（包含结束值）。

26. ufunc的reduce 方法：和Python的reduce函数类似，沿着一个指定的轴将运算符插入这根轴上所有子数组或者元素当中，并返回结果，比如：np.add.reduce([[1,2,3],[4,5,6]], axis=1)，其中axis为1，指的是第2个维度上的元素，即[1,2,3]和[4,5,6]。

27. accumulate与reduce类似，不同之处在于返回的数组里包括所有中间结果，致使其shape与输入数组的相同。比如：np.add.accumulate([[1,2,3],[4,5,6]], axis=1)，返回array([[1,3,6],[4,9,15]])，其中1（1），3（1+2），4（4），9（4+5）是中间结果。

28. reduceat方法通过indices参数（列表形式）指定一系列插入reduce的起始和终止位置（reduce at嘛），规则如下：如果indice中某元素小于其后元素，则相应结果为对以这两个元素为位置产生的slice里的数组元素进行reduce；否则结果是这个元素对应的数组元素。对于最后一个元素，因为其后再没元素，结果为对所有元素进行reduce。例子：np.add.reduceat(np.array([1,2,3,4]),indices=[0,1,0,2,0,3,0])，返回array([1,2,3,3,6,4,10])。

29. outer方法，.outer(a,b)方法的计算等同于如下程序：a.shape += (1,)*b.ndim    (a,b)    a = a.squeeze()。例子：np.multiply.outer([1,2,3,4,5],[2,3,4])。在这个例子中，首先，a.shape += (1,)*1，即a.shape变成了(5,1)，然后multiply(a,b)，然后，a = a.squeeze()剔除数组a中长度为1的轴，也就是恢复a的shape，把a给还原了。

30. 终于进入矩阵部分了，NumPy和Matlab不一样，对于多维数组的运算，缺省情况下并不使用矩阵运算。但因为NumPy中同时存在ndarray和matrix类，用户很容易将两者弄混，有违the Zen of Python（import this！），因此matrix类的优先状态是不启用。

31. 矩阵乘法：如果需要将一维数组当作列矢量或者行矢量进行矩阵运算时，推荐先reshape之。

32. dot函数：对于一维数组，计算内积，就是两个数组对应下标元素的乘积（应该就是点积）；对于二维数组计算矩阵乘积（嘛意思？）；对于更高维数组，规则是：dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m])，即结果数组中的每个元素都是数组a的最后一维上的所有元素与数组b的倒数第二维上的所有元素的乘积之和。

33. np.alltrue()，check是否数组每个元素均为True。

34. inner : 对于两个一维数组计算内积；对于多维数组，它计算的结果数组中的每个元素都是：数组a和b的最后一维的内积，因此数组a和b的最后一维的长度必须相同，即：inner(a, b)[i,j,k,m] = sum(a[i,j,:]*b[k,m,:])。

35. outer : 只按照一维数组进行计算，如果传入参数是多维数组，则先将此数组展平为一维数组之后再进行运算，结果是列向量和行向量的矩阵乘积。感觉跟multiply.outer没啥区别啊，见29。

36. np.linalg库还有很多线性代数领域的函数和功能，比如以前用过的np.linalg.det。这本书里提到了可求逆矩阵的inv函数，和求解多元一次方程组的solve函数。后者用法：np.linalg.solve(a,b)，其中a是一个N*N的二维数组，而b是一个长度为N的一维数组，solve函数找到一个长度为N的一维数组x，使得a和x的矩阵乘积正好等于b，数组x就是多元一次方程组的解。

37. 现在进入文件存取了。使用tofile可以方便地将数组中数据以二进制的格式写进文件。tofile输出的数据不带格式，shape和dtype都没有，且使用C语言格式输出，不管数组本身的格式。因此，用fromfile读回来的时候需要自己格式化数据。

38. tofile方法还可以输出为文本格式，此时需要指定sep参数。用fromfile读取时也需要指定此参数。

39. 较为方便的方法是用load和save函数进行数据存取，文件名后缀是.npy，格式是NumPy专用的二进制类型，不用再管shape和dtype这些，不过代价是很难与其他语言的程序协作。通用和专有的矛盾！

40. savez：多个数组保存在同一文件中，第一个参数是输出文件名，后面是其后的参数都是需要保存的数组，自动起名为arr_0, arr_1, ...，也可以使用关键字参数为数组起一个名字。输出文件扩展名为.npz，是一个压缩包，其中每个文件都是.npy，内含一个数组，文件名即内涵数组名。

41. 对于.npz文件，load函数自动识别并且返回一个类似于字典的对象，可以通过数组名作为关键字获取数组的内容。

42. savetxt和loadtxt分别存取txt文件，与tofile不同的是数字此时不是以二进制的方式保存的。用法：np.savetxt("a.txt", a, fmt="%d", delimiter=",")，fmt默认为'%.18e'，delimiter默认为' '。显然如果以非默认设置输出了读取时也要修改相应参数。

43. load和save函数还可以对已经打开的文件对象进行操作，比如：f = file("result.npy", "wb")    np.save(f, a)    f.close()，其中a是一个数组。读取时也是先定义一个句柄f，然后多次np.load(f)，可以顺序（save时的顺序）读取数组。