NP1 Introduction

· Tech

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 loop

1.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 loop

2. 可读性 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写的明显没啥可读性,但就很快。

参考

  1. From Python to NumPy

    https://www.labri.fr/perso/nrougier/from-python-to-numpy/

cicada@blog:~