Songyuli's blog
0%

向量化与面向对象

Posted on Disqus:

在这一篇文章中,我们将讨论两个话题:向量化和面向对象

第一个话题大家应该已经有相关的知识了,但是”可能”还缺乏一个非常清晰的理解;
而第二个话题大家可能听说的不多,虽然像”面向对象”这样的话题对于计算机程序设计非常重要,但是这对于我们做科学计算实际上意义有限,但是基本概念又不能不知道

向量化

这一部分主要是关于numpy中ndarray类的

虽然我的意思是,关于科学计算,真正需要”掌握”的就是ndarray,其他查一下就能现学现用

原理

在numpy中存在这样的操作,将一个数组与另一个数组(或者标量)进行数学运算时,这些数据会被”广播”进行计算

1
2
3
4
import numpy as np
a=np.array([1,2,3,4,5])
b=1
a+b
array([2, 3, 4, 5, 6])

这些计算更多的可以被看做是对于线性代数的一种扩展,多了一种标量加矢量的运算

效果

在numpy中是使用C编写的代码,使用了大量的底层运算,其速度远快于构造python中的循环
下面就是一个例子来比较两种代码的效果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def use_list():
    x=[1]*1000 
    for i in range(len(x)):
        x[i]+=1

def use_ndarray():
    x=np.ones(1000)
    x+=1

#计时
%timeit use_list()
%timeit use_ndarray()


67.7 µs ± 9.24 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
3.97 µs ± 193 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

可以非常容易地看出,使用python列表(前者)的效率会远低于使用ndarray(后者)

shape

前文提到,向量化的操作与线性代数有千丝万缕的联系.因此,数组的形状就显得非常重要
可以直接访问shape属性获得关于数组的形状信息

1
2
3
4
a=np.array([[1,2,3],[4,5,6]])
b=np.array([7,8])
print("The shape of a is {}".format(a.shape))
print("The shape of b is {}".format(b.shape))
The shape of a is (2, 3)
The shape of b is (2,)

而对于二维这种最为常见的情况来说,a[i,j]是先行后列.

而为了实现计算,需要修改数组的形状,例如将b变为一个行向量可以进行矩阵乘法

1
2
c=b.reshape(1,2)
np.dot(c,a)
array([[39, 54, 69]])

这种reshape在处理解决实际问题的时候往往非常有用,在我个人的理解中,如果需要批量的计算,往往会把数组变到二维便于进行计算(更高的维度也不是不可以,但是会比较复杂).然后等到计算完成之后再变到自己所需的shape
同时,reshape也支持使用-1缺省

这里有一个需要注意的地方,reshape返回的是数组的浅拷贝,因此如果对reshape之后的元素进行修改,那么会影响到原来的数组.

1
2
c[0,1]=1
b
array([7, 1])

axis

如果执行一些操作(例如说sum),其就会向你提供一个可选参数axis,是指这个运算沿着哪个轴

下面就展示一下沿着不同的轴进行操作的效果

1
2
print("Along axis 0 is {}".format(a.sum(axis=0)))
print("Along axis 1 is {}".format(a.sum(axis=1)))
Along axis 0 is [5 7 9]
Along axis 1 is [ 6 15]

可以将axis的效果简单地理解为如果沿着某个轴进行一些聚合操作(例如求和),那么执行操作之后数组的shape的这个元素就会变为1

布尔数组与布尔索引

在向量化中,还有一个重要的工具,记为布尔数组
我们可以对一个数组进行比较运算,得到布尔数组

1
a>4
array([[False, False, False],
       [False,  True,  True]])

而这样的布尔数组一方面可以用作输入实现分类,另外一方面也可以作为索引来提取满足条件的元素
例如下面的代码展示如何将a中大于4的元素提取出来,以及按照(-1,1)分类

1
2
print(a[a>4])
print(np.where(a>4,1,-1))
[5 6]
[[-1 -1 -1]
 [-1  1  1]]

使用这样的功能,我们可以实现一些原本看起来难以实现向量化的运算
例如下面一段代码模拟了抛针法估算$\pi$的过程

