python - 改进cython代码的高效矩阵向量结构

Question

我必须在 SODE 系统上运行一些模拟。因为我需要使用随机图，所以我认为使用 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 对 cython 源代码进行颜色编码。如果您单击一行，它会显示相应的 C 源代码。根据经验，您不希望内循环中有任何黄色内容。

您的代码中会出现一些问题:

• x[j][i] 在每次调用时为 x[j] 创建一个临时数组，因此使用 x[j, i] 代替。
• 而不是 cdef ndarray x 更好地提供维度和 dtype (cdef ndarray[ndim=2, dtype=float]) 或 --- 最好 ---使用类型化的内存 View 语法: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 顺序排列(最后一个维度变化最快)，将内存 View 声明为 double[:,::1] 并迭代数组，最后一个索引变化最快。

编辑:参见http://cython.readthedocs.io/en/latest/src/userguide/memoryviews.html对于typed memoryview语法double[:,::1]等