用於改進Cython代碼的高效矩陣向量結構

Question

我必須在SODEs系統上運行一些模擬。因爲我需要使用隨機圖，所以我認爲使用python生成圖的相鄰矩陣然後用C生成模擬是一個好主意。所以我轉向了cython。

爲了儘可能提高速度，我編寫了我的代碼，提示cython documentation。但知道我真的不知道我的代碼是否好。我也運行cython toast.pyx -a，但我不明白這些問題。

• 在cython中用於向量和矩陣的最佳結構是什麼？我應該如何在我的代碼上使用np.array或double來定義例如bruit？請注意，我將比較矩陣（0或1）的元素以便進行總和或不是。結果將是一個矩陣NxT，其中N是系統的維數，T是我想用於模擬的時間。
• 我在哪裏可以找到double[:]的文檔？
• 如何在函數的輸入中聲明向量和矩陣（下面的G，W和X）？我怎麼能聲明一個向量double？

不過我讓我的代碼聊我：

from __future__ import division 
import scipy.stats as stat 
import numpy as np 
import networkx as net 

#C part 
from libc.math cimport sin 
from libc.math cimport sqrt 

#cimport cython 
cimport numpy as np 
cimport cython 
cdef double tau = 2*np.pi #http://tauday.com/ 

#graph 
def graph(int N, double p): 
    """ 
    It generates an adjacency matrix for a Erdos-Renyi graph G{N,p} (by default not directed). 
    Note that this is an O(n^2) algorithm and it gives an array, not a (sparse) matrix. 

    Remark: fast_gnp_random_graph(n, p, seed=None, directed=False) is O(n+m), where m is the expected number of edges m=p*n*(n-1)/2. 
    Arguments: 
     N : number of edges 
     p : probability for edge creation 
    """ 
    G=net.gnp_random_graph(N, p, seed=None, directed=False) 
    G=net.adjacency_matrix(G, nodelist=None, weight='weight') 
    G=G.toarray() 
    return G 


@cython.boundscheck(False) 
@cython.wraparound(False) 
#simulations 
def simul(int N, double H, G, np.ndarray W, np.ndarray X, double d, double eps, double T, double dt, int kt_max): 
    """ 
    For details view the general description of the package. 
    Argumets: 
     N : population size 
     H : coupling strenght complete case 
     G : adjiacenty matrix 
     W : disorder 
     X : initial condition 
     d : diffusion term 
     eps : 0 for the reversibily case, 1 for the non-rev case 
     T : time of the simulation 
     dt : increment time steps 
     kt_max = (int) T/dt 
    """ 
    cdef int kt 
    #kt_max = T/dt to check 
    cdef np.ndarray xt = np.zeros([N,kt_max], dtype=np.float64) 
    cdef double S1 = 0.0 
    cdef double Stilde1 = 0.0 
    cdef double xtmp, xtilde, x_diff, xi 

    cdef np.ndarray bruit = d*sqrt(dt)*stat.norm.rvs(N) 
    cdef int i, j, k 

    for i in range(N): #setting initial conditions 
     xt[i][0]=X[i] 

    for kt in range(kt_max-1): 
     for i in range(N): 
      S1 = 0.0 
      Stilde1= 0.0 
      xi = xt[i][kt] 

      for j in range(N): #computation of the sum with the adjiacenty matrix 
       if G[i][j]==1: 
        x_diff = xt[j][kt] - xi 
        S2 = S2 + sin(x_diff) 

      xtilde = xi + (eps*(W[i]) + (H/N)*S1)*dt + bruit[i] 

      for j in range(N): 
       if G[i][j]==1: 
        x_diff = xt[j][kt] - xtilde 
        Stilde2 = Stilde2 + sin(x_diff) 

      #computation of xt[i] 
      xtmp = xi + (eps*(W[i]) + (H/N)*(S1+Stilde1)*0.5)*dt + bruit 
      xt[i][kt+1] = xtmp%tau 

    return xt 

非常感謝您！

更新

我改變了變量定義的順序，對np.array和double向xt[i][j]和xt[i,j]的矩陣long long。代碼現在非常快，HTML文件中的黃色部分就在聲明的周圍。謝謝！