1
2
3
4
5
6
7
8
9
10
11
def calc_pi(sample_number:int,ran_num:int):
    ran=np.random.RandomState(ran_num)
    x=ran.uniform(-0.5,0.5,sample_number)#生成中心坐标
    theta=ran.uniform(0,2*np.pi,sample_number)#生成角度
    r1=x-0.5*np.cos(theta)
    r2=x+0.5*np.cos(theta)
    k=r1*r2
    num=np.sum(np.where(k<=0,1,0))
    return 2*sample_number/num

calc_pi(1000000,1)
3.139934312574181

实际上,关于向量化,大家还有很多需要学,我这里展示的仅为冰山一角

“流畅地编程”:面向对象和一些特性

我一开始并不打算讲面向对象,一方面在于面向对象本质太过复杂,完全可以涉及到另外一个完整而复杂的学科,而另外一方面则更为重要,在于面向对象对于我们这种”非专业”的程序设计者尤其如果你仅仅希望将python看成一个分析数据的脚本的人,面向对象是一个完全不需要涉及的内容;并且如果滥用里面的一些东西,无疑会造成很大灾难(不信的可以看看我Github上的几个屎山仓库)

不过稍微学习一下相关内容也会有用,至少你不使用,但是在别人使用的时候你要看懂发生了什么

另外,我将介绍一些python中的特性,这些特性和前面说的向量化一样,往往会被人忽视(因为没有它也能工作)

但是请不要过于依赖一些偏门特性,正如python作者Tim Peters所说

“Here’s the plan: When someone uses a feature you don’t understand, simply shoot them. This is easier than learning something new, and before too long the only living coders will be writing in an easily understood, tiny subset of Python 0.9.6 .”

面向对象:它是什么,会给我们带来什么影响

Object-oriented programming (OOP) is a computer programming model that organizes software design around data, or objects, rather than functions and logic. An object can be defined as a data field that has unique attributes and behavior.

上面这段话摘自这个网页,其说明了OOP的一个基本特点:如何处理事物是次要的,主要是将事物抽象成一个个对象
这样的抽象有好处也有坏处
首先讲好处,其可以将更多东西抽象成一个,不仅仅是函数那样输入输出,这可以极大提升代码规整度并且在实现或改进新功能时尽可能不影响到用户(这里的用户是更加高层的调用)的调用,将各个地方实现”脱耦”
但是相应的存在缺陷,一个大问题就在于面向对象编程基本的可读性会变差(对于小程序),同时可能导致开销增大
实际上python自己就是非常依赖对象进行工作的,关于什么是属性,什么是方法,大家可以去参考这里
至于为何我不讲,一方面是因为我懒,另外一方面这个东西我将可能会误人子弟.大家看到类的基本特征就行了,像继承多态,装饰器等东西大家想看就看,不想看也不勉强.

python中的一些技巧

这里面收录了一些我看到的优雅的代码,就将其收集下来,以供参考

使用迭代器生成形状

在我们进行科学绘图的时候,总是会遇见一些让你眼前一黑的情况,比如要对于未知个特征的数据进行绘图,不同绘图的线形又要单独指定.这样的场景并不少见,但每一次碰到都会让人非常痛苦,下面是我看到的一个有意思的做法,可以用来参考

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import matplotlib.pyplot as plt
def line_gen():  
    """line_gen 产生不超过九个线形

    Yields
    ------
    str
        产生得到线形的字符串
    """    
    line_shapes=["-","--","-."]
    line_cols=["r","g","b"]
    for line_shape in line_shapes:
        for line_col in line_cols:
            yield line_col+line_shape
            
indata=np.linspace(1,5,100)
pdatas=[indata**i for i in range(5)]
ls=line_gen()#创建一个生成器实例
for pdata in pdatas:
    plt.plot(indata,pdata,next(ls))
plt.show()

​
png
​

可以看到,迭代器完美地处理了这样的情况.
其原理可以大概这样解释:

  1. 解释器注意到函数中的yield关键字,因此在调用line_gen()时不会运行这个函数,而是返回一个生成器对象(generator object)
  2. 之后每次调用next中以这个生成器对象作为参数,那么函数就会运行到下一个yield处并返回(需要注意到的是,再一次用同一生成器调用next的时候,函数保留了原来的”记忆”)
  3. 如果在next时函数以及执行完了,那么就会发生StopIteration错误

    yield的作用与其说”返回”,不如说是”产生”,因为函数并没有终止