Python科学计算

$$\mathrm{《Python科学计算》学习笔记}$$

[TOC]

Numpy

数组入门

数组创建

np.array：数组
a.shape =：原地变换.
b = a.reshape((...))：变换后，内容共享！
a.dtype
np.arange(始,末不含,步长)
np.linspace(),[endpoint=False]
logspace([base = x])
a = np.fromstring(s,dtype = int)
a = np.fromfunction(func,(10,))
对于array的切片赋值，内容共享（这和list列表不同）

下标方法

切片赋值，内容共享（这和list列表不同）
整数列表存取，不共享
bool数组，不共享
- 不能用bool列表，否则等同于0，1的整数列表存取

多维数组

二维数组：a = b(列数组) + c(行数组)
数组下标是一个元组！（不是slice）

多维数组存取

特例：整数,整数：返回得到内容数字，不是array！
整数看成切片，只有整数/切片，内容共享（返回array）
- slice(a,b,c)返回切片
- np.s_[::2,2:]生成切片
只要有元组/列表，bool数组任何一个，内容不共享
- 如果元组/列表只有一个数字，当作数字->切片处理，但依然不共享
- 元组/列表非单个数字部分必须所有维数长度一致，用于一一对应（不同于切片作用）
剩余维数对应下标:

结构数组

persontype = np.dtype({'names':['name1','name2',...], 'formats':['S32','i']}, [align = True]) ←内存对齐
np.dtype([('name1','|S32'),('name2','i'),...])
np.dtype({'name1':(类型,地址偏移),...})

a = np.array([("HKK",3),("KHH",4)],dtype = persontype)
- |忽视顺序；<低位字节在前；>高位字节在前
a[0]['name1']
a.tofile('test.dat')写入二进制文件
b = a['name1']内容共享
元素为一个数组的定义方法2：('name','type',(维度长度))

数组的存储细节

stride各维地址偏移量，dimensions各维长度，dim count维数
np.array([order="F"])：数组按Fortan方式存储，即列优先
as_strided(a,shape = , strides = )
查询
- a.flags的OWNDATA，C_CONTIGUOUS，F_CONTIGUOUS
- id(b.base) == id(a)

函数

ufunc

np.sin(x,[out = x])
math.sin()（单独运算快）和np.sin()（整个数组一起运算快）
a.item(*,*) ~ a[*,*]但返回标准Python类型

数学运算

y = x1 + x2 ->np.add(x1,x2[,y])
subtract，multiply
divide（/），true_divide（/返回精确），floor_divide（//返回取整）
- _future_.division？
negative（-），power（**）
remainder或mod
表达式不能太复杂：自动的中间变量

关系运算

equal，not_equal，less，less_equal，greater，greater_equal
均已重载

布尔运算

无法重载原符号
np.logical_and(a==b , a>b) 相当于a>=b
np.logical_or，np.logical_not，np.logical_xor

其它

位运算：np.bitwise_and, ...相当于&,|,~,^
- 对布尔数组位运算=关系运算！（但注意优先级变了）
any()，all()

自定义ufunc函数

一般函数：somefun(x,c1)
可对数组操作函数：ufunc = np.frompyfun(somefun,2,1) （输入个数，输出个数）
y = ufunc(x,c1)（x可以是数组）
y = y.astype(np.float)（必须改类型）（astype产生新数组）

广播

维度数目不同，少的在低维添长为1的维
对应维运算，某维长为1的运算当作复制为需要长度；长度不为1的，如果长度不同，报错
新产生广播（结果有用维1维）：x,y = np.ogrid[:5,:5]
- np.ogrid[:1:4j,:1:3j]表示0-1中4个值/3个值（相当linspace）
（结果有用维多维）：x,y = np.mgrid[:5,:5]
已有数组产生广播（结果有用维1维）：gy,gx = np.ix_(y,x)

ufunc函数方法

<operation>.reduce(array [,axis = 0, dtype = None])
- 数组元素连续运算，只给最后结果
<op>.accumulate(...)：类似reduce但是，结果与输入维数相同，给中间结果
<op>.reduceat(array, indices = [...])：结果和indices长度相同
- 前小后大连续算；前大后小给那个值；最后一个连续算到[-1]
<op>.out(array1, array2)：产生二维数组，是一个组合运算结果表

多维数组下标存取

下标是元组，应尽量显式使用元组表示下标
下标：整数；切片；整数数组（列表/元组则转）；布尔数组
只有整数数组：
- 各维数组广播，取对应维的值，结果与广播后维度、长度相同
- ind0, ind1, ind2 = np.broadcast_arrays(i0,i1,i2)与b = a[i0,i1,i2]
整数数组+切片
- 切片的那一维造成的维度全部往后放，广播造成的维度靠前
- 整数、单个元素的整数数组造成的效果与切片[..., 2:3,...]效果不同，不产生任何结果维度上的影响，只是该维下标钉死
布尔数组下标
- 相当于使用np.nonzero(b)（b有几维，结果有几个数组组成的元组，表示True元素下标的对应维的值）

函数库

基本函数

