ReachabilityAnalysis.jl
ReachabilityAnalysis.jl copied to clipboard
Errors during parallel computation (Threads.@threads)
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?
Without checking, I think the TMJets
algorithms may be troublesome when used in parallel because the Taylor implementation uses global variables.
I see. Thanks for your explanation! Is there any workaround to solve several @ivp
s in parallel? OR, is there any possibility for your group to improve this problem in the future?
Did you try _solve_distributed
? https://github.com/JuliaReach/ReachabilityAnalysis.jl/blob/master/src/Continuous/solve.jl#L106
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.
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 !
I see. Thank you so much! I'll have a try tomorrow!
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 typejulia> a = solve(IVP(S, X0s), T=0.1, threading=true)
it will print the same errors. However, the second/third/fourth....time I typejulia> 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)