Python语言程序设计MOOC-10周数学库的使用及其他

2020-06-21

python

在MOOC上，学习了北京理工大学嵩天老师的Python语言程序设计，此为学习笔记记录，备查，备翻阅复习。

numpy库的使用

numpy库

Numeric python 是用python 实现的科学计算包。Numpy 比python自身的嵌套列表要高效得多。使用时，需要先安装numpy库。

该库包括，强大的N维数组对象array，函数库，实用的线性代数、傅里叶变换和随机数生成函数。

Numpy的主要对象是，同种元素的多维数组。

数组的维度dimensions，叫做轴axes
轴的个数叫做秩rank
例如，在3D空间的一个点的坐标[1, 2, 3]是一个秩为1的数组，因为它只有一个轴，轴的长度是3.
例如，在[[1., 0., 0.], [0., 1., 2.]]中，数组的秩是2，第一个维度长度是2，第二个维度长度是3.

Numpy的数组类，被称作是ndarray,通常被称作数组。注意numpy.array 和标准python库类array.array并不相同，后者只处理一维数组和提供少量功能。

ndarray 对象的属性

ndarray.shape 数组的维度，这是一个指示数组在每个维度上大小的整数数组
ndarray.size 数组元素的总个数，等于shape属性中元组元素的乘积。
ndarray.dtype 一个用来描述数组中元素类型的对象。
ndarray.itemsize 数组中每个元素字节的大小。
ndarray.data 包含实际数组元素的缓冲区，通常我们通过索引引用数组元素，所以不使用这个属性。

属性调用示例：

from numpy import *
a = arange(15).reshape(3, 5)
a
Out[4]: 
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14]])

a.shape
Out[5]: (3, 5)

a.ndim
Out[6]: 2

a.dtype.name
Out[7]: 'int32'

a.itemsize
Out[8]: 4

a.size
Out[9]: 15

type(a)
Out[10]: numpy.ndarray

b = array([6, 7, 8])

b
Out[12]: array([6, 7, 8])

type(b)
Out[13]: numpy.ndarray

numpy创建数组

方法1，使用array函数，利用常规的python列表和元组创造数组，所创建的数组类型由原序列中的元素类型决定。
方法2，使用占位符创建数组。例如，zeros, ones, empty。
方法3， numpy提供arange函数来返回一个数组。

# 方法1
from numpy import *
a = array([2, 3, 4]) # 从整型数据列表创建
a
Out[14]: array([2, 3, 4])

a.dtype
Out[15]: dtype('int32')

b = array([1.2, 3.5, 5.1]) # 从浮点数据列表创建

b.dtype
Out[17]: dtype('float64')

c = array([[1, 2], [3, 4]], dtype = complex) # 从嵌套列表创建，并制定array数据类型
c
Out[18]: 
array([[1.+0.j, 2.+0.j],
       [3.+0.j, 4.+0.j]])
# 常见的一个错误是，提供的不是列表，而是列表的元素
a = array(2, 3, 4)
a
Out[18]: ValueError: only 2 non-keyword arguments accepted

# 方法2
zeros((3, 4)) # 创建一个全是0的数组
Out[20]: 
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])

ones((2, 3, 4), dtype = int16) # 创建一个全是1的数组，dtype 可以指定
Out[21]: 
array([[[1, 1, 1, 1],
        [1, 1, 1, 1],
        [1, 1, 1, 1]],

       [[1, 1, 1, 1],
        [1, 1, 1, 1],
        [1, 1, 1, 1]]], dtype=int16)

empty((2 ,3)) # 创建一个内容随机并且依赖于内存状态的数组
Out[22]: 
array([[2.32e-321, 0.00e+000, 0.00e+000],
       [0.00e+000, 0.00e+000, 0.00e+000]])

# 方法3
arange(10, 30, 5) # 即返回一个数组
Out[23]: array([10, 15, 20, 25])

