Hello,
I get a strange error when setting parallel = true in abc_inference.
So I start by opening a julia REPL, enable julia to use multiple cores and launch jupyterlab:
|julia> ENV["JULIA_NUM_THREADS"] = 5
5
|julia> using IJulia
|julia> jupyterlab( dir = "./")
The model is a SIRD with age classes, which compiles and solves correctly:
#jupyterlab(detached = true, dir = "./")
using DifferentialEquations, Distributions, DiffEqBayes
# Model parameters
Ξ² = 0.01# infection rate
Ξ»_R = 0.05 # inverse of transition time from infected to recovered
Ξ»_D = 0.83 # inverse of transition time from infected to dead
iβ = 0.075 # fraction of initial infected people in every age class
π« = vcat([Ξ², Ξ»_R, Ξ»_D]...)
# regional contact matrix and regional population
## regional contact matrix
regional_all_contact_matrix = [3.45536 0.485314 0.506389 0.123002 ; 0.597721 2.11738 0.911374 0.323385 ; 0.906231 1.35041 1.60756 0.67411 ; 0.237902 0.432631 0.726488 0.979258] # 4x4 contact matrix
## regional population stratified by age
N= [723208 , 874150, 1330993, 1411928] # array of 4 elements, each of which representing the absolute amount of population in the corresponding age class.
# Initial conditions
Iβ = repeat([iβ],4)
Sβ = N.-Iβ
Rβ = [0.0 for n in 1:length(N)]
Dβ = [0.0 for n in 1:length(N)]
D_totβ = [0.0 for n in 1:length(N)]
β¬ = vcat([Sβ, Iβ, Rβ, Dβ, D_totβ]...)
# Time
final_time = 20
π― = (1.0,final_time);
function SIRD_ac!(du,u,p,t)
# Parameters to be calibrated
Ξ², Ξ»_R, Ξ»_D = p
# initializa this parameter (death probability stratified by age, taken from literature)
Ξ΄β, Ξ΄β, Ξ΄β, Ξ΄β = [0.003/100, 0.004/100, (0.015+0.030+0.064+0.213+0.718)/(5*100), (2.384+8.466+12.497+1.117)/(4*100)]
Ξ΄ = vcat(repeat([Ξ΄β],1),repeat([Ξ΄β],1),repeat([Ξ΄β],1),repeat([Ξ΄β],4-1-1-1))
C = regional_all_contact_matrix
# State variables
S = @view u[4*0+1:4*1]
I = @view u[4*1+1:4*2]
R = @view u[4*2+1:4*3]
D = @view u[4*3+1:4*4]
D_tot = @view u[4*4+1:4*5]
# Differentials
dS = @view du[4*0+1:4*1]
dI = @view du[4*1+1:4*2]
dR = @view du[4*2+1:4*3]
dD = @view du[4*3+1:4*4]
dD_tot = @view du[4*4+1:4*5]
# Force of infection
Ξ = Ξ²*[sum([C[i,j]*I[j]/N[j] for j in 1:size(C)[1]]) for i in 1:size(C)[2]]
# System of equations
@. dS = -Ξ*S
@. dI = Ξ*S - ((1-Ξ΄)*Ξ»_R + Ξ΄*Ξ»_D)*I
@. dR = Ξ»_R*(1-Ξ΄)*I
@. dD = Ξ»_D*Ξ΄*I
@. dD_tot = dD[1]+dD[2]+dD[3]+dD[4]
end;
#data w.r.t. calibrate
piedmont_cumulative_deaths = [0,2,4,5,5,13,17,21,26,46,59,81,111,133,154,175,209,238,283,315]
# create problem and solve it to check it's allright
problem = ODEProblem(SIRD_ac!, β¬, π―, π«)
solution = solve(problem, saveat = 1:final_time);
If I run abc_inference serially on 1 processor it works:
t = collect(1.:final_time);
priors = [Uniform(0.0,1.0), Uniform(0.0,1.0) , Uniform(0.0,1.0)]
results = abc_inference(problem, Tsit5(), t, piedmont_cumulative_deaths, priors, save_idxs=[4*4+2], Ο΅=0.00001, progress = true)
Preparing to run in serial on 1
Running ABC rejection algorithm... 0%|β | ETA: 0:33:36
processor
Running ABC rejection algorithm...100%|ββββββββββββββββββββββββββββββ| Time: 0:00:05
Preparing to run in serial on 1 processor
ABC SMC population 1, new Ο΅: 612.89...100%|ββββββββββββββββββββββββββββββ| Time: 0:00:01
ABC SMC population 2, new Ο΅: 521.17...100%|ββββββββββββββββββββββββββββββ| Time: 0:00:04
ABC SMC population 3, new Ο΅: 381.8...100%|ββββββββββββββββββββββββββββββ| Time: 0:00:18
ABC SMC population 4, new Ο΅: 358.02...100%|ββββββββββββββββββββββββββββββ| Time: 0:00:55
Total number of simulations: 1.33e+05
Cumulative number of simulations = [1456, 4393, 12297, 43248, 133091]
Acceptance ratio: 3.76e-03
Tolerance schedule = [612.89, 521.17, 381.8, 358.02, 348.86]
Median (95% intervals):
Parameter 1: 0.33 (0.15,0.78)
Parameter 2: 0.65 (0.03,0.98)
Parameter 3: 0.65 (0.00,0.98)
But if I launch the parallelized version, i get an error:
t = collect(1.:final_time);
priors = [Uniform(0.0,1.0), Uniform(0.0,1.0) , Uniform(0.0,1.0)]
results = abc_inference(problem, Tsit5(), t, piedmont_cumulative_deaths, priors, save_idxs=[4*4+2], Ο΅=0.00001, progress = true, parallel = true)
Preparing to run in parallel on 5 processors
Preparing to run in parallel on 5 processors
BoundsError: attempt to access 309-element Array{ApproxBayes.ParticleSMC,1} at index [1:500]
Stacktrace:
[1] throw_boundserror(::Array{ApproxBayes.ParticleSMC,1}, ::Tuple{UnitRange{Int64}}) at .\abstractarray.jl:541
[2] checkbounds at .\abstractarray.jl:506 [inlined]
[3] getindex(::Array{ApproxBayes.ParticleSMC,1}, ::UnitRange{Int64}) at .\array.jl:815
[4] runabc(::ApproxBayes.ABCSMC, ::Array{Int64,1}; verbose::Bool, progress::Bool, parallel::Bool) at C:\Users\claud\.julia\packages\ApproxBayes\RZWxX\src\ABCalgorithm.jl:225
[5] #abc_inference#15 at C:\Users\claud\.julia\packages\DiffEqBayes\DNrGu\src\abc_inference.jl:36 [inlined]
[6] top-level scope at In[8]:3
[7] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091
What am I missing?
Thank you very much