Cython：创建数组抛出“常量表达式中不允许”

Question

我尝试将复杂的函数从 Python 重写为 Cython 以大幅加快速度，但遇到以下问题：在使用

编译我的函数 hh_vers_vector.pyx 时

setup(
    ext_modules=cythonize("hh_vers_vector.pyx"),
)

它会抛出以下错误

    cdef int numSamples = len(Iext);

    # initial values
    cdef float v[numSamples]
                      ^
    ------------------------------------------------------------

    hh_vers_vector.pyx:47:27: Not allowed in a constant expression

但是，如果我将“numSamples”作为数字输入函数，则没有问题。我不明白，因为我认为 len(Iext) 也会返回一个数字 10.000。那么为什么 Cython 对这个表达式有问题呢？

cdef float v[numSamples] # ERROR
cdef float v[10000]      # NO ERROR

到目前为止，我的函数看起来像这样：

from math import exp
import time

def hhModel(*params, Iext, float dt, int Vref):

    ## Unwrap params argument: these variables are going to be optimized
    cdef float ENa = params[0]
    cdef float EK  = params[1]
    cdef float EL  = params[2]
    cdef float GNa = params[3]
    cdef float GK  = params[4]
    cdef float GL  = params[5]

    ## Input paramters
    # I  : a list containing external current steps, your stimulus vector [nA/1000]
    # dt : a crazy time parameter [ms]
    # Vref : reference potential [mV]
    # T  : Total simulation time in [ms]

    def alphaM(float v, float vr):       return 0.1 * (v-vr-25) / ( 1 - exp(-(v-vr-25)/10) )
    def betaM(float v, float vr):        return 4 * exp(-(v-vr)/18)
    def alphaH(float v, float vr):       return 0.07 * exp(-(v-vr)/20)
    def betaH(float v, float vr):        return 1 / ( 1 + exp( -(v-vr-30)/10 ) )
    def alphaN(float v, float vr):       return 0.01 * (v-vr-10) / ( 1 - exp(-(v-vr-10)/10) )
    def betaN(float v, float vr):        return 0.125 * exp(-(v-vr)/80)

    ## steady-state values and time constants of m,h,n

    def m_infty(float v, float vr):      return alphaM(v,vr) / ( alphaM(v,vr) + betaM(v,vr) )
    def h_infty(float v, float vr):      return alphaH(v,vr) / ( alphaH(v,vr) + betaH(v,vr) )
    def n_infty(float v, float vr):      return alphaN(v,vr) / ( alphaN(v,vr) + betaN(v,vr) )

    ## parameters
    cdef float Cm, gK, gL, INa, IK, IL, dv_dt, dm_dt, dh_dt, dn_dt, aM, bM, aH, bH, aN, bN
    cdef float Smemb = 4000    # [um^2] surface area of the membrane
    cdef float Cmemb = 1       # [uF/cm^2] membrane capacitance density
    Cm = Cmemb * Smemb * 1e-8  # [uF] membrane capacitance

    gNa = GNa * Smemb * 1e-8   # Na conductance [mS]
    gK  = GK  * Smemb * 1e-8   # K conductance [mS]
    gL  = GL  * Smemb * 1e-8   # leak conductance [mS]

    ######### HERE IS THE PROBLEM ##############
    cdef int numSamples = len(Iext);
    ######### HERE IS THE PROBLEM ##############

    # initial values
    cdef float v[numSamples]
    cdef float m[numSamples]
    cdef float h[numSamples]
    cdef float n[numSamples]

    v[0]  = Vref                    # initial membrane potential
    m[0]  = m_infty(v[0], Vref)     # initial m
    h[0]  = h_infty(v[0], Vref)     # initial h
    n[0]  = n_infty(v[0], Vref)     # initial n

    ## calculate membrane response step-by-step
    for j in range(0, numSamples-1):

        # ionic currents: g[mS] * V[mV] = I[uA]
        INa = gNa * m[j]*m[j]*m[j] * h[j] * (ENa-v[j])
        IK = gK * n[j]*n[j]*n[j]*n[j] * (EK-v[j])
        IL = gL * (EL-v[j])

        # derivatives
        # I[uA] / C[uF] * dt[ms] = dv[mV]
        dv_dt = ( INa + IK + IL + Iext[j]*1e-3) / Cm;

        aM = 0.1 * (v[j]-Vref-25) / ( 1 - exp(-(v[j]-Vref-25)/10))
        bM = 4 * exp(-(v[j]-Vref)/18)
        aH = 0.07 * exp(-(v[j]-Vref)/20)
        bH = 1 / ( 1 + exp( -(v[j]-Vref-30)/10 ) )
        aN = 0.01 * (v[j]-Vref-10) / ( 1 - exp(-(v[j]-Vref-10)/10) )
        bN = 0.125 * exp(-(v[j]-Vref)/80)

        dm_dt = (1-m[j])* aM - m[j]*bM
        dh_dt = (1-h[j])* aH - h[j]*bH
        dn_dt = (1-n[j])* aN - n[j]*bN

        # calculate next step
        v[j+1] = (v[j] + dv_dt * dt)
        m[j+1] = (m[j] + dm_dt * dt)
        h[j+1] = (h[j] + dh_dt * dt)
        n[j+1] = (n[j] + dn_dt * dt)

    return v

Answer 1

顾名思义，Cython 中的

cdef 用于 C 定义。因此代码

cdef float v[10000]

被翻译成下面的C代码

float v[10000];

意思是 static 大小为 10k 的浮点数组，在编译时定义。

另一方面

cdef float v[numSamples]

将转换为 C 代码

int numSamples = <..>;
float v[numSamples];

这是试图编译包含可变数量元素的静态数组，这不是有效的 C。因此会出现“常量表达式错误”。

要么使用 python 列表来存储浮点数，要么使用dynamically allocate C arrays via malloc/free。