arange(0, 2, 0.3) # arange也接受浮点型的参数
Out[24]: array([0. , 0.3, 0.6, 0.9, 1.2, 1.5, 1.8])

打印数组

打印数组时，numpy以类似嵌套列表的形式显示
一维数组被打印成行
二维数组被打印成矩阵
三维数组被打印成矩阵列表

a = arange(6) # 1d array
print(a)
[0 1 2 3 4 5]

b = arange(12).reshape(4, 3) # 2d array
print(b)
[[ 0  1  2]
 [ 3  4  5]
 [ 6  7  8]
 [ 9 10 11]]

c = arange(24).reshape(2, 3, 4) # 3d array
print(c)
[[[ 0  1  2  3]
  [ 4  5  6  7]
  [ 8  9 10 11]]

 [[12 13 14 15]
  [16 17 18 19]
  [20 21 22 23]]]

基本运算

数组的算术运算是按元素进行。Numpy中的乘法运算符*是指示按元素计算
矩阵乘法，可以使用dot函数或创建矩阵对象实现

a = array([20, 30, 40, 50])
b = arange(4)
b
Out[30]: array([0, 1, 2, 3])

c = a - b
c
Out[32]: array([20, 29, 38, 47])

b ** 2
Out[33]: array([0, 1, 4, 9], dtype=int32)

10 * sin(a)
Out[34]: array([ 9.12945251, -9.88031624,  7.4511316 , -2.62374854])

a < 35
Out[35]: array([ True,  True, False, False])

非数组的运算，可以利用ndarray类方法实现。
numpy提供常见的通用函数，如sin, cos, exp. 在numpy中这些函数作用按数组的元素运算，产生一个数组作为输出。

a = random.random((2, 3))
a
Out[36]: 
array([[0.62669257, 0.62536583, 0.69025192],
       [0.54679281, 0.29197729, 0.14721231]])

a.sum()
Out[37]: 2.928292721634076

a.min()
Out[38]: 0.1472123098871817

a.max()
Out[39]: 0.6902519182474993

b = arange(12).reshape(3, 4)
b
Out[40]: 
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

b.sum(axis = 0) # 每一列的求和
Out[41]: array([12, 15, 18, 21])

b.min(axis = 1) # 每一行的最小值
Out[42]: array([0, 4, 8])

b.cumsum(axis = 1) # 每一行的累计和
Out[43]: 
array([[ 0,  1,  3,  6],
       [ 4,  9, 15, 22],
       [ 8, 17, 27, 38]], dtype=int32)

索引切片迭代

数组可以被索引切片和迭代。

a = arange(10) ** 3
a
Out[53]: array([  0,   1,   8,  27,  64, 125, 216, 343, 512, 729], dtype=int32)

a[2]
Out[54]: 8

a[2:5]
Out[55]: array([ 8, 27, 64], dtype=int32)

a[0:6:2] = -1000 # 等价于a[:6:2] = -1000，即从pos 0到6，每隔1个修改元素为-1000

a
Out[57]: 
array([-1000,     1, -1000,    27, -1000,   125,   216,   343,   512,
         729], dtype=int32)

a[::-1] # reverse a
Out[58]: 
array([  729,   512,   343,   216,   125, -1000,    27, -1000,     1,
       -1000], dtype=int32)

矩阵运算

numpy对于多维数组的运算，缺省情况下并不适用矩阵运算，对数组进行矩阵运算，可调用相应的函数。
numpy库也提供了matrix类，使用matrix类创建的是矩阵对象，它们的加减乘除运算采用矩阵方式计算。

a = matrix([[1, 2, 3], [5, 5, 6], [7, 9, 9]])

a * a ** -1
Out[62]: 
matrix([[ 1.00000000e+00,  1.66533454e-16, -1.11022302e-16],
        [ 0.00000000e+00,  1.00000000e+00, -4.44089210e-16],
        [ 4.44089210e-16,  5.55111512e-17,  1.00000000e+00]])

