我必须在 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
最佳答案
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
对于类型化内存 View 语法double[:, ::1]
等等
关于python - 用于改进 cython 代码的高效矩阵向量结构,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/39789262/