# Abc_inference with parallel = true gives error in DifferentialEquations.jl

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
``````

What am I missing?

Thank you very much

1 Like

@Vaibhavdixit02 did you look into this when you added the kwarg forwarding?

2 Likes

Replied on the GitHub issue, it looks to be from inside ApproxBayes to me

2 Likes