Cython: créer un tableau throws "non autorisé dans une expression constante"

Question

J'essaie de réécrire une fonction complexe de Python à Cython pour l'accélérer massivement et je rencontre le problème suivant: lors de la compilation de ma fonction hh_vers_vector.pyx en utilisant

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

il jette l'erreur suivante

cdef int numSamples = len(Iext); 

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

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

cependant, si je donne « numSamples » comme numéro dans la fonction, il n'y a pas de problème. Je ne comprends pas cela, parce que je pensais que len (Iext) retournera un nombre de 10.000 aussi bien. Alors pourquoi Cython a-t-il un problème avec cette expression?

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

Ma fonction ressemble à ceci jusqu'à présent:

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

cdef dans Cython est pour les définitions de C, comme son nom l'indique. Par conséquent, le code

cdef float v[10000] 

se traduit au code C suivant

float v[10000]; 

Signification un tableau de flottaison statique de taille 10k, définie à temps compilation.

Ceci, d'autre part

cdef float v[numSamples] 

se traduirait par le code C

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

Ceci est une tentative de compiler un tableau statique contenant un nombre variable d'éléments, ce qui est pas valide C. D'où l'erreur d'expression constante. Utilisez soit une liste Python pour stocker les flottants, soit dynamically allocate C arrays via malloc/free.