from __future__ import division 
import scipy.stats as stat 
import numpy as np 
import networkx as net 

#C part 
from libc.math cimport sin 
from libc.math cimport sqrt 

#cimport cython 
cimport numpy as np 
cimport cython 
cdef double tau = 2*np.pi #http://tauday.com/ 

#graph 
def graph(int N, double p): 
    """ 
    It generates an adjacency matrix for a Erdos-Renyi graph G{N,p} (by default not directed). 
    Note that this is an O(n^2) algorithm and it gives an array, not a (sparse) matrix. 

    Remark: fast_gnp_random_graph(n, p, seed=None, directed=False) is O(n+m), where m is the expected number of edges m=p*n*(n-1)/2. 
    Arguments: 
     N : number of edges 
     p : probability for edge creation 
    """ 
    G=net.gnp_random_graph(N, p, seed=None, directed=False) 
    G=net.adjacency_matrix(G, nodelist=None, weight='weight') 
    G=G.toarray() 
    return G 


@cython.boundscheck(False) 
@cython.wraparound(False) 
#simulations 
def simul(int N, double H, long long [:, ::1] G, double[:] W, double[:] X, double d, double eps, double T, double dt, int kt_max): 
    """ 
    For details view the general description of the package. 
    Argumets: 
     N : population size 
     H : coupling strenght complete case 
     G : adjiacenty matrix 
     W : disorder 
     X : initial condition 
     d : diffusion term 
     eps : 0 for the reversibily case, 1 for the non-rev case 
     T : time of the simulation 
     dt : increment time steps 
     kt_max = (int) T/dt 
    """ 
    cdef int kt 
    #kt_max = T/dt to check 
    cdef double S1 = 0.0 
    cdef double Stilde1 = 0.0 
    cdef double xtmp, xtilde, x_diff 

    cdef double[:] bruit = d*sqrt(dt)*np.random.standard_normal(N) 

    cdef double[:, ::1] xt = np.zeros((N, kt_max), dtype=np.float64) 
    cdef double[:, ::1] yt = np.zeros((N, kt_max), dtype=np.float64) 
    cdef int i, j, k 

    for i in range(N): #setting initial conditions 
     xt[i,0]=X[i] 

    for kt in range(kt_max-1): 
     for i in range(N): 
      S1 = 0.0 
      Stilde1= 0.0 

      for j in range(N): #computation of the sum with the adjiacenty matrix 
       if G[i,j]==1: 
        x_diff = xt[j,kt] - xt[i,kt] 
        S1 = S1 + sin(x_diff) 

      xtilde = xt[i,kt] + (eps*(W[i]) + (H/N)*S1)*dt + bruit[i] 

      for j in range(N): 
       if G[i,j]==1: 
        x_diff = xt[j,kt] - xtilde 
        Stilde1 = Stilde1 + sin(x_diff) 

      #computation of xt[i] 
      xtmp = xt[i,kt] + (eps*(W[i]) + (H/N)*(S1+Stilde1)*0.5)*dt + bruit[i] 
      xt[i,kt+1] = xtmp%tau 

    return xt 

Answer 1

cython -a color codes the cython source。如果你點擊一行，它會顯示相應的C源代碼。作爲一個經驗法則，你不需要任何內部循環中的黃色。

幾件事情在你的代碼跳出：

• x[j][i]創建x[j]在每次調用一個臨時數組，所以使用x[j, i]代替。
• 而不是cdef ndarray x更好提供維度和dtype（cdef ndarray[ndim=2, dtype=float]）或---最好---使用類型化的memoryview語法：cdef double[:, :] x。

例如，而不是

cdef np.ndarray xt = np.zeros([N,kt_max], dtype=np.float64) 

更好地利用

cdef double[:, ::1] xt = np.zeros((N, kt_max), dtype=np.float64) 

• 確保您在緩存友好型內存訪問。例如，確保你的數組是按C順序排列的（上一個維度變化最快），將內存視圖聲明爲double[:, ::1]，並迭代數組，最後一個索引變化最快。

編輯：見http://cython.readthedocs.io/en/latest/src/userguide/memoryviews.html 爲類型memoryview語法double[:, ::1]等