可视化计算——数组计算与曲线绘图（上）

(//▲//)大计第三次笔记bingo!

可视化计算——数组计算与曲线绘图（上）

①如何用代码绘制图形——mathplotlib

( ightarrow)mathplotlib是python最著名的绘图库，文档相当完备，提供了一整套和matlab相似的命令API.另外，Gallery页面中有几百幅缩略图，打开之后都有源程序，如果需要绘制某种类型的图，只需要在http://matplotlib.org/gallery.html中浏览/复制/粘贴一下，基本上都能搞定

现在用最简单的坐标点方式绘制sin(x)图像：

import matplotlib.pyplot as plt #导入绘图库并命名为plt
from math import sin,radians    #redians是输入角度转化为弧度的函数，sin运算值应为弧度值
x=[0,30,60,90,120,150,180,210,240,270,300,330,360]
y=[]
for i in x:
    m=sin(radians(i))
    y.append(m)                 #y.append(m)函数将数m添加到列表y的最后
plt.plot(x,y,'.')               #绘图
plt.show()                      #打印结果

可以得到结果：

我们只需要将相邻x的间隔缩小，就可以得到一个近似连续的函数图形：

import matplotlib.pyplot as plt 
from math import sin,radians
x=[]
y=[]
for i in range(0,360+1,1):                #取步长为1，范围在[0,361)
    x.append(i)
y=[sin(radians(m))for m in x ]            #遍历x中的所有元素
y=[sin(radians(n))for n in range(len(x))] #遍历所有下标
plt.plot(x,y,'.')
plt.show()

其中有两种方法取完0，1，...，360的每一个值，一个方法是遍历x列表中所有元素，另一个是利用函数len(x)计算出列表x里面的元素个数有361个，然后用range(361)取完从0到360的每一个数.相当于给x中每个元素编号i0，i1，i2，...，i360，然后遍历所有的下标.

然而对于更一般的函数，可以函数抽象，封装细节，让其更具通用性，也就是运用def定义函数.例如我们得到的任意一个表达式：

[f(x)=sin(frac{x}{5})+frac{sin(frac{x}{2})×cos(frac{x}{3})}{1.2+sin(frac{x}{2})} ]

我们可以先将它定义，然后绘图：

import matplotlib.pyplot as plt 
from math import sin,radians,cos
def f(x):
    result=sin(x/5)+sin(x/2)*cos(x/3)/(1.2+sin(x/2))
    return result
x=[]
for i in range(361):
    x.append(i)
y=[f(radians(e))for e in x]
plt.plot(x,y,'.')
plt.show()

得到如下结果：

如此我们就能对一般的函数进行粗略的作图了

②完善绘图参数

我们可以通过更改plot.()里面的参数来更好地达到绘图的需求

修改图线颜色的参数：plt.plot(x,y,'b.')
其中b代表blue为蓝色，其余还有y-黄色；r-红色；g-绿色；c-蓝绿色；m-洋红色；k-黑色；w-白色

修改图标形状的参数：plt.plot(x,y,'y-')
其中'-'代表的是直线，其余还有：

另外，我们还可以加上其他的一些元素：

plt.grid('on')           #打开网格
plt.xlabel('x')          #在x轴上添加标记
plt.ylabel('y')          #在y轴上添加标记
plt.axis([0,360,-1,1])   #设置坐标的取值范围
plt.legend(['sin(x)','cos(x)'],loc='upper center')   #添加图例，并且上部居中

例如，我们绘制角度在[0,360)的sin(x),cos(x)及tan(x)的图像：

import matplotlib.pyplot as plt 
from math import sin,radians,cos,tan
x=[]
for i in range(360):
    x.append(i)
y1=[sin(radians(p))for p in x]
y2=[cos(radians(q))for q in x]
plt.plot(x,y1,'y-')
plt.plot(x,y2,'b--')
x.remove(90)
x.remove(270)
y3=[tan(radians(m))for m in x]
plt.plot(x,y3,'r-.')
plt.axis([0,360,-2,2])
plt.grid('on')
plt.legend(['sin(x)','cos(x)','tan(x)'],loc='upper center')
plt.show()

可以得到图形：

在这个例子中我们需要非常注意的是，如果不限制y轴范围，则图形会因为x趋近于90度时tan(x)趋近于无穷而变得很大，展示出来的sin和cos就被压缩的很小（几乎变成了一条直线）；除此之外，因为tan(x)的定义域不包括90度和270度两个点，所以在画tan(x)的图像的时候需要用x.remove()函数清除掉列表x中的90和270两个值.

③数组计算----Numpy

在上面的例子中，我们运用range(0,360,1)取步长为1绘制坐标点，但如果我们需要画出[0,1)之间的图像，或者需要更加精确（有更多的坐标点）地描绘出图像，我们的步长就必须要缩短成浮点数.但是经过实践发现range()不支持浮点数类型的步长，这就必须引进numpy里面的arange函数，引入数组的概念.

数组（Array）是有序的元素序列。若将有限个类型相同的变量的集合命名，那么这个名称为数组名。组成数组的各个变量称为数组的分量，也称为数组的元素，有时也称为下标变量。用于区分数组的各个元素的数字编号称为下标。数组是在程序设计中，为了处理方便， 把具有相同类型的若干元素按有序的形式组织起来的一种形式。这些有序排列的同类数据元素的集合称为数组。
数组是用于储存多个相同类型数据的集合。
——————百度百科

list：可以作数组使用
元素可以是任何对象，因此，保存的是对象的指针
保存简单的[1,2,3]，需要有3个指针和3个整数对象，比较浪费内存和CPU计算时间

array对象：直接保存数值（不是指针）array(‘i’, [1,2,3])和C语言的一维数组类似
不支持多维，也没有运算支持，不适合数值处理

numpy：针对多维数组的扩展包，提供两种基本对象：ndarray（N-dimensional array object）存储单一数据类型的多维数组；ufunc（universal function object）对数组进行处理的函数
如：

import matplotlib.pyplot as plt
from numpy import array, sin, radians
x = array(range(-5*360, 5*360))
y = sin(radians(x))
plt.plot(x, y, 'b-')
plt.show()

其中我们用array产生了ndarray类型数组
用y=sin(radians(x))代替了列表y=[sin(radians(i))for i in x]

生成ndarray数组：

①arange(开始值,终值，步长)
ex:arange(0,360,0.1)
②linspace(开始值，终值，指定个数)
ex:linspace(0,360,3601)（即从0到360取相邻两数间隔相等的3601个数，包括0与360）

转换成ndarray数组：

①array()
将list或者tuple类型转换为ndarray类型
一维：array([1,2,3])
二维：array([ [1,2,3],[4,5,6] ])
②reshape()
改变数组形状
array([ [1,2,3],[4,5,6] ]).reshape(3,2)
(reshape(x,y)即将数组内元素变成x个元素中有y个元素)

例如：

x=array([1,2],[3,4],[5,6],[7,8]).reshape(2,4)
print(x)

PS C:python> python .	emp.py
[[1 2 3 4]
 [5 6 7 8]]

ndarray中一些特殊数组的生成：

①全0数组zeros((2,3))，结果是 [ [ 0. 0. 0.][ 0. 0. 0.] ]
②全1数组ones(3)，结果是[ 1. 1. 1.]
③单位方阵identity(3)，结果是3阶单位方阵
④对角矩阵eye(3)，结果是 [ [ 1.  0.  0.][ 0.  1.  0.][ 0.  0.  1.] ]
eye(N,M=None, k=0, dtype=<type 'float'>)
第k条对角线全1，其它元素为0的矩阵
(ps:从左上角到右下角是标号为0的对角线，向右上角移动一次标号+1，向左下角移动一次标号-1)

ndarray数组的属性：

例如a=zeros((2,2,2))
①a.ndim----表示数组维数，结果为3
②a.shape----表示数组每一维的大小，结果为(2,2,2)
③a.size----表示数组元素总个数，结果为8
④a.dtype----表示元素的类型，结果为float64
⑤a.itemsize----表示每个元素所占字节数，结果为8

ndarray数组的简单索引与切片：

例如a=array([ [2,3,4],[5,6,7] ])
①a[1,2]----表示的是索引第1行第2个数（特别注意行和列都是从0开始算的），结果为7
②a[1,0:2]----表示的是索引第一行整行从标号为0到标号为1的数（取不到2，左闭右开），结果为[5,6]
③a[1,:] = [8,9,10]----表示将第一行整个换成[8,9,10]，print(a)的结果为[ [2,3,4],[8,9,10] ]

ndarray数组的加减乘除计算：

例如a = ones((2,2))b = eye(2)
a×b----每个位置的数相应相乘
a±b----每个位置的数相应加减
a/b----每个位置的数相应相除（ps：若其中某项出现除数为0，会报错）
a×2----将每个位置的数扩大两倍
a.sum()----计算数组里所有元素的和
a.max()----找到数组里所有元素中的最大值

数组的拷贝与深度拷贝：

尝试一下两组代码：

from numpy import *
a=array([[1,2],[3,4]])
b=a
b[0,1]=3
print(a,b,sep='
')

from numpy import *
a=array([[1,2],[3,4]])
c=a.copy()
c[0,1]=3
print(a,c,sep='
')

发现第一组用赋值的方法得来的数组b变换后数组a也跟着进行了变换，第二组用函数copy得来的数组c变换后并没有影响数组a.可以简单地解释为b由a赋值得来但都具有相同的指针，而c和a具有不同的指针.因为指针还没有深度学习所以暂时没有展开讲，相信学习以后能够更好地解释.

利用数组更加详实地绘制图像：

我们绘制步长更小，描点更多的正弦函数图像，用到数组相关的知识：

import numpy as np
from math import radians
import matplotlib.pyplot as plt
x=np.arange(-360,360,0.001)       #利用数组可以取到小数的步长
radians_vec=np.vectorize(radians) #将函数radians利用vectorize向量化
y=np.sin(radians_vec(x))
plt.plot(x,y,'b-')
plt.grid('on')
plt.show()

就可以得到相应的函数图像.

④计算效率

再上一个例子中我们用到了函数矢量化的方法，是因为numpy提供了向量化函数，支持数组的向量运算，而普通的列表仅能够标量计算.我们可以通过time函数来比较标量计算和矢量计算的时间长短.

import time
from math import *
start=time.clock()            #第一次调用，返回程序运行的实际时间
x=[]
for i in range(0,36000000):
    x.append(i)
y=[sin(e)for e in x]
t=time.clock()-start          #第二次调用，返回自第一次调用后，到这次调用的时间间隔
print('math.sin标量计算时间为%.5f'%t)

PS C:python> python .	emp.py
math.sin标量计算时间为7.91506

import time
import numpy as np
start=time.clock()
x=np.arange(0,36000000)
y=np.sin(x)
t=time.clock()-start
print('numpy.sin矢量计算时间为%.5f'%t)

PS C:python> python .	emp.py
numpy.sin矢量计算时间为0.42856

可以看出矢量计算比标量计算的时间快了10倍以上，这是因为numpy矢量计算的是C语言内部实现.但是numpy进行标量计算时的运行速度要比math标量计算慢很多：

mport time
import numpy as np
import math

x=np.linspace(0, 10, 1000001)

start = time.clock()
y1 = np.sin(x)
t1 = time.clock()-start
print ('numpy.sin矢量计算耗时：%.4f'%t1, '秒')

start = time.clock()
y2 = [math.sin(e) for e in x]
t2 = time.clock()-start
print (' math.sin标量计算耗时：%.4f'%t2, '秒')

start = time.clock()
y3 = [np.sin(e) for e in x]
t3 = time.clock()-start
print ('numpy.sin标量计算耗时：%.4f'%t3, '秒')

PS C:python> python .	emp.py
numpy.sin矢量计算耗时：0.0014 秒
 math.sin标量计算耗时：0.1709 秒
numpy.sin标量计算耗时：0.8727 秒

⑤函数向量化方法

在“利用数组更加详实地绘制图像”中提到了vectorize向量化函数.其实向量化函数分为“自动向量化”和“手动向量化”两种方法.

在自动向量化中，如果没有if检测，则直接将from math import *改为from numpy import *;如果有if检测，则需要使用vectorize()进行自动向量化

在手动向量化中，可以在使用where()替换if检测；也可以使用循环专门处理数组参数

例如：

from math import log
from numpy import linspace,vectorize
import matplotlib.pyplot as plt
log_vec=vectorize(log)       #将math库里的log函数向量化
def H(x):
    if x>0:
        m=log_vec(x)
        return m
    elif x<0:
        p=log_vec(-x)
        return p
x=linspace(-5,5,100)        
H_vec=vectorize(H)           #将自己定义的函数H(x)向量化
y=H_vec(x)
plt.plot(x,y,'b-')
plt.grid('on')
plt.show()

得到下面的函数图像：

当然有一些简化代码的方法，比如：

def f(x):
    return (0 if x>=0 else 1)   #避免写两遍的return

def f(x):
    return where(x>=0,0,1)      #在x>=0时输出0，否则输出1
#其中用where的好处是此定义函数可以不用vectorize来向量化
#另外还有些简单的向量化写法：
y=vectorize(f)(x)  #在向量化函数的同时计算了函数值，避免写两行代码

另外还可以运用循环专门处理数组参数：

#便于理解，举例一个数组x=[1,2,3,4,5]
def H_loop(x):
    r = zeros(len(x))   #产生一个长度为5的全0数组r=[0,0,0,0,0]
    for i in range(len(x)):   #对于range(5),即[0,1,2,3,4]的列表
        r[i] = H(x[i])   #r[0]=H(x[0])=H(1),......
    return r

len()是一个常用的创建列表或者数组对应下标的方法.

这次笔记就记到这里啦！ENDo(￣▽￣)ブ