矩阵中更高级的运算，可以在numpy的线性代数子库linalg中找到，
例如， inv函数，计算逆矩阵
solve函数，求解多元一次方程组

卫星定位（GPS）程序的numpy实例

GPS定位的基本原理，是根据高速运动卫星的瞬间位置作为已知的起算数据。采用空间距离-后方交会的方法，确定待测点的位置。

假设t时刻在地点待测点上安置GPS接收机，可以测定GPS信号到达接收机的时间delta t，再加上接收机所接收到的卫星星历等其他数据，就可以确定一个方程组来对位置信息进行求解。

举例：

设地球上一个点R ，同时收到66颗卫星（S1， S2， S3， S4， S5， S6）发射的信号：

信号源	x y z位置	收到信号的时间戳
S1	3，2，3	10010.00692286
S2	1，3，1	10013.34256381
S3	5，7，4	10016.67820476
S4	1，7，3，	10020.01384571
S5	7，6，7	10023.34948666
S6	1，4，9	10030.02076857

设 (x, y, z)是R的当前位置，t 是发出信号的时间戳，信号是以光速传输的。

则对于信号源S1满足， (x - 3) ^2 + (y - 2) ^2 + (z - 3)^2 = [(10010.00692286 - t) * c] ^2

即欧式空间距离与信号传输距离相等。

对于其他5个信号源，也可以得到类似的方程，进而建立一个方程组。虽然对于每个方程都是，一元二次方程，但是因为二次项的系数都是1，所以容易消去二次项。从而得到一元一次方程组。

\[\left(\begin{array}{cccc} 4 & -4 & -12 & 3.59751 \\ 0 & -2 & -16 & 2.99792 \\ 8 & 6 & -10 & 2.39834 \\ 0 & 6 & -12 & 1.79875 \\ 12 & 4 & -4 & 1.19917 \end{array}\right)\left(\begin{array}{c} x \\ y \\ z \\ t \end{array}\right)=\begin{array}{|c|} 35971.1 \\ 29957.2 \\ 24031.4 \\ 17993.5 \\ 12059.7 \end{array}\]

即 A * X = B

所以X = A的逆 * B

所以上述过程即可以通过python来求解完成。

GPS定位计算问题的IPO模式

输入， 6颗卫星的欧式坐标和信号
处理，GPS定位算法
输出，GPS接收设备的地理坐标和当前时间

算法就是：假设第i颗卫星的坐标和时间戳表示为 (xi, yi, zi, ti)

