小编典典

numpy复发

python

有没有办法在numpy中不使用for来进行复发?

使用np.addwithout关键字可以解决问题dtype="int"

import numpy as np
N = 100

fib = np.zeros(N, dtype=np.int)
fib[:2] = 1.
np.add(fib[:-2], fib[1:-1], out=fib[2:])

print(fib[:10])

输出: [ 1 1 2 3 5 8 13 21 34 55]

但是,如果dtype更改为np.float

import numpy as np
N = 100

fib = np.zeros(N, dtype=np.float)
fib[:2] = 1.
np.add(fib[:-2], fib[1:-1], out=fib[2:])

print(fib[:10])

输出: [ 1. 1. 2. 1. 3. 1. 4. 1. 5. 1.]

有人可以告诉我为什么吗?或任何其他方式进行递归?


阅读 219

收藏
2021-01-20

共1个答案

小编典典

这是使用scipy过滤器的一种相当有效的方法:

import numpy as np
from scipy import signal

def scipy_fib(n):
    x = np.zeros(n, dtype=np.uint8)
    x[0] = 1

    sos = signal.tf2sos([1], [1, -1, -1])
    return signal.sosfilt(sos, x)

(请注意,lfilter由于信号处理这是一个不稳定的过滤器,导致Python解释器崩溃,因此我们不能使用它,因此我们必须转换为二阶节形式。)

过滤方法的问题在于我不确定是否可以将其推广到解决ODE的实际问题。

另一个解决方案是简单地用Python编写循环并使用JIT编译器对其进行加速:

import numba as nb

@nb.jit(nopython=True)
def numba_fib(n):
    y = np.empty(n)
    y[:2] = 1

    for i in range(2, n):
        y[i] = y[i-1] + y[i-2]

    return y

JIT方法的优点是您可以实现各种花哨的东西,但是如果坚持使用简单的循环和分支而不调用任何(非JITted)Python函数,则效果最佳。

计时结果:

# Baseline without JIT:
%timeit numba_fib(10000)  # 100 loops, best of 3: 5.47 ms per loop

%timeit scipy_fib(10000)  # 1000 loops, best of 3: 719 µs per loop
%timeit numba_fib(10000)  # 10000 loops, best of 3: 33.8 µs per loop

# For fun, this is how Paul Panzer's solve_banded approach compares:
%timeit banded_fib(10000)  # 1000 loops, best of 3: 1.33 ms per loop

scipy过滤器方法比纯Python循环快,但比JITted循环慢。我猜想filter函数是相对通用的,它可以完成我们不需要的东西,而JITted循环则可以编译成很小的循环。

2021-01-20