Assert in evaluate

[lbenet]

Closes #73

[lbenet]

@mforets Could you please check if the changes made in this branch do not break the integrations? While it is a sanity check, it may break a couple of things...

[mforets]

I just ran the ReachabilityAnalysis tests with this branch, and "TMJets algorithm: linear IVPs" tests fail, see full stacktrace below. (side note: i'm preparing a larger test with all the ARCH-COMP2020 NLN models.. to be added as a standalone script in ReachabilityBenchmarks).

$ julia --color=yes test/runtests.jl
Test Summary:      | Pass  Total
Flowpipe interface |   27     27
Test Summary:                      |
Solution interface: initial states | No tests
Test Summary:                 | Pass  Total
Solution interface: time span |    5      5
Test Summary:       |
Concrete projection | No tests
Test Summary:                            | Pass  Broken  Total
Time-triggered hybrid automaton (HACLD1) |   12       4     16
Test Summary:                  | Pass  Total
Time-triggered solve (EMBrake) |    2      2
Test Summary:                |
Lazy projection of reachsets | No tests
Test Summary:                    |
Concrete projection of reachsets | No tests
Test Summary:                    |
Concrete time-shift of reachsets | No tests
Test Summary:         | Pass  Total
Flowpipe constructors |    9      9
Test Summary: | Pass  Total
INT algorithm |   10     10
Test Summary: | Pass  Total
BOX algorithm |    6      6
Test Summary:    | Pass  Total
GLGM06 algorithm |    5      5
Test Summary:   | Pass  Total
ASB07 algorithm |   18     18
Test Summary:                           | Pass  Total
ASB07 with time-triggered hybrid system |    2      2
Test Summary:    | Pass  Total
TMJets algorithm |    4      4
TMJets algorithm: linear IVPs: Error During Test at /home/mforets/.julia/dev/ReachabilityAnalysis/test/algorithms/TMJets.jl:24
  Got exception outside of a @test
  AssertionError: iscontained(a, tm)
  Stacktrace:
   [1] evaluate at /home/mforets/.julia/dev/TaylorModels/src/evaluate.jl:8 [inlined]
   [2] _broadcast_getindex_evalf at ./broadcast.jl:631 [inlined]
   [3] _broadcast_getindex at ./broadcast.jl:604 [inlined]
   [4] getindex at ./broadcast.jl:564 [inlined]
   [5] macro expansion at ./broadcast.jl:910 [inlined]
   [6] macro expansion at ./simdloop.jl:77 [inlined]
   [7] copyto! at ./broadcast.jl:909 [inlined]
   [8] copyto! at ./broadcast.jl:864 [inlined]
   [9] copy at ./broadcast.jl:840 [inlined]
   [10] materialize at ./broadcast.jl:820 [inlined]
   [11] evaluate(::Array{TaylorModels.TaylorModel1{TaylorN{Float64},Float64},1}, ::IntervalArithmetic.Interval{Float64}) at /home/mforets/.julia/dev/Ta
ylorModels/src/evaluate.jl:23
   [12] overapproximate(::TaylorModelReachSet{Float64}, ::Type{Zonotope}) at /home/mforets/.julia/dev/ReachabilityAnalysis/src/Flowpipes/reachsets.jl:6
35
   [13] _is_intersection_empty(::TaylorModelReachSet{Float64}, ::HalfSpace{Float64,Array{Float64,1}}, ::ReachabilityAnalysis.ZonotopeEnclosure) at /hom
e/mforets/.julia/dev/ReachabilityAnalysis/src/Flowpipes/setops.jl:663
   [14] validated_integ!(::Array{TaylorModelReachSet{Float64},1}, ::Function, ::Array{Float64,1}, ::IntervalBox{1,Float64}, ::Float64, ::Float64, ::Int
64, ::Int64, ::Float64, ::Int64, ::HalfSpace{Float64,Array{Float64,1}}, ::ReachabilityAnalysis.ZonotopeEnclosure, ::Float64, ::Bool, ::Nothing; parse_e
qs::Bool, check_property::Function) at /home/mforets/.julia/dev/ReachabilityAnalysis/src/Algorithms/TMJets/reach.jl:116
   [15] validated_integ! at /home/mforets/.julia/dev/ReachabilityAnalysis/src/Algorithms/TMJets/reach.jl:23 [inlined] (repeats 2 times)
   [16] _validated_integ!(::Array{TaylorModelReachSet{Float64},1}, ::Function, ::Array{Float64,1}, ::IntervalBox{1,Float64}, ::Float64, ::Float64, ::In
t64, ::Int64, ::Float64, ::Int64, ::HalfSpace{Float64,Array{Float64,1}}, ::ReachabilityAnalysis.ZonotopeEnclosure, ::Float64, ::Bool) at /home/mforets/
.julia/dev/ReachabilityAnalysis/src/Algorithms/TMJets/reach.jl:11
   [17] post(::TMJets{Float64,ReachabilityAnalysis.ZonotopeEnclosure}, ::InitialValueProblem{ConstrainedLinearContinuousSystem{Float64,IdentityMultiple
{Float64},HalfSpace{Float64,Array{Float64,1}}},Interval{Float64,IntervalArithmetic.Interval{Float64}}}, ::IntervalArithmetic.Interval{Float64}; time_sh
ift::Float64, kwargs::Base.Iterators.Pairs{Symbol,Tuple{Float64,Float64},Tuple{Symbol},NamedTuple{(:tspan,),Tuple{Tuple{Float64,Float64}}}}) at /home/m
forets/.julia/dev/ReachabilityAnalysis/src/Algorithms/TMJets/post.jl:34
   [18] solve(::InitialValueProblem{ConstrainedLinearContinuousSystem{Float64,IdentityMultiple{Float64},HalfSpace{Float64,Array{Float64,1}}},Interval{F
loat64,IntervalArithmetic.Interval{Float64}}}, ::TMJets{Float64,ReachabilityAnalysis.ZonotopeEnclosure}; kwargs::Base.Iterators.Pairs{Symbol,Tuple{Floa
t64,Float64},Tuple{Symbol},NamedTuple{(:tspan,),Tuple{Tuple{Float64,Float64}}}}) at /home/mforets/.julia/dev/ReachabilityAnalysis/src/Continuous/solve.
jl:61
   [19] top-level scope at /home/mforets/.julia/dev/ReachabilityAnalysis/test/algorithms/TMJets.jl:31
   [20] top-level scope at /home/mforets/Tools/julia-1.4.1-linux-x86_64/julia-1.4.1/share/julia/stdlib/v1.4/Test/src/Test.jl:1113
   [21] top-level scope at /home/mforets/.julia/dev/ReachabilityAnalysis/test/algorithms/TMJets.jl:25
   [22] include(::String) at ./client.jl:439
   [23] top-level scope at /home/mforets/.julia/dev/ReachabilityAnalysis/test/runtests.jl:25
   [24] include(::Module, ::String) at ./Base.jl:377
   [25] exec_options(::Base.JLOptions) at ./client.jl:288
   [26] _start() at ./client.jl:484

Test Summary:                 | Pass  Error  Total
TMJets algorithm: linear IVPs |    1      1      2
ERROR: LoadError: LoadError: Some tests did not pass: 1 passed, 0 failed, 1 errored, 0 broken.
in expression starting at /home/mforets/.julia/dev/ReachabilityAnalysis/test/algorithms/TMJets.jl:24
in expression starting at /home/mforets/.julia/dev/ReachabilityAnalysis/test/runtests.jl:25

[lbenet]

Thanks a lot @mforets.

The assert condition that I included in the function iscontained is defined here. Note that it checks if a is within the re-centered domain, i.e., subtracting the expansion point. Maybe this is the cause of the problems, or something I should correct. The reason I proceeded like this, is because, if you use x0 as the expansion point, the series is written in terms of h=x-x0.

[mforets]

I had a closer look at this, and i think it's just a floating-point issue:

# here R is a Taylor model reach-set 
julia> R
TaylorModelReachSet{Float64}(TaylorModel1{TaylorN{Float64},Float64}[  0.33830061458435673 + 0.037588957176039633 x₁ + ( - 0.33830061458435673 - 0.037588957176039633 x₁) t
+ ( 0.16915030729217836 + 0.018794478588019817 x₁) t² + ( - 0.05638343576405945 - 0.006264826196006605 x₁) t³ + ( 0.014095858941014863 + 0.0015662065490016513 x₁) t⁴ + ( -
 0.0028191717882029726 - 0.00031324130980033026 x₁) t⁵ + ( 0.0004698619647004954 + 5.220688496672171e-5 x₁) t⁶ + ( - 6.712313781435649e-5 - 7.458126423817387e-6 x₁) t⁷ + (
 8.39039222679456e-6 + 9.322658029771734e-7 x₁) t⁸ + [-9.82648e-14, 1.4178e-14]], [0.285312, 0.432405])

julia> t0 = tstart(R)
0.28531269029503514

julia> Δt = tspan(R)
[0.285312, 0.432405]

julia> Δt - t0
[0, 0.147093]

julia> domain(set(R)[1])
[0, 0.147093]

julia> @show Δt - t0 ⊆ domain(set(R)[1])
Δt - t0 ⊆ domain((set(R))[1]) = false
false

This operation is done in functions that overapproximate a Taylor model with other sets such as a zonotope, see this function.

Is it possible to set a tolerance for the containment check? (In LazySets we can set the floating-point tolerance for set operations, see comparisons.jl).

[mforets]

... on the other hand, in that code we can just extract the full domain and pass it to evaluate, saving this difference Δt - t0 which produces the error. In general, this feels a bit restrictive, but if you think this is the way to go, i've opened https://github.com/JuliaReach/ReachabilityAnalysis.jl/pull/189 with the fix on RA.

[mforets]

The reason I proceeded like this, is because, if you use x0 as the expansion point, the series is written in terms of h=x-x0.

IIUC your point was that evaluate(tm, domain(tm)) does not evaluate the tm on its domain. Or it does, but with the definition of shifted domain. I'll think further about this change, maybe it's wrong.

[mforets]

I'll think further about this change, maybe it's wrong.

It is correct, because such expansions in time are always centered at zero (cf. this line, where zI is the zero interval) :relieved:

That said, I think we should nevertheless merge (https://github.com/JuliaReach/ReachabilityAnalysis.jl/pull/189) because that's in some sense a more precise evaluation.

---

Let me try to give a MWE for the issue described in the previous comments.

julia> ii = interval(0.0, 1.3)
[0, 1.30001]

julia> tm = TaylorModel1(2, zero(ii), ii)
 [1, 1] t + [0, 0]

julia> d = domain(tm)
[0, 1.30001]

julia> t0 = inf(d) + 0.05
0.05

julia> Δt = d + 0.05
[0.05, 1.35001]

julia> evaluate(tm, d)
[0, 1.30001]

julia> evaluate(tm, Δt - t0)
ERROR: AssertionError: element [0, 1.30001] is not contained in  [1, 1] t + [0, 0]
Stacktrace:
 [1] evaluate(::TaylorModel1{Interval{Float64},Float64}, ::Interval{Float64}) at /home/mforets/.julia/dev/TaylorModels/src/evaluate.jl:8
 [2] top-level scope at REPL[84]:1

julia> inf(d)
0.0

julia> inf(Δt - t0)
0.0

julia> inf(Δt - t0) ⊆ d
true

julia> sup(Δt - t0) ⊆ d
false

julia> sup(Δt - t0)
1.3000000000000003

[lbenet]

Let me try to explain why I did it the way I did it.

Usually, when expanding a function f(t) around t0 we write: f(t) = f_0 + f_1 (t-t0) + f_2 (t-t_0)^2 + ... (1) which we usually write as: f(t0 + h) = f_0 + f_1 h + f_2 h^2 + ... (2) with h = t-t0. Equation (2) corresponds to evaluate(f, h) or f(h) in TaylorSeries.

The evaluation of f(t1) therefore corresponds to f(t0 + (t1-t0)), or evaluate(f, t1-t0), which is the reason we always subtract the expansion point. That is the reason that for the evaluation of the TaylorModels I check if t1 is inside D-t0.

[lbenet]

As for your example, if you want to evaluate tm at t0, you need to evaluate(tm, t0-tm.x0):

julia> evaluate(tm, t0-tm.x0)
[0.05, 0.0500001]

[mforets]

Yes, that makes total sense and i'm aware it's the convention used in this package :+1:

The only thing the MWE above shows is that, because of rounding errors, evaluate(tm, Δt - t0) doesn't work as before even with a "tiny" difference in the domain. In the reachability application, we keep track of consecutive time-spans like this: [0, t1], [t1, t2], [t2, t3], ... so sometimes you want to evaluate the taylor model reach-set in its whole domain and we were using arithmetic operations on the [ti, tj] to deduce the domain... this fails now, because of rounding errors. In fact, that's just:

julia> using IntervalArithmetic

julia> d = 0 .. 1.3
[0, 1.30001]

julia> sup((d + 0.05) - (inf(d) + 0.05)) ⊆ d
false

Of course, this is the way TaylorModels + interval domains is expected to work.. so i'm fine with this change, i'll fix what needs to be fixed in the other package :smile:

[lbenet]

@dpsanders I think this is ready for review.

[lbenet]

Merging...