Cython没有明显的速度:如何处理Cython中的Numpy广播



我是Cython的初学者。我正在尝试加快称为多次的功能,并且大幅提高将非常有帮助。

我的原始版本此功能大量使用了大型多维阵列上的numpy广播。当尝试使该功能含糊不清时,我最初尝试将这些数组转换为Cython MemoryViews。但是,我遇到了编译错误,因为该函数使用这些数组执行算术,此后我了解到这不受MemoryViews的支持。

因此,我修改了我的冷冻代码,将数组声明为numpy ndarrays,而不是内存视图。现在,Cythonized功能有效,但与原始的Pure Python/numpy版本相比,速度没有任何明显的差异。

因此,我现在可以:克服需要执行的算术。

要继续选项b。(,我希望在此功能的for循环中解决算术。我可能会编写其他Cython函数,以逐元素的基础执行数组乘法和添加,我知道这将允许我使用MemoryViews。但是,鉴于此代码的复杂数量广播,我希望这可能涉及大量的努力(而且我不一定知道如何开始...(。此外,考虑到按元素进行元素(肯定在C中公认的(,我不确定这种努力是否会取得成果。

实际上可能不会更快。

我非常欢迎任何建议或支持。谢谢。

from numpy import pi, exp, sin, cos
import numpy as np
cimport numpy as np
cimport cython
cdef np.ndarray[np.double_t, ndim=3] foo(bar, double dt, f, xi):
    cdef int nt, nh, nj
    cdef Py_ssize_t t
    nt = bar.shape[0]
    nh = bar.shape[1]
    if len(bar.shape) < 3:
        bar = bar[:, np.newaxis, :]
    cdef np.ndarray[np.double_t, ndim=3] bar_c = bar
    nj = len(f)
    k = (2 * pi * f) ** 2
    wn = k ** 0.5
    wd = (wn * (1 - xi ** 2) ** 0.5)
    cdef np.ndarray[np.double_t, ndim=3] u = np.zeros((nt, nj, nh))
    cdef np.ndarray[np.double_t, ndim=3] v = np.zeros((nt, nj, nh))
    C1 = exp(-xi * wn * dt)
    C2 = sin(wd * dt)
    C3 = cos(wd * dt)
    cdef np.ndarray[np.double_t, ndim=2] A11 = C1 * (C3 + (xi / (1 - xi ** 2) ** 0.5) * C2)
    cdef np.ndarray[np.double_t, ndim=2] A12 = (C1 / wd) * C2
    cdef np.ndarray[np.double_t, ndim=2] A21 = (-wn / (1 - xi ** 2) ** 0.5) * C1 * C2
    cdef np.ndarray[np.double_t, ndim=2] A22 = C1 * (C3 - (xi / (1 - xi ** 2) ** 0.5) * C2)
    cdef np.ndarray[np.double_t, ndim=2] B11 = C1 * (
        (((2 * xi ** 2 - 1) / (dt * wn ** 2)) + xi / wn) * C2 / wd
        + ((2 * xi / (dt * wn ** 3)) + (1 / wn ** 2)) * C3
    ) - 2 * xi / (dt * wn ** 3)
    cdef np.ndarray[np.double_t, ndim=2] B12 = (
        -C1
        * (
            ((2 * xi ** 2 - 1) / (dt * wn ** 2)) * C2 / wd
            + ((2 * xi) / (dt * wn ** 3)) * C3
        )
        - (1 / wn ** 2)
        + 2 * xi / (dt * wn ** 3)
    )
    cdef np.ndarray[np.double_t, ndim=2] B21 = -A12 - ((A11 - 1) / (dt * wn ** 2))
    cdef np.ndarray[np.double_t, ndim=2] B22 = -B21 - A12
    for t in range(0, nt - 1):
        u[t + 1, :, :] = (
            A11 * u[t, :, :]
            + A12 * v[t, :, :]
            + B11 * bar_c[t, :, :]
            + B12 * bar_c[t + 1, :, :]
        )
        v[t + 1, :, :] = (
            A21 * u[t, :, :]
            + A22 * v[t, :, :]
            + B21 * bar_c[t, :, :]
            + B22 * bar_c[t + 1, :, :]
        )
    cdef np.ndarray[np.double_t, ndim=3] out = -2 * xi * wn * v - (wn ** 2) * u - bar_c
    return out

道歉不减小代码的大小。鉴于查询的性质,我很难确定一个最小的,可再现的例子。

谢谢@codesurgeon和 @9000的答案。您已经确认/建议,在某些情况下,C代码确实可以使Numpy操作的速度提高 - 即使考虑到Cython MemoryViews需要逐件执行C操作,而我的原始代码是使用在大型阵列上广播的Numpy操作。

这促使我进一步探索了事情。在元素基础上,在内存视图上操作,该代码现在更快地介于3〜40倍之间(非常依赖于输入数组的大小,这将会有所不同(。

修改基本上看起来像这样:

    # Memoryview declarations
    cdef double[:] u_mv = u
    cdef double[:] v_mv = v
    cdef double[:, :, :] out_mv = out
    cdef double[:, :, :] bar_mv = bar
    for j in range(nj):
        ...same calculation of constants C1, C2 etc. as before...
        for h in range(nh):
            for t in range(0, nt - 1):
                u_mv[t + 1] = (
                    A11 * u_mv[t]
                    + A12 * v_mv[t]
                    + B11 * bar_mv[t, 0, h]
                    + B12 * bar_mv[t + 1, 0, h]
                )
                v_mv[t + 1] = (
                    A21 * u_mv[t]
                    + A22 * v_mv[t]
                    + B21 * bar_mv[t, 0, h]
                    + B22 * bar_mv[t + 1, 0, h]
                )
                out_mv[t + 1, j, h] = (
                    -2 * xi_j * wn * v_mv[t + 1] - (wn ** 2) * u_mv[t + 1]
                    - bar_mv[t + 1, 0, h]
                )
            out_mv[0, j, h] = -bar_mv[0, 0, h]
    return out

再次感谢。

最新更新