# Recording entry and exit time of state space sets with DifferentialEquations.jl
**URL:** https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029
**Category:** Modelling & Simulations
**Tags:** question
**Created:** [September 3, 2020, 10:44pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029 "2020-09-03T22:44:16Z")
**Posts on this page:** 15
**Page:** 1
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [September 3, 2020, 10:44pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/1 "2020-09-03T22:44:16Z")
I am trying to estimate entry and exit times in sets of the state space, as functionality for DynamicalSystems.jl. I have a working implementation for [discrete systems](https://github.com/JuliaDynamics/ChaosTools.jl/pull/123/files#diff-a3440daae39e67f640a670f2c06a7e67) but for continuous it turns out it is much harder. Or hopefully not, perhaps you guys have a better idea.
The figure bellow describes my situation:
I have an ODE that I evolve in time and it gives the continuous black line (in theory) and the blue points in reality. I define some state u\_0, and I need to record the times that the trajectory enters a ball of size ε\_1 centered in u\_0, and also the times that the trajectory exits this ball. The figure shows several u\_0 purely for visualization, in my application only 1 u\_0 will exist.
I also want to do this computation for several ball sizes ε\_1, ε\_2, …, at once, for performance reasons: if the trajectory cannot enter ball of size ε\_1 (e.g it is in state m+1 in the figure), then it absolutely cannot enter ε\_2 and thus the check can be skipped.
Now, I am really not sure how to go about this. Initially I wanted to use the callback functionality of diff eq, to minimize the amount of work I need to do. But I am not sure how to extract t\_1 and t\_2 (and save each in a **different container** , because the first is entry times, the second is exit times).
EDIT HERE:
> The callbacks provide directly the interpolated time
and therefore I can directly store this time in some array. Therefore a simple way to do it (but probably not most performant) is to use `ContinuousCallback`.
---
### Author:  [@jamesjscully](https://discourse.julialang.org/u/jamesjscully)
#### Post date: [September 3, 2020, 11:10pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/2 "2020-09-03T23:10:12Z")
You could use linear interpolation between each u in the trajectory and check to see if the minimum distance between the interpolated line and the center of the ball is equal to the radius of the ball.
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [September 3, 2020, 11:19pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/3 "2020-09-03T23:19:07Z")
This was my original idea. I am a bit afraid because I don’t know whether the integrators can do interpolation like so:
i.e. that there is a loop that could be in principle be interpolated, but falls out of the linear connection.
---
### Author:  [@jamesjscully](https://discourse.julialang.org/u/jamesjscully)
#### Post date: [September 3, 2020, 11:26pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/4 "2020-09-03T23:26:03Z")
The computed trajectory is discrete though, and if the true trajectory has additional features that the integrator doesn’t capture, you just need a better integrator.
---
### Author:  [@jamesjscully](https://discourse.julialang.org/u/jamesjscully)
#### Post date: [September 3, 2020, 11:27pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/5 "2020-09-03T23:27:03Z")
The loop in your picture would be populated by points that you could interpolate between
---
### Author:  [@mforets](https://discourse.julialang.org/u/mforets)
#### Post date: [September 4, 2020, 12:30am UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/6 "2020-09-04T00:30:53Z")
> [@Datseris](#):
>
> Or hopefully not, perhaps you guys have a better idea.
A different approach is to formulate the problem as a hybrid system with modes at the ball’s positions and identity reset maps. Then solve it using set-based integration and count the amount of time spent on each mode. This method is by construction mathematically rigorous and exhaustive (ie. covers all trajectories of the continuous system).
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [September 4, 2020, 12:38am UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/7 "2020-09-04T00:38:57Z")
Although I got to admit that I understood almost nothing of what you said, I should point out that I am interested in the _return_ time to the ε-ball, not the time spent inside it.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 4, 2020, 3:03am UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/8 "2020-09-04T03:03:49Z")
Treating as a hybrid system is using ContinuousCallbacks like we discussed on the #diffeq-bridged channel. You’d just have an upcrossing affect! that fills one array and a downcrossing affect! that fills a different one. That should be a lot cheaper than the interval-based method, but the interval-based method has more of a guarantee (though in this kind of system I wouldn’t see it missing points).
> [@Datseris](#):
>
> This was my original idea. I am a bit afraid because I don’t know whether the integrators can do interpolation like so:
Integrators already handle this: this is exactly the reason for the `interp_points` argument in callbacks, which defaults to 10 meaning it checks at 10 intermediate points of the trajectory for a zero crossing and thus would catch this (and most people never change this argument).
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [September 4, 2020, 7:54pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/9 "2020-09-04T19:54:47Z")
Alright, I am now very pleased to say that I have a working version for a single radius, that goes like this:
```julia
using DynamicalSystemsBase
using OrdinaryDiffEq: Tsit5
using DynamicalSystemsBase.DiffEqBase: ContinuousCallback, ODEProblem, solve
using Distances
εdistance(u, u0, ε::Real) = euclidean(u, u0) - ε
function transit_time_statistics(ds::ContinuousDynamicalSystem, u0, εs, T;
alg = Tsit5(), diffeq...
)
eT = eltype(ds.t0)
check_εs_sorting(εs, length(u0))
exits = [eT[] for _ in 1:length(εs)]
entries = [eT[] for _ in 1:length(εs)]
# Make the magic callback:
crossing(u, t, integ) = ChaosTools.εdistance(u, u0, εs[1])
negative_affect!(integ) = push!(entries[1], integ.t)
positive_affect!(integ) = push!(exits[1], integ.t)
cb = ContinuousCallback(crossing, positive_affect!, negative_affect!;
save_positions = (false, false)
)
prob = ODEProblem(ds, (eT(0), eT(T)); u0 = u0)
sol = solve(prob, alg;
callback=cb, save_everystep = false, dense = false,
save_start=false, save_end = false, diffeq...
)
return exits, entries
end
```
This works so far without problem. You will notice that `εs` is a vector of ball radii, so I need to extend this framework to use a `CallbackSet`. So I change my callback generating code to:
```julia
callbacks = ContinuousCallback[]
for i in eachindex(εs)
crossing(u, t, integ) = ChaosTools.εdistance(u, u0, εs[i])
negative_affect!(integ) = push!(entries[i], integ.t)
positive_affect!(integ) = push!(exits[i], integ.t)
cb = ContinuousCallback(crossing, positive_affect!, negative_affect!;
save_positions = (false, false)
)
push!(callbacks, cb)
end
cb = CallbackSet(callbacks...)
```
which, sometimes, works excellently. There are edge cases that I haven’t been able to identify yet where the integrator and/or interpolation seems to get stuck when exiting the ε-ball and or entering the inner most ε-ball and the callback `positive_affect!` gets called non-stop.
For example if I have `εs = sort!([1.0, 0.1, 0.01]; rev=true)`, I can get resulting in output arrays like (_fully reproducible code at the end_):
```julia
julia> exits[2]
885033-element Array{Float64,1}:
0.009521130925491618
105.1285060629573
227.79608608298875
456.0488702006206
561.549220491485
954.0609666798003
1610.4113223742547
1645.4402235743025
1680.468316963271
1949.7593549776773
2658.4544078395106
⋮
27017.53815060816
27017.53815060816
27017.53815060816
27017.53815060816
27017.53815060816
27017.53815060816
```
while the corresponding `entries[2]`, which, by the definition of continuity of the trajectory, can have at most 1 more entry that `exits[2]`, is:
```julia
julia> entries[2]
80-element Array{Float64,1}:
105.10947207509459
227.7840886145673
456.03680102588226
561.5319922039274
⋮
23077.10816055895
24143.055169575884
24992.380689658865
25009.884596962012
25062.475812273526
26355.939325594318
26508.239300245827
26783.405955508315
27017.519471225707
```
So I guess something happens when exiting one of the inner balls. @ChrisRackauckas do you have any suggestions? In the example above I attempted to integrate for total of `T = 100000.0` so the algorithm shouldn’t stop at `27017.519471225707` ( I guess it hit the max number of steps and thus terminated)
I can also say that I’ve noticed that when this stuck happens, it was the first time that my trajectory crossed the innermost ball of radius `εs[3]`, which has the smallest radius.
* * *
**Reproducible code** (but unfortunatley, needs `fixodprob` branch of `DynamicalSystemsBase`)
```julia
using DynamicalSystemsBase
using OrdinaryDiffEq: Tsit5
using DynamicalSystemsBase.DiffEqBase: ContinuousCallback, ODEProblem, solve
using Distances
function transit_time_statistics(ds::ContinuousDynamicalSystem, u0, εs, T;
alg = Tsit5(), diffeq...
)
eT = eltype(ds.t0)
check_εs_sorting(εs, length(u0))
exits = [eT[] for _ in 1:length(εs)]
entries = [eT[] for _ in 1:length(εs)]
# Make the magic callbacks:
callbacks = ContinuousCallback[]
for i in eachindex(εs)
crossing(u, t, integ) = ChaosTools.εdistance(u, u0, εs[i])
negative_affect!(integ) = push!(entries[i], integ.t)
positive_affect!(integ) = push!(exits[i], integ.t)
cb = ContinuousCallback(crossing, positive_affect!, negative_affect!;
save_positions = (false, false)
)
push!(callbacks, cb)
end
cb = CallbackSet(callbacks...)
prob = ODEProblem(ds, (eT(0), eT(T)); u0 = u0)
sol = solve(prob, alg;
callback=cb, save_everystep = false, dense = false,
save_start=false, save_end = false, diffeq...
)
return exits, entries
end
ro = Systems.roessler(ones(3))
tr = trajectory(ro, 5000; Ttr = 100)
u0 = trajectory(ro, 11; Ttr = 100)[end] # return center
εs = sort!([1.0, 0.1, 0.01]; rev=true)
T = 100000.0 # maximum time
exits, entries = transit_time_statistics(ro, u0, εs, T)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 4, 2020, 7:56pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/10 "2020-09-04T19:56:57Z")
> [@Datseris](#):
>
> So I guess something happens when exiting one of the inner balls. @ChrisRackauckas do you have any suggestions? In the example above I attempted to integrate for total of `T = 100000.0` so the algorithm shouldn’t stop at `27017.519471225707` ( I guess it hit the max number of steps and thus terminated)
It should throw a warning. You might want to increase `maxiters`. You can add that to the `ODEProblem` you generate so it’s higher by default, since yes with that high of a `T` you might accidentally hit it with the defaults.
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [September 4, 2020, 7:58pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/11 "2020-09-04T19:58:29Z")
> [@Datseris](#):
>
> ```julia
> 27017.53815060816
> 27017.53815060816
> 27017.53815060816
> 27017.53815060816
> 27017.53815060816
> 27017.53815060816
>
> ```
That is the problem: the integrator gets stuck at that point and never progresses. But I am not sure why. Increasing `maxiters` doesn’t help here, I clearly shouldn’t be getting the value `27017.53815060816` more than once.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 4, 2020, 8:01pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/12 "2020-09-04T20:01:58Z")
Oh it might be hitting double callback… what DiffEqBase version do you have?
(@kanav99 is fixing some of these issues)
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [September 4, 2020, 8:26pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/13 "2020-09-04T20:26:53Z")
```julia
(@v1.5) pkg> st DiffEqBase
Status `C:\Users\datse\.julia\environments\v1.5\Project.toml`
[2b5f629d] DiffEqBase v6.44.3
```
In the mean time, could you please point me to the point of the source where the interpolation and calling of Roots.jl happens? I want to see if it is possible for me to somehow incorporate this information in my callbacks, so that I let DiffEq handle the first callback (outer most ε ball), but do the inner ones manually. Perhaps this will be both more stable and more performant?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 4, 2020, 8:30pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/14 "2020-09-04T20:30:39Z")
> [@Datseris](#):
>
> `DiffEqBase v6.44.3`
You’ll want to update. This issue is likely from us using `eps(typeof(t))` instead of `eps(t)` which was fixed in [Callback bug fixes by pengwyn · Pull Request #574 · SciML/DiffEqBase.jl · GitHub](https://github.com/SciML/DiffEqBase.jl/pull/574) . That should fix it.
The source for this is all in [https://github.com/SciML/DiffEqBase.jl/blob/master/src/callbacks.jl#L650-L673](https://github.com/SciML/DiffEqBase.jl/blob/master/src/callbacks.jl#L650-L673) . Note that @kanav99 is currently writing a new rootfinder to fix the performance issues of Roots.jl
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [September 4, 2020, 11:05pm UTC](https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029/15 "2020-09-04T23:05:55Z")
Hm, I did update to .47 but the problem remained the same. I’ll clean it up to a pure DiffEqBase MWE and open up an issue.