Cython 1D 归一化滑动互相关优化

Question

我有以下代码，它执行标准化互相关，在 python 中寻找两个信号的相似性：

def normcorr(template,srchspace):
template=(template-np.mean(template))/(np.std(template)*len(template)) # Normalize template
CCnorm=srchspace.copy()
CCnorm=CCnorm[np.shape(template)[0]:] # trim CC matrix
for a in range(len(CCnorm)):
    s=srchspace[a:a+np.shape(template)[0]]
    sp=(s-np.mean(s))/np.std(s)
    CCnorm[a]=numpy.sum(numpy.multiply(template,sp))
return CCnorm

但是你可以想象它太慢了。查看 cython 文档，在原始 python 中执行循环时，速度会大大提高。因此，我尝试编写一些 cython 代码，其中包含如下所示的变量数据类型：

from __future__ import division
import numpy as np
import math as m
cimport numpy as np
cimport cython
def normcorr(np.ndarray[np.float32_t, ndim=1] template,np.ndarray[np.float32_t, ndim=1]  srchspace):
    cdef int a
    cdef np.ndarray[np.float32_t, ndim=1] s
    cdef np.ndarray[np.float32_t, ndim=1] sp
    cdef np.ndarray[np.float32_t, ndim=1] CCnorm
    template=(template-np.mean(template))/(np.std(template)*len(template))
    CCnorm=srchspace.copy()
    CCnorm=CCnorm[len(template):]
    for a in range(len(CCnorm)):
        s=srchspace[a:a+len(template)]
        sp=(s-np.mean(s))/np.std(s)
        CCnorm[a]=np.sum(np.multiply(template,sp))
    return CCnorm

但是一旦我编译它，代码实际上比纯 python 代码运行得慢。我在这里 (How to call numpy/scipy C functions from Cython directly, without Python call overhead?) 发现从 cython 调用 numpy 可能会显着减慢代码速度，这是我的代码的问题吗，在这种情况下，我必须定义内联函数来替换对 np 的所有调用，或者还有什么我我做错了吗？

Answer 1

因为你在 cython 循环中调用 numpy 函数，所以不会有速度提升。

如果你使用pandas，你可以在numpy中使用roll_mean()和roll_std()和convolve()进行非常快的计算，代码如下：

import numpy as np
import pandas as pd

np.random.seed()

def normcorr(template,srchspace):
    template=(template-np.mean(template))/(np.std(template)*len(template)) # Normalize template
    CCnorm=srchspace.copy()
    CCnorm=CCnorm[np.shape(template)[0]:] # trim CC matrix
    for a in range(len(CCnorm)):
        s=srchspace[a:a+np.shape(template)[0]]
        sp=(s-np.mean(s))/np.std(s)
        CCnorm[a]=np.sum(np.multiply(template,sp))
    return CCnorm

def fast_normcorr(t, s):
    n = len(t)
    nt = (t-np.mean(t))/(np.std(t)*n)
    sum_nt = nt.sum()
    a = pd.rolling_mean(s, n)[n-1:-1]
    b = pd.rolling_std(s, n)[n-1:-1]
    b *= np.sqrt((n-1.0) / n)
    c = np.convolve(nt[::-1], s, mode="valid")[:-1]
    result = (c - sum_nt * a) / b    
    return result

n = 100
m = 1000
t = np.random.rand(n)
s = np.random.rand(m)

r1 = normcorr(t, s)
r2 = fast_normcorr(t, s)
assert np.allclose(r1, r2)

您可以检查结果r1 和r2 是否相同。这是timeit 测试：

%timeit normcorr(t, s)
%timeit fast_normcorr(t, s)

输出：

10 loops, best of 3: 59 ms per loop
1000 loops, best of 3: 273 µs per loop

速度提高了 200 倍。

Answer 2

如果您使用cython -a 编译代码并查看 HTML 输出，您会发现 Python 开销很大。

@cython.boundscheck(False)
@cython.cdivision(True)   # Don't check for divisions by 0
def normcorr(np.ndarray[np.float32_t, ndim=1] template,np.ndarray[np.float32_t, ndim=1]  srchspace):
    cdef int a
    cdef int N = template.shape[0]
    cdef NCC = srchspace.shape[0] - N
    cdef np.ndarray[np.float32_t, ndim=1] s
    cdef np.ndarray[np.float32_t, ndim=1] sp
    cdef np.ndarray[np.float32_t, ndim=1] CCnorm
    template=(template - template.mean()) / (template.std() * N)
    CCnorm=srchspace[N:].copy()    # You don't need to copy the whole array
    for a in xrange(NCC):  # Use xrange in Python2
        s=srchspace[a:a+N]
        sp=(s-np.mean(s)) / np.std(s)
        CCnorm[a]= (template * sp).sum()
    return CCnorm

为了提高性能，可以优化最后两行：

@cython.boundscheck(False)
@cython.cdivision(True)
cdef multiply_by_normalised(np.ndarray[np.float32_t, ndim=1] template, np.ndarray[np.float32_t, ndim=1] s):
    cdef int i
    cdef int N = template.shape[0]
    cdef float_32_t mean, std, out = 0

    mean = s.mean()
    std = s.std()

    for i in xrange(N):
        out += (s[i] - mean) / std * template[i]
    return out

如果还需要挤更多的时间，可以使用bottleneck的mean和std函数，比Numpy快。