I have been using Julia to simulate calculations involving directed networks. Monte Carlo simulations(100 iterations) with a network with 100 nodes used to take 10mins. But the code doesn’t scale well, for a directed network with 500 nodes the same set of operations take 90 minutes. I am a new migrant to Julia, could anyone suggest how to fix such scalability issues? My aim is to be able to use the same set of manipulations on data obtained from real networks.
Can you supply a piece of code that shows this behaviour?
It could be a memory issue. If you use the package BenchmarkTools.jl to run your main function it should show you how much memory was allocated. Perhaps for 500 nodes this value is larger than your physical memory?
1 Like
I tried to optimize it in all ways possible (that I know of). The code is below, and any inputs on how to improve the same are welcome.
using Pkg
using Distributed
using BenchmarkTools
using JLD
using Plots
@everywhere using Distributions
@everywhere using StatsBase
@everywhere using LightGraphs
@everywhere using LinearAlgebra
####################################################################################################################################################
@everywhere mutable struct AllIndParam
N::Int64 ;
Nᵣ::Int64 ;
Nₗ::Int64 ;
μ::Int64 ;
σ::Int64 ;
choice::Int64 ;
num1::Int64 ;
num2::Int64 ;
Cmax::Int64 ;
K::Int64 ;
h::Int64 ;
tf::Int64 ;
disconnected::Int64 ;
loopin::Int64 ;
mc::Int64 ;
mcf::Int64 ;
rᵣ::Float64 ;
rₗ::Float64 ;
P::Float64 ;
Pᵢ::Float64 ;
Pᵣ::Float64 ;
vs_avg::Float64 ;
vc_avg::Float64 ;
rowsum::Float64 ;
g::SimpleDiGraph ;
edg::Array{Int64,1} ;
cedg::Array{Int64,1} ;
unconn::Array{Int64,1} ;
temp1::Array{Int64,1} ;
end
#########################################################################################################################################################
@everywhere mutable struct AllDepParam
n::Int64 ;
agents::Array{Int64,1} ;
con::Array{Int64,1} ;
lib::Array{Int64,1} ;
o::Array{Int64,1} ;
upagents::Array{Int64,1} ;
upcount::Array{Int64,1} ;
Sthresh::Array{Int64,1} ;
R::Array{Float64,1} ;
θ::Array{Float64,1} ;
θtemp::Array{Float64,1} ;
C::Array{Float64,1} ;
Asep::Array{Int64,2} ;
A::Array{Int64,2} ;
Ain::Array{Int64,2} ;
allag::Array{Int64,2} ;
NC::Array{Int64,2} ;
score::Array{Float64,2} ;
W::Array{Float64,2} ;
end
################################################################################################################################################################
@everywhere mutable struct AllTimParam
Maxc::Array{Int64,1} ;
uptime::Array{Int64,2} ;
θₜ::Array{Float64,2} ;
conn::Array{Bool,1} ;
end
#################################################################################
#(N,Nᵣ,Nₗ,μ,σ,choice,num1,num2,Cmax,K,h,tf,disconnected,loopin,mc,mcf,rᵣ,rₗ,P,Pᵢ,Pᵣ,vs_avg,vc_avg,rowsum,g,edg,cedg,unconn,temp1)
@everywhere ip = AllIndParam(1000,500,500,90,10,0,0,0,0,1000,1,1500,1,1,1,1,round(10*pi/180,digits=4),round(40*pi/180,digits=4),0.5,0.5,0.25,0.0,
0.0,0.0,SimpleDiGraph(Int64),[],[],[],[]) ;
#(n,agents,con,lib,o,upagents,upcount,Sthresh,R,θ,θtemp,C,Asep,A,Ain,allag,NC,temp,score,W)
@everywhere dp = AllDepParam(Int64((ip.tf + 1)/ip.h),collect(1:1:ip.N),zeros(Int64,ip.Nᵣ),zeros(Int64,ip.Nₗ),
zeros(Int64,ip.N),[],zeros(Int64,ip.N),zeros(Int64,ip.N),zeros(Float64,ip.N),zeros(Float64,ip.N),zeros(Float64,ip.N),
zeros(Float64,ip.N),zeros(Int64,ip.N,ip.N),zeros(Int64,ip.N,ip.N),zeros(Int64,ip.N,ip.N),zeros(Int64,ip.N,ip.N),
zeros(Int64,ip.N,ip.N),zeros(Float64,ip.N,ip.N),zeros(Float64,ip.N,ip.N)) ;
#(Maxc,uptime,θₜ,conn)
@everywhere tp = AllTimParam(zeros(Int64,dp.n),zeros(Int64,dp.n,ip.N),zeros(Float64,dp.n,ip.N),zeros(Bool,dp.n)) ;
#################################################################################
@everywhere function fixrigid(ip::AllIndParam,dp::AllDepParam)
ip.choice == 1 ? dp.con = sortperm(dp.θ,rev=true)[1:ip.Nᵣ] :
ip.choice == 2 ? dp.con = sortperm(dp.θ)[1:ip.Nᵣ] :
ip.choice == 3 ? dp.con = intersect(findall(dp.θ.<((ip.μ+ip.σ)*pi/180)),findall(dp.θ.>((ip.μ-ip.σ)*pi/180)))[1:ip.Nᵣ] :
ip.choice == 4 ? (rigid[1:ip.num1] = sortperm(dp.θ,rev=true)[1:ip.num1] ; rigid[ip.num1+1:ip.Nᵣ] = sortperm(dp.θ)[ip.num1+1:ip.Nᵣ] ;) :
dp.con = sample(dp.agents, ip.Nᵣ, replace = false) ;
end
#################################################################################
@everywhere function initial(ip::AllIndParam,dp::AllDepParam)
dp.θ = round.(rand(Normal(ip.μ,ip.σ),ip.N)*pi/180,digits=4) ;
fixrigid(ip,dp) ;
dp.lib = filter!(e->e∉dp.con,dp.agents) ;
dp.R[dp.con] .= ip.rᵣ ;
dp.R[dp.lib] .= ip.rₗ ;
end
#################################################################################
@everywhere function separation(ip::AllIndParam,dp::AllDepParam)
@inbounds for i = 1:1:ip.N
for j = 1:1:ip.N
(dp.θ[i] - dp.θ[j])>-dp.R[i] && (dp.θ[i] - dp.θ[j])<dp.R[i] && (i!=j) ? dp.Asep[i,j] = 1 : dp.Asep[i,j] = 0 ;
end
end
end
##################################################################################
@everywhere function correctdist(ip::AllIndParam,dp::AllDepParam)
ip.disconnected = 1 ;
while ip.disconnected==1
initial(ip,dp) ;
separation(ip,dp) ;
is_connected(SimpleDiGraph(dp.Asep)) && length(findall(sum(dp.Asep,dims=2).==0))==0 ? ip.disconnected = 0 : ip.disconnected = 1 ;
end
end
#################################################################################
@everywhere function desiredgraph(ip::AllIndParam,dp::AllDepParam)
ip.disconnected = 1 ;
while ip.disconnected==1
@inbounds for i = 1:1:ip.N
ip.edg = findall(dp.Asep[i,:].!=0) ;
ip.cedg = findall(rand(length(ip.edg)).<ip.Pᵣ) ;
dp.A[i,ip.edg[ip.cedg]].= 1 ;
end
ip.unconn = findall(sum(dp.A,dims=2).==0) ;
@inbounds for j = 1:1:length(ip.unconn)
ip.edg = findall(dp.Asep[ip.unconn[j],:].!=0) ;
ip.loopin = 1 ;
while (ip.loopin==1)
ip.cedg = findall(rand(length(ip.edg)).<ip.Pᵣ) ;
length(ip.cedg)>0 ? ip.loopin = 0 : ip.loopin = 1 ;
end
dp.A[ip.unconn[j],ip.edg[ip.cedg]] .= 1 ;
end
ip.g = SimpleDiGraph(dp.A) ;
is_connected(ip.g) ? ip.disconnected = 0 : ip.disconnected = 1 ;
end
dp.o = reshape(sum(dp.A,dims=2),ip.N) ;
end
#################################################################################
@everywhere function updatingagents(ip::AllIndParam,dp::AllDepParam)
dp.upagents = setdiff(findall(rand(ip.N).<ip.P),findall(dp.o.==0)) ;
dp.upcount[dp.upagents] = dp.upcount[dp.upagents].+ 1 ;
end
#################################################################################
@everywhere function kcentrality(ip::AllIndParam,dp::AllDepParam)
dp.C = broadcast(abs,katz_centrality(ip.g)) ;
dp.C = dp.C/sum(dp.C) ;
ip.Cmax = argmax(dp.C) ;
end
#################################################################################
@everywhere function initialweights(ip::AllIndParam,dp::AllDepParam)
@inbounds for i = 1:1:ip.N
for j = 1:1:ip.N
dp.A[i,j] == 1 && i!=j ? dp.W[i,j] = rand(1)[1] : dp.W[i,j] = 0 ;
end
ip.rowsum = sum(dp.W[i,:]) ;
ip.rowsum!=0 ? dp.W[i,:] = (dp.W[i,:]*(1-dp.C[i]))/ip.rowsum : dp.W[i,:] = dp.W[i,:] ;
dp.W[i,i] = dp.C[i] ;
end
end
#################################################################################
@everywhere function ttneighbors(ip::AllIndParam,dp::AllDepParam)
Twhop = dp.A^2 ;
Twhop[diagind(Twhop)] .= 0 ;
Thhop = dp.A^3 ;
Thhop[diagind(Thhop)] .= 0 ;
dp.W[diagind(dp.W)] .= 0 ;
TwWt = dp.W^2 ;
TwWt[diagind(TwWt)] .= 0 ;
ThWt = dp.W^3 ;
ThWt[diagind(ThWt)] .= 0 ;
@inbounds for i = 1:1:length(dp.upagents)
for j = 1:1:ip.N
if(dp.A[dp.upagents[i],j]==0 && dp.Asep[dp.upagents[i],j]==1 && rand(1)[1]<ip.Pᵢ)
if(Twhop[dp.upagents[i],j]!=0 && Thhop[dp.upagents[i],j]==0)
dp.allag[dp.upagents[i],j] = j ;
dp.score[dp.upagents[i],j] = TwWt[dp.upagents[i],j]/Twhop[dp.upagents[i],j] ;
dp.NC[dp.upagents[i],j] = dp.NC[dp.upagents[i],j] + 1 ;
elseif(Twhop[dp.upagents[i],j]==0 && Thhop[dp.upagents[i],j]!=0)
dp.allag[dp.upagents[i],j] = j ;
dp.score[dp.upagents[i],j] = ThWt[dp.upagents[i],j]/Thhop[dp.upagents[i],j] ;
dp.NC[dp.upagents[i],j] = dp.NC[dp.upagents[i],j] + 1 ;
elseif(Twhop[dp.upagents[i],j]!=0 && Thhop[dp.upagents[i],j]!=0)
dp.allag[dp.upagents[i],j] = j ;
dp.score[dp.upagents[i],j] = (TwWt[dp.upagents[i],j]+ThWt[dp.upagents[i],j])/(Twhop[dp.upagents[i],j]+Thhop[dp.upagents[i],j]) ;
dp.NC[dp.upagents[i],j] = dp.NC[dp.upagents[i],j] + 1 ;
end
elseif(dp.A[dp.upagents[i],j]==1 && dp.upagents[i]!=j)
dp.allag[dp.upagents[i],j] = j ;
dp.score[dp.upagents[i],j] = dp.W[dp.upagents[i],j] ;
dp.NC[dp.upagents[i],j] = 0 ;
end
end
dp.score[dp.upagents[i],:] = (dp.score[dp.upagents[i],:]*(1-dp.C[dp.upagents[i]]))/sum(dp.score[dp.upagents[i],:]) ;
dp.score[dp.upagents[i],dp.upagents[i]] = dp.C[dp.upagents[i]] ;
dp.allag[dp.upagents[i],dp.upagents[i]] = dp.upagents[i] ;
dp.W[dp.upagents[i],dp.upagents[i]] = dp.C[dp.upagents[i]] ;
end
end
#################################################################################
@everywhere function asyncupdate(ip::AllIndParam,dp::AllDepParam)
dp.θtemp = copy(dp.θ) ;
dp.Asep[diagind(dp.Asep)] .= 1 ;
@inbounds for i = 1:1:length(dp.upagents)
ip.vs_avg = 0.0 ;
ip.vc_avg = 0.0 ;
ip.temp1 = findall(dp.allag[dp.upagents[i],:].!=0) ;
for j = 1:1:length(ip.temp1)
if(dp.Asep[dp.upagents[i],ip.temp1[j]]==1)
ip.vc_avg = ip.vc_avg + dp.score[dp.upagents[i],ip.temp1[j]]*cos(dp.θ[ip.temp1[j]]) ;
ip.vs_avg = ip.vs_avg + dp.score[dp.upagents[i],ip.temp1[j]]*sin(dp.θ[ip.temp1[j]]) ;
end
end
dp.θtemp[dp.upagents[i]] = round(atan(ip.vs_avg,ip.vc_avg),digits=4) ;
end
dp.Asep[diagind(dp.Asep)] .= 0 ;
dp.θ = copy(dp.θtemp) ;
end
#################################################################################
@everywhere function networkevol(ip::AllIndParam,dp::AllDepParam)
dp.Sthresh = round.(ip.K*dp.C) ;
@inbounds for i = 1:1:ip.N
for j = 1:1:ip.N
if(dp.allag[i,j]!=0 && (dp.NC[i,j]>dp.Sthresh[i]))
dp.NC[i,j] = 0 ;
dp.A[i,j] = 1 ;
elseif(dp.Asep[i,j]==0 && dp.A[i,j]==1)
dp.A[i,j] = 0 ;
end
end
end
end
#################################################################################
@everywhere function directweights(ip::AllIndParam,dp::AllDepParam)
@inbounds for i = 1:1:ip.N
for j = 1:1:ip.N
if(dp.A[i,j]!=0)
if(dp.W[i,j]!=0)
dp.W[i,j] = dp.W[i,j] ;
elseif(dp.W[i,j]==0)
dp.W[i,j] = dp.score[i,j] ;
end
end
end
ip.rowsum = sum(dp.W[i,:]) ;
if(ip.rowsum!=0)
dp.W[i,:] = (dp.W[i,:]*(1-dp.C[i]))/ip.rowsum ;
end
dp.W[i,i] = dp.C[i] ;
end
end
#################################################################################
@everywhere function MonteCarlo(ip::AllIndParam,dp::AllDepParam,tp::AllTimParam)
correctdist(ip,dp) ;
desiredgraph(ip,dp) ;
kcentrality(ip,dp) ;
initialweights(ip,dp) ;
@inbounds for t = 1:1:dp.n
tp.Maxc[t] = ip.Cmax ;
tp.θₜ[t,:] = dp.θ' ;
updatingagents(ip,dp) ;
tp.uptime[t,:] = dp.upcount' ;
ttneighbors(ip,dp) ;
asyncupdate(ip,dp) ;
separation(ip,dp) ;
networkevol(ip,dp) ;
ip.g = SimpleDiGraph(dp.A) ;
dp.o = reshape(sum(dp.A,dims=2),ip.N) ;
kcentrality(ip,dp) ;
directweights(ip,dp) ;
end
@show tp.conn[t] = is_connected(ip.g) ;
end
#################################################################################
g::SimpleDiGraph ;
Make this g::SimpleDiGraph{Int64} or g::SimpleDiGraph{Int32}. Otherwise the type is not concrete.
I recommend optimizing a single-threaded / sequential implementation first, and I recommend taking a look at Profile.@profile before posting here, and isolating the core issues of your problem. Then you are more likely to receive meaningful help.
1 Like
I have tried to use @profile, but as a beginner I find the output of the Profile.print() extremely hard to understand. I face the same issue with the benchmark tools. I have written the code after referring to several Julia resources for best coding practises.
Note that both Juno and ProfileView.jl can help you visualize your profiling results and make them more accessible.
1 Like
I apologize for maybe sounding annoyed. In order to get better help here, I suggest you refactor a little. First, remove all @everywhere / multithreading / distributed annotations; single threaded performance needs to be optimized first. Second, you need to refactor the code you post here in order to allow people to help you. Optimally, you would post two code sections: One that defines utility functions (the actual meat), and a second one that drives you computation. That would look like
test_data = #.... generate random test data with similar performance characteristics as your real problem
@btime MonteCarlo(test_data)
#outbut of btime, you must post this!
If possible, also include some
function summarize_results(result)
...
end
that pretty prints 3 lines of summary statistics of the result of you computation. This is such that people can quickly check whether their modifications still produce sensible results for some tiny examples. Next, you also explain what you want to compute in one paragraph, e.g. “I want to simulate the process described in [arxiv link].” or whatever you want.
Then people can look at your MonteCarlo function, try one of the steps after the other, and propose simple modifications. People get often 100 x improvements within a couple of hours of this process, if they spend the efffort of formatting their question right. This is often type stability plus complexity class improvement, due to caching some values instead of recomputing them in a hot loop. The way your code is currently structured, it is very hard to read (which means that few people will bother, and at most skim for obvious issues like the non-concrete graph type I mentioned).
1 Like