NP1 Introduction
1. 示例:从 python 到 numpy
numpy的核心是向量化,这和我们通常采取的面向过程编程的思维不太一样。我们将会大量的面对 "vectors", "arrays", "views" or "ufuncs"。
引用
举个简单的示例来说明我们的思维应该采取怎样的转变——随机漫步。一种可能的面向对象做法是定义一个 RandomWalker 类,然后写一个 walk 方法,每走一步(随机的一步)就返回当前的位置。
1.1 面向对象方法
yield 是 python 中用于定义生成器的关键字。
class RandomWalker:
def __init__(self):
self.position = 0
def walk(self, n):
self.position = 0
for i in range(n):
yield self.position
self.position += 2*random.randint(0, 1) - 1
walker = RandomWalker()
walk = [position for position in walker.walk(1000)]最终 walk 数组保存了随机游走了 0 ~ 999 次后的位置。我们测试其性能:
>>> from tools import timeit
>>> walker = RandomWalker()
>>> timeit("[position for position in walker.walk(n=10000)]", globals())
10 loops, best of 3: 15.7 msec per loop1.2 程序化方法
对于这么简单的问题,我们大概可以省去类的定义,只专注于 walk 方法,也就是计算每次随机移动后的连续位置。
def random_walk(n):
position = 0
walk = [position]
for i in range(n):
position += 2*random.randint(0, 1)-1
walk.append(position)
return walk
walk = random_walk(1000)这个新方法省了一些 CPU 周期,但也没省多少,因为这个函数跟面向对象那套写法基本上一模一样,省下来的那点周期多半还是来自 Python 内部面向对象机制的开销。
1.3 向量化方法
我们用 Python 的 itertools 模块可以做得更好,它提供了一组函数,能创建迭代器来实现高效循环。如果我们注意到随机游走其实就是一步步累积起来的结果,就可以改写这个函数,通过累积和从而避免循环:
def random_walk_faster(n=1000):
from itertools import accumulate
# Only available from Python 3.6
steps = random.choices([-1,+1], k=n)
return [0]+list(accumulate(steps))
walk = random_walk_faster(1000)上述方法实则将函数向量化了:我们没有用循环来一步步挑出连续的步长再加到当前位置上,而是一次性生成所有步长,然后用 accumulate 函数来算出所有位置。这样换取到了速度的提升:
>>> from tools import timeit
>>> timeit("random_walk_faster(n=10000)", globals())
10 loops, best of 3: 2.21 msec per loop我们很容易将上述写法迁移到 numpy 写法中,而且速度提高了 500 倍:
def random_walk_fastest(n=1000):
# No 's' in numpy choice (Python offers choice & choices)
steps = np.random.choice([-1,+1], n)
return np.cumsum(steps)
walk = random_walk_fastest(1000)
>>> from tools import timeit
>>> timeit("random_walk_fastest(n=10000)", globals())
1000 loops, best of 3: 14 usec per loop2. 可读性 v.s. 速度
def function_1(seq, sub):
return [i for i in range(len(seq) - len(sub)) if seq[i:i+len(sub)] == sub]
def function_2(seq, sub):
target = np.dot(sub, sub)
candidates = np.where(np.correlate(seq, sub, mode='valid') == target)[0]
check = candidates[:, np.newaxis] + np.arange(len(sub))
mask = np.all((np.take(seq, check) == sub), axis=-1)
return candidates[mask]上述两个函数都尝试在序列 seq 中寻找子序列 sub 出现的所有起始位置。但是第二个numpy写的明显没啥可读性,但就很快。
参考
From Python to NumPy
cicada@blog:~