2017-10-20 93 views
0

我嘗試重寫一個複雜的函數形式Python到Cython以加速它的大規模,我遇到以下問題:編譯我的函數hh_vers_vector.pyx使用Cython:創建一個數組拋出「不允許在一個常量表達式」

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 
+0

[Cython「不允許在常量表達式中」,boundscheck False不起作用]的可能重複(https://stackoverflow.com/questions/46306737/cython-not-allowed-in-a-constant-expression -boundscheck-false-doesnt-work) – DavidW

+0

我看到那篇文章,但它並沒有真正幫助.. –

+0

我認爲@ danny對這個問題的回答有點更好解釋,但它的確說大體上是相同的東西(所以我我不知道爲什麼這個答案有幫助,另一個職位沒有......)。使用另一個問題中給出的內存視圖的建議是最明智的解決方案(比malloc/free更容易) – DavidW

回答

3

cdef在用Cython是對C的定義,顧名思義。因此,代碼

cdef float v[10000] 

被轉換爲下面的C代碼

float v[10000]; 

含義大小10K的靜態 float數組,在編譯時定義。

此,在另一方面

cdef float v[numSamples] 

將轉化成C代碼

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

這是編譯包含可變數目的元件的靜態數組,這是無效的嘗試C.因此,「常量表達式錯誤」。

要麼使用python列表來存儲浮點數,要麼使用dynamically allocate C arrays via malloc/free

+0

感謝您的評論!我猜pythonl列表會比C數組慢嗎?我以前從來沒有用過malloc也沒有free,你知道代碼會如何看起來像嗎? –

+0

順便說一句,使用'DEF numSamples = 200000'也可以,但是我不明白你的解釋了。 –

+1

靜態數組在編譯時不能分配一個可變數字。它們是靜態的。像上面這樣靜態定義的數字是靜態/常量。在運行時更改的int變量不是。 – danny

相關問題