abc_inference with parallel = true gives error

[ClaudMor]

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

[Vaibhavdixit02]

@claudio20497 it looks like something coming up inside ApproxBayes, the n_particles is not being obtained on some threads and leading to the error from what I understand. I'll pull in @marcjwilliams1 to get this confirmed.

[marcjwilliams1]

Apologies, I never got to the bottom of this. If I'm being honest, I'm unlikely to have time to look in detail for the next month or so, but happy to take a look at a PR if anyone else feels like investigating.