効率的な行列ベクトル構造のcythonコードを改善する

Question

私はSODEのシステム上でいくつかのシミュレーションを実行する必要があります。ランダムグラフを使用する必要があるので、グラフの隣接行列を生成するためにPythonを使用し、次にシミュレーションのためにCを使用することをお勧めします。だから私はcythonに目を向ける。

私はできるだけ速くするために、cython documentationのヒントに沿って自分のコードを書いています。しかし、私のコードが良いかどうかは本当に分かりません。私もcython toast.pyx -aを実行するが、私は問題を理解していない。

• cythonのベクトルと行列にはどのような構造が最適ですか？たとえば、をnp.arrayまたはdoubleというコードでどのように定義する必要がありますか？合計を行うかどうかを判断するために、行列の要素（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 

はどうもありがとうございました！

更新

私はxt[i,j]にdoubleとxt[i][j]にnp.array、変数の定義の順序を変更し、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カラーコードcythonソース。行をクリックすると、対応するCソースが表示されます。経験則として、内側のループには黄色のものは必要ありません。

物事のカップルは、あなたのコードに飛び出す：

• x[j][i]は、各呼び出しでx[j]のための一時的な配列を作成し、その代わりにx[j, i]を使用しています。
• cdef ndarray xの代わりに、次元性とdtype（cdef ndarray[ndim=2, dtype=float]）のどちらかを指定するか、---好ましくは、型付きメモリビュー構文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]と宣言し、最後のインデックスが最も速く変化する配列を反復処理します。

EDIT：型付けmemoryview構文double[:, ::1]など