Python科学计算学习之高级数组（二）

1. 代码性能和向量化

  背景：python是一种解释型的编程语言，基本的python代码不需要任何中间编译过程来得到机器代码，而是直接执行。而对于c、c++等编译性语言就需要在执行代码前将其编译为机器指令。 但是，解释型代码的速度比编译型代码要慢，为了使得python代码更快，最好尽可能的使用numpy和scipy包中的函数编写部分代码。（注意：numpy和scipy是诸如c、c++等编译型语言编写实现的）

例如：python语言的numpy向量化语句为什么比for快？

python之类语言的for循环，和其它语言相比，额外付出了什么。

python是解释执行的。举例来说，执行 x = 1234+5678 ，对编译型语言，是从内存读入两个short int到寄存器，然后读入加法指令，通知cpu内部的加法器动作，最后把加法器输出存储到x对应的内存单元（实质上，最后这个动作几乎总会被自动优化为“把加法器输出暂存到寄存器而不是内存单元，因为访问内存的时间消耗常常是访问寄存器的几十倍”）。一共2~4条指令（视不同cpu指令集而定）。

而换了解释性语言，它得先把“x = 1234+5678”当成字符串，逐个字符比对以分析语法结构——不计空格这也是11个字符，至少要做11个循环；每个循环至少需要执行的指令有：取数据（如读'x'这个字符）、比较数据、根据比较结果跳转（可能还得跳转回来）、累加循环计数器、检查循环计数器是否到达终值、根据比较结果跳转。这就是至少6条指令，其中包含一次内存读取、至少两次分支指令（现代cpu有分支预测，若命中无额外消耗，否则……）。总计66条指令，比编译型语言慢至少17倍（假设每条指令执行时间相同。但事实上，访存/跳转类指令消耗的时间常常是加法指令的十倍甚至百倍）。这还只是读入源码的消耗，尚未计入“语法分析”这个大头；加上后，起码指令数多数百倍（消耗时间嘛……我猜起码得多数千倍吧）。

1. 向量化：

     为提升代码的性能（运行时间），通常需要将代码向量化。使numpy包的切片、运算符和函数来替代代码中的for循环以及运行速度较慢的代码片段，可以显著提高代码的性能。

规则：尽可能避免使用for循环而采用向量化形式，善用python的numpy库中的内置函数。例如：np.exp ,np.log ,np.maxmum(v,0) 等。

简单实例进行说明：

import numpy as np

import time

 a = np.random.rand(1000000)

b = np.random.rand(1000000)

tic = time.time()       

c = np.dot(a, b)       #向量化运算

toc = time.time()

print("c: %f" % c)

print("vectorized version:" + str(1000*(toc-tic)) + "ms")

#采用for循环语句进行编程

c = 0

tic = time.time()

for i in range(1000000):

    c += a[i] * b[i]

toc = time.time()

print("c: %f" % c)

print("for loop:" + str(1000*(toc-tic)) + "ms")

运行结果：

c: 250099.479223

vectorized version:32.00173377990723ms

c: 250099.479223

for loop:1680.09614944458ms  #可见，向量化的实现代码速度上有飞速提升，而且代码看起#来更加简洁。

##说明，无论有多长的数据列表并且需要对他们进行数学转换，考虑将这些python数据

结构转换为numpy.ndarray对象并使用固有的矢量化功能。

1. python广播

     当两个数组中每个元素都进行相应的运算的时候，需要两个数组的形状相同，如果形状不同，则使python的广播机制进行处理。例如,当一个向量（一维数组）和一个标量（零维数组）相加时，为了能够执行加法，标量需扩展为向量，这种通用机制称为广播。

3.1广播数组：

“广播”的一个工作原则是：两个数组的维度应该相同（即要对一个二维数组进行广播，那么用来广播的数组也应该是二维的），并且只能有一个维度的长度允许不一样，且那个不一样的维度在用来广播的数组里面的长度应该为1（比如，对于一个（3,4）的二维数组，那么用来广播的数组必须是（3,1）或（1,4）；比如对于一个三维的数组（3,4,5），用来广播的数组必须是（1,4,5）或（3,1,5）或（3,4,1）），这样子，我们才说两个数组是广播兼容的。广播会在沿着长度为1的那个维度进行扩散进行。（广播原则：如果两个数组的后缘维度（即：从末尾算起的维度）的轴长相符或者其中的一方长度为1，则认为广播兼容，广播在缺失和长度为1的轴上进行）

如下实例：说明广播是如何操作的：重塑、扩展

import numpy as np

a=np.arange(0,60,10).reshape(-1,1) #建立一个二维数组，形状数（6，1）

print(a.shape)

print(a)

b=np.arange(0,5)   #建立一个一维数组b（向量），形状为（5，）

print(b.shape)

print(b)

c=a+b          #注意：此处向量需要被广播，第一运算步骤为：重塑，将向量的形状从（5，）转换为（1，5）。第二步运算是扩展，将向量的形状从（1，5）转换为（6，5）。              #注意：形状（n,）不能自动广播到向量（m,n）

print(c.shape)

print(c)

运行结果：

(6, 1)

[[ 0]

 [10]

 [20]

 [30]

 [40]

 [50]]

(5,)

[0 1 2 3 4]

(6, 5)

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

解释：

首先b.shape=(1,5)  #由于a与b的维数不一样，首先需让b的维度（shape

#属性性）向a对齐，即向量变为矩阵

print(b.shape)

print(b)

其次，加法的两个输入数组属性分别为（6，1）和（1，5），输出数组的各个轴的长度为输入数组各个轴的长度的最大值，则输出数组的属性为（6，5）；将b在第0轴进行复制，a在第一轴上进行复制。

结果为：

a=a.repeat(5,axis=1)      

print(a)

b=b.repeat(6,axis=0)

print(b)

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

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

注意：numpy内部不会使用repeat进行数据扩展，而是使用内部集成的函数ogrid（创建广播预算用的数组）和mgrid函数（返回是进行广播后的数组）

3.2 python的广播方便与计算：

①　一维向量+常量

import numpy as np

vector=np.arange(4)

b=vector+1.

print(b.shape)

print(b)    #result为：(4,)    向量[1. 2. 3. 4.]

②　多维向量+常数

③　多维向量+行向量

④　多维向量+列向量

a=np.array([[1,2,3],[4,5,6]])

b=[1,2,3]

c=[[4],[5]]

print(a)

print(a+1)

print(a+b)

print(a+c)     #运行结果：

[[1 2 3]

 [4 5 6]]

[[2 3 4]

 [5 6 7]]

[[2 4 6]

 [5 7 9]]

[[ 5  6  7]

 [ 9 10 11]]