求和：np.sum(a [,axis = None, dtype = ***, out = ])
- axis不填则为整个求和，dtype数据较多时要用精度高的
平均值：np.mean(a [,axis = None, dtype = ***, out = ])
- 整数默认用双精度浮点数算
std()：标准差；var()：方差

最值与排序

产生随机数组：a = np.random.randint(0,10,size=(4,5))
min() max() ptp()（极差）
- axis out参数
np.argmax() np.argmin() 第一次出现的位置
- 不指定axis参数，给平坦化后位置
- 平坦化：a.ravel()（不改变a）
- 从平坦化位置转成多维下标：idx = np.unravel_index(num, a.shape)
从小到大排序：np.sort(a [, axis = -1])
- 返回新数组
- 平坦化排序：axis = None
排序下标：idx = np.argsort(a [, axis = -1])
- 返回值表示：原来在该值为下标的那个数字放到现在这个值的位置上，见下
- a[idx] == np.sort(a)
中位数：np.median(a [,axis = -1])

多项式

一元多项式：p = np.poly1d(a)，a是数组，高次项在左边
- p()像函数一样用
p + [1,2] p*p p/[1,1]（返回商式和余式）
求导：p.deriv()；积分：p.integ()
求根：r = np.roots(p)：返回多个根放在一个数组里，p(r) == array([0,0,0,...])
从根转回多项式系数：np.poly(r)，返回系数数组
拟合：a = np.polyfit(x,y,deg)（deg最高次指数），返回系数数组
- 通过系数数组计算多项式值：np.polyval(a,x)

分段函数

x = where(condition, y, z)
- x condition y z是同样维数数组
- C语言级别，比frompyfunc()快
select(condlist, choicelist, default = 0)
- 找condlist中最先满足的
- np.select([x>=c, x<c0], [0, x/c0*hc], default = (c-x)/(c-c0)*hc)
- 内存开销很大
piecewise(x, condlist, funclist)
- funclist比condlist长1，最后表示其它情况
lambda x : x**x：创建简单小函数

统计函数

找出数组中不同值，从小到大排：x = np.unique(a)
- return_index = True：多返回idx,a[idx] == x
- return_inverse = True：返回a中元素在x中的下标
np.bincount(a [, weights = w]) 非负整数数组，元素出现次数统计
- w权重，和a长度相同
直方图统计：hist, bin_edges = np.histogram(a, bins = , range = None, normed = False, weights = None)
- bins区间数，或给定的每个区间边界列表（长度区间数+1）
- range数据范围，默认(最小,最大)
- normed：False统计个数，否则概率化
- hist：统计结果数组；bin_edges：区间边界数组（长度区间数+1）

线性代数

Matrix对象

a = np.matrix([[...],[...],[...]])
直接计算：a* a**-1
后面不再用Matrix对象

dot乘积运算

dot(x,y)
一维乘一维：相当于行列相乘
请不要用一维乘二维：x一维，y二维行向量不行；x二维列，y一维不行；x二维行向量，y一维居然可以
多维乘多维x[i1,i2,l,m], y[j1,j2,j3,m,n]，结果z[i1,i2,l,j1,j2,j3,n]

inner乘积inner(x,y)

多维数组最后一维相乘求和（本质就是一维乘一维，其它维不过相当于列表）
x[i1,i2,m], y[j1,j2,j3,m]，结果z[i1,i2,j1,j2,j3]

outer乘积outer(x,y)

列向量乘行向量！展成二维
输入多维，先展平为一维

多元一次方程组

x = np.linalg.solve(a,b)，a是N*N二维数组，b是N长一维数组，输出N的一维数组
x = np.linalg.lstsq(a,b)，a是X*N数组，得最小二乘解
正态随机数：np.random.standard_normal(len)

文件存取

底层二进制存取

a.tofile("name.bin")
b = np.fromfile("name.bin",dtype = np.int32)
需要指定dtype，读入后b.shape = (*,*)
指定sep参数，文本输出/输入，sep是分隔符

高级二进制存取

np.save("a.npy", a)（单）
- C语言格式
c = np.load("a.npy")（单）
np.savez("result.npz", a, b, newname = c)（多）
- 关键字参数给参数起名，非关键字默认：arr_0,arr_1,...
r = np.load("result.npz")（多）
- r["arr_0"] = r["newname"] =

高级文本存取

np.savetxt("a.txt", a [, fmt = '%.18e', delimiter = " "])
np.loadtxt("b.csv", delimiter = ",", dtype = np.***)
结构存取dtype = persontype

简单的底层文本存取

f = file("test.csv")
f.readline()
data = np.loadtxt(f, delimiter = ",")
f = file("result.npy", "wb") f = file("result.npy", "rb")
np.save(f,a) np.load(f)
- save load可连续使用，在文件中存取多个变量

内存映射数组

从文件读取的变量和文件内容挂钩，改变量也可以改文件
a = np.memmap(filename, deype = uint8, mode = "**", offset = 0, shape = None, order = "C")
- dtpyeNumpy类型；offset读文件起始位置偏移，字节单位；orderC/Fortran
- mode：r只读；c修改但不写入文件；r+数组可读写，结果写入文件；w+创建或覆盖已有文件（可以缺省filename）
a.flush()把a的值按mode规定的规则写入（有的就不能写入）