用于改进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]等