\[\left(\begin{array}{llll} x_{1}-x_{2} & y_{1}-y_{2} & z_{1}-z_{2} & c^{2}\left(t_{1}-t_{2}\right) \\ x_{1}-x_{3} & y_{1}-y_{3} & z_{1}-z_{3} & c^{2}\left(t_{1}-t_{3}\right) \\ x_{1}-x_{4} & y_{1}-y_{4} & z_{1}-z_{4} & c^{2}\left(t_{1}-t_{4}\right) \\ x_{1}-x_{5} & y_{1}-y_{5} & z_{1}-z_{5} & c^{2}\left(t_{1}-t_{5}\right) \\ x_{1}-x_{6} & y_{1}-y_{6} & z_{1}-z_{6} & c^{2}\left(t_{1}-t_{6}\right) \end{array}\right)\left(\begin{array}{l} x \\ y \\ z \\ t \end{array}\right)=\left(\begin{array}{l} {\left[\left(t_{2}^{2}-t_{1}^{2}\right) c^{2}-\left(x_{2}^{2}-x_{1}^{2}+y_{2}^{2}-y_{1}^{2}+z_{2}^{2}-z_{1}^{2}\right)\right] / 2} \\ {\left[\left(t_{3}^{2}-t_{1}^{2}\right) c^{2}-\left(x_{3}^{2}-x_{1}^{2}+y_{3}^{2}-y_{1}^{2}+z_{3}^{2}-z_{1}^{2}\right)\right] / 2} \\ {\left[\left(t_{4}^{2}-t_{1}^{2}\right) c^{2}-\left(x_{4}^{2}-x_{1}^{2}+y_{4}^{2}-y_{1}^{2}+z_{4}^{2}-z_{1}^{2}\right)\right] / 2} \\ {\left[\left(t_{5}^{2}-t_{1}^{2}\right) c^{2}-\left(x_{5}^{2}-x_{1}^{2}+y_{5}^{2}-y_{1}^{2}+z_{5}^{2}-z_{1}^{2}\right)\right] / 2} \\ {\left[\left(t_{6}^{2}-t_{1}^{2}\right) c^{2}-\left(x_{6}^{2}-x_{1}^{2}+y_{6}^{2}-y_{1}^{2}+z_{6}^{2}-z_{1}^{2}\right)\right] / 2} \end{array}\right.\]

程序中会使用的numpy 的函数包括：

numpy.dot(a, b)计算矩阵a 和矩阵 b 的点积

numpy.linalg.inv(a) 求矩阵a 的逆矩阵

程序如下：

from numpy import *
i = 1
c = 0.299792458 # 光速，单位km/us
x = zeros((6, 4)) # 存储6个卫星的（x, y, z, t）参数
while i <= 6:
    print("%s %d" % ("please input (x, y, z, t) of group", i))
    temp = input()
    x[i - 1] = temp.split()
    j = 0
    while j < 4:
        x[i - 1][j] = float(x[i - 1][j])
        j += 1
    i += 1
a = zeros((4, 4)) # 系数矩阵
b = zeros((4, 1)) # 常数项
j = 0
while j <4:
    a[j][0] = 2 * (x[5][0] - x[j][0])
    a[j][1] = 2 * (x[5][1] - x[j][1])
    a[j][2] = 2 * (x[5][2] - x[j][2])
    a[j][3] = 2 * c * c *(x[j][3] - x[5][3])
    b[j][0] = x[5][0] * x[5][0] - x[j][0] * x[j][0] + \
              x[5][1] * x[5][1] - x[j][1] * x[j][1] + \
              x[5][2] * x[5][2] - x[j][2] * x[j][2] + \
              c * c * (x[j][3] * x[j][3] - x[5][3] * x[5][3])
    j += 1
a = linalg.inv(a) # 系数矩阵求逆
print(dot(a, b))

运行的时候，就是先提示输入6个点

please input (x, y, z, t) of group 1

3 2 3 10010.00692286
please input (x, y, z, t) of group 2

1 3 1 10013.34256381
please input (x, y, z, t) of group 3

5 7 4 10016.67820476
please input (x, y, z, t) of group 4

1 7 3 10020.01384571
please input (x, y, z, t) of group 5

7 6 7 10023.34948666
please input (x, y, z, t) of group 6

1 4 9 10030.02076857

再往下运行就是输出待测点的坐标和时间戳了：

print(dot(a, b))
[[5.e+00]
 [3.e+00]
 [1.e+00]
 [1.e+04]]

Matplotlib库

一个数据可视化的函数库。使用前需要安装。其官方手册达2600页，这里介绍主要功能。

最常用的就是子库，pyplot。

提供2D图标制作的基本函数
在图形上创建画图区域
在画图区域上画线
在线上标注等功能

举例：简单散点图的绘制

import matplotlib.pyplot as plt
plt.plot([1, 2, 3, 4, 5, 6], 'ro') # 只有一行数据时，默认是y值，并自动x从0赋值。
# ro 表示点的颜色为红色，形状为圆形
plt.ylabel("Description of y value") # ylabel 是设置纵坐标轴名称
plt.show() # 图形参数设置好之后，show 显示图形

Python语言程序设计MOOC-10周 数学库的使用及其他