Errors during parallel computation (Threads.@threads)

[Xinyi-Yu]

There are some errors when using Threads.@threads to slove @ivp. The specific codes in test.jl (adopted from page) and errors are as follows,

using ReachabilityAnalysis, Plots

@taylorize function seir2!(du,u,p,t)
  S, E, I, R, α, β, γ = u
  βIS = β * (I * S)
  αE = α*E
  γI = γ*I
  du[1] = -βIS      # dS
  du[2] = βIS - αE  # dE
  du[3] = -γI + αE  # dI
  du[4] = γI        # dR
  local zerou = zero(u[1])
  du[5] = zerou
  du[6] = zerou
  du[7] = zerou
end

E₀ = 1e-4
u₀ = [1-E₀, E₀, 0, 0]
param_error = 0.01
α = 0.2 ± param_error
γ = 0.5 ± param_error


Threads.@threads for i = 1:10
    β = i
    p = [α, β, γ]
    X0 = IntervalBox(vcat(u₀, p));
    prob = @ivp(x' = seir2!(x), dim=7, x(0) ∈ X0);
    sol = solve(prob, tspan=(0.0, 200.0), alg=TMJets21a(orderT=7, orderQ=1));
end

.............$ julia -t 4 test.jl
ERROR: LoadError: TaskFailedException
Stacktrace:
 [1] wait
   @ ./task.jl:334 [inlined]
 [2] threading_run(func::Function)
   @ Base.Threads ./threadingconstructs.jl:38
 [3] top-level scope
   @ ./threadingconstructs.jl:97

    nested task error: AssertionError: order ≤ get_order() && lencoef ≤ num_coeffs
    Stacktrace:
      [1] resize_coeffsHP!(coeffs::Vector{Float64}, order::Int64)
        @ TaylorSeries ~/.julia/packages/TaylorSeries/vTIIX/src/auxiliary.jl:39
      [2] HomogeneousPolynomial
        @ ~/.julia/packages/TaylorSeries/vTIIX/src/constructors.jl:100 [inlined]
      [3] HomogeneousPolynomial
        @ ~/.julia/packages/TaylorSeries/vTIIX/src/constructors.jl:106 [inlined]
      [4] HomogeneousPolynomial(#unused#::Type{Float64}, nv::Int64)
        @ TaylorSeries ~/.julia/packages/TaylorSeries/vTIIX/src/constructors.jl:132
      [5] #TaylorN#27
        @ ~/.julia/packages/TaylorSeries/vTIIX/src/constructors.jl:202 [inlined]
      [6] TaylorN
        @ ~/.julia/packages/TaylorSeries/vTIIX/src/constructors.jl:202 [inlined]
      [7] set_variables(::Type{Float64}, names::Vector{String}; order::Int64)
        @ TaylorSeries ~/.julia/packages/TaylorSeries/vTIIX/src/parameters.jl:146
      [8] set_variables(::Type{Float64}, names::String; order::Int64, numvars::Int64)
        @ TaylorSeries ~/.julia/packages/TaylorSeries/vTIIX/src/parameters.jl:166
      [9] #set_variables#15
        @ ~/.julia/packages/TaylorSeries/vTIIX/src/parameters.jl:171 [inlined]
     [10] post(alg::TMJets21a{Float64, ZonotopeEnclosure}, ivp::InitialValueProblem{BlackBoxContinuousSystem{typeof(seir2!)}, IntervalBox{7, Float64}}, timespan::IntervalArithmetic.Interval{Float64}; Δt0::IntervalArithmetic.Interval{Float64}, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:tspan, :alg), Tuple{Tuple{Float64, Float64}, TMJets21a{Float64, ZonotopeEnclosure}}}})
        @ ReachabilityAnalysis ~/.julia/packages/ReachabilityAnalysis/HmeX4/src/Algorithms/TMJets/TMJets21a/post.jl:29
     [11] solve(::InitialValueProblem{BlackBoxContinuousSystem{typeof(seir2!)}, IntervalBox{7, Float64}}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:tspan, :alg), Tuple{Tuple{Float64, Float64}, TMJets21a{Float64, ZonotopeEnclosure}}}})
        @ ReachabilityAnalysis ~/.julia/packages/ReachabilityAnalysis/HmeX4/src/Continuous/solve.jl:64
     [12] macro expansion
        @ ~/YXY/8.0/case4/test.jl:30 [inlined]
     [13] (::var"#51#threadsfor_fun#1"{UnitRange{Int64}})(onethread::Bool)
        @ Main ./threadingconstructs.jl:85
     [14] (::var"#51#threadsfor_fun#1"{UnitRange{Int64}})()
        @ Main ./threadingconstructs.jl:52
in expression starting at /home/dso/YXY/8.0/case4/test.jl:25

If I delete Threads.@threads or use julia -t 1 test.jl instead of -t 4, it works. If I delete sol = solve(prob, tspan=(0.0, 200.0), alg=TMJets21a(orderT=7, orderQ=1)); or write something else there, it works too.
Any ideas about this problems?

[schillic]

Without checking, I think the TMJets algorithms may be troublesome when used in parallel because the Taylor implementation uses global variables.

[Xinyi-Yu]

I see. Thanks for your explanation! Is there any workaround to solve several @ivps in parallel? OR, is there any possibility for your group to improve this problem in the future?

[mforets]

Did you try _solve_distributed? https://github.com/JuliaReach/ReachabilityAnalysis.jl/blob/master/src/Continuous/solve.jl#L106

[mforets]

See https://github.com/JuliaReach/ReachabilityAnalysis.jl/blob/master/src/Continuous/solve.jl#L46 to use it from the main solve; needs the initial states (X0) to be an array.

[mforets]

Here is a test: https://github.com/JuliaReach/ReachabilityAnalysis.jl/blob/master/test/algorithms/TMJets.jl#L32
Let us know if it worked or not !

[Xinyi-Yu]

I see. Thank you so much! I'll have a try tomorrow!

[Xinyi-Yu]

Hello, @mforets and @schillic. I have tried _solve_distributed. Sometimes it works but sometimes it doesn't, which is the same results as the above codes I wrote.

• If I directly use $ julia -t 4 test.jl, it printed the same error as above;
• If I use $ julia -t 4 and then type codes line by line, the first time I type julia> a = solve(IVP(S, X0s), T=0.1, threading=true) it will print the same errors. However, the second/third/fourth....time I type julia> a = solve(IVP(S, X0s), T=0.1, threading=true), it works.

I have tried on two computers (one is linux and one is windows) and the results are the same. This time, the test.jl is as follows,

using ReachabilityAnalysis, Plots
@taylorize function seir2!(du,u,p,t)
  S, E, I, R, α, β, γ = u
  βIS = β * (I * S)
  αE = α*E
  γI = γ*I
  du[1] = -βIS      # dS
  du[2] = βIS - αE  # dE
  du[3] = -γI + αE  # dI
  du[4] = γI        # dR
  local zerou = zero(u[1])
  du[5] = zerou
  du[6] = zerou
  du[7] = zerou
end
E₀ = 1e-4
u₀ = [1-E₀, E₀, 0, 0]
param_error = 0.01
α = 0.2 ± param_error
γ = 0.5 ± param_error
β = 1
p = [α, β, γ]
X0 = IntervalBox(vcat(u₀, p));
prob = @ivp(x' = seir2!(x), dim=7, x(0) ∈ X0);
X0s = split(X0, [2,1,1,1,1,1,1])
S = system(prob)

a = solve(IVP(S, X0s), T=0.1, threading=true)