# Intersection of \`SciMLBase.EnsembleSolution\`s
**URL:** https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390
**Category:** Numerics
**Tags:** question, package, diffeq
**Created:** [August 31, 2021, 12:32am UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390 "2021-08-31T00:32:32Z")
**Posts on this page:** 11
**Page:** 1
### Author:  [@cadojo](https://discourse.julialang.org/u/cadojo)
#### Post date: [August 31, 2021, 12:32am UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/1 "2021-08-31T00:32:32Z")
I have two `EnsembleSolutions` which represent tubes of orbital trajectories in space (invariant manifolds). I’d like to find the intersection of these two `EnsembleSolutions`. This functionality comes in two levels…
1. Given two `EnsembleSolution`s, find the intersection between them
2. In reality, these `EnsembleSolution`s are in different coordinate frames, so I’d really like to apply a coordinate frame transformation to each solution before finding the intersection
**Is there any off-the-shelf way to find the intersection between two `SciMLBase.EnsembleSolution` instances?** I’m aware of the [K nearest neighbors](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/5) algorithm, and other methods for generally computing the intersection between two multidimensional discrete sets, but I was wondering if there was an off the shelf way which took advantage of the interpolates in each `EnsembleSolution`.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 31, 2021, 12:45am UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/2 "2021-08-31T00:45:20Z")
Why not just use `saveat` so that the intersection is the whole thing?
---
### Author:  [@cadojo](https://discourse.julialang.org/u/cadojo)
#### Post date: [August 31, 2021, 1:09am UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/3 "2021-08-31T01:09:37Z")
I might be misunderstanding, but if I `saveat` I’ll have more discrete data. There’s still no guarantee that there will be an intersection between the two discrete sets of data – in fact, there’s really never an intersection at a discrete data point. I _could_ set `saveat` to something really small, but the vast majority of these manifolds are nowhere near each other.
Doesn’t `saveat` allow a user to save the solution at every `saveat` timestep? Or am I misunderstanding?
Edit: I have a minimum working example below.
```julia
#
# Manifolds!
#
# Environment setup
import Pkg
Pkg.activate(; temp=true)
Pkg.add("Plots")
Pkg.add("Unitful")
Pkg.add("DifferentialEquations")
Pkg.add(; url="https://github.com/cadojo/GeneralAstrodynamics.jl", rev="main")
# Load dependencies
using Plots, Unitful
using DifferentialEquations
using GeneralAstrodynamics
# Find a periodic orbit near Earth (orbit, period)
Oₑ, Pₑ = halo(SunEarth; Az=100_000u"km", L=2)
# Find a periodic orbit near Jupiter (orbit, period)
Oⱼ, Pⱼ = halo(SunJupiter; Az=300_000u"km", L=1)
# Compute an unstable manifold departing Earth
Mₑ = manifold(Oₑ, Pₑ; trajectories=50, saveat=0.01, duration=3Pₑ, eps=1e-6, Trajectory=Val{false}, direction=Val{:unstable})
# Compute a stable manifold arriving at Jupiter
Mⱼ = manifold(Oⱼ, Pⱼ; trajectories=50, saveat=0.01, duration=3Pⱼ, eps=-1e-6, Trajectory=Val{false}, direction=Val{:stable})
# Assume there's some function that takes a state vector,
# and transforms it to a common reference frame
# magic_transformation(state::AbstractVector) = ... # another state vector
# What's the intersection of Mₑ and Mⱼ?
plot(Mₑ, vars=:xy)
plot!(Mⱼ; vars=:xy, palette=:greens)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 31, 2021, 9:54am UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/4 "2021-08-31T09:54:52Z")
> [@cadojo](#):
>
> Doesn’t `saveat` allow a user to save the solution at every `saveat` timestep? Or am I misunderstanding?
No, you can tell the solver to save at pre-defined times and intervals.
---
### Author:  [@cadojo](https://discourse.julialang.org/u/cadojo)
#### Post date: [September 1, 2021, 2:39am UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/5 "2021-09-01T02:39:42Z")
Understood about how `saveat` works. I’m still not understanding how saving at time intervals would help to identify times when two separate solutions intersect. If I don’t know the times of intersection ahead of time, how can I use `saveat` to find the intersection between two solutions?
---
### Author:  [@mforets](https://discourse.julialang.org/u/mforets)
#### Post date: [September 1, 2021, 3:54am UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/6 "2021-09-01T03:54:19Z")
If you don’t want to use flowpipe construction methods / reachability analysis (which implement conservative time discretization, hence including all behaviors between time points), one simple approximation may be to “convexify” the solution sets in Mₑ and Mⱼ, then intersect them in space.
---
### Author:  [@cadojo](https://discourse.julialang.org/u/cadojo)
#### Post date: [September 1, 2021, 4:22am UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/7 "2021-09-01T04:22:33Z")
I was thinking about looking into reachability analysis, especially after seeing the JuliaCon 2021 workshop / content.
I understand that set propagation is scalable for linear systems, and numerical integration is scalable for nonlinear systems. These dynamics are nonlinear, and while I don’t know if they formally fit the definition of “chaotic dynamics”, they are certainly extremely sensitive to initial conditions. I’m a bit concerned that I’ll need a level of precision that set propagation won’t be able to provide (but like I said, I’m definitely going to try, if for no other reason than as an excuse to play around with `ReachabilityAnalysis.jl`).
The point about convexification is really interesting — I wasn’t aware that capability existed / was used. @mforets are you aware of any Julia implementations which “convexify” a discrete set? I searched a bit, but haven’t found any.
---
### Author:  [@mforets](https://discourse.julialang.org/u/mforets)
#### Post date: [September 1, 2021, 7:15am UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/8 "2021-09-01T07:15:03Z")
The discussion about reachability seems like a side topic to this thread, but I couldn’t resist to comment about it on the first place ;).
> are you aware of any Julia implementations which “convexify” a discrete set?
The Julia package to do this sort of thing is LazySets.jl.
I could reproduce your MWE, very nice! Where can I see the dynamic equations you used? For a conceptual proof of the idea in my previous comment, I think it’s necessary to reason by grouping the results of different trajectories on the same time bins. Is that something that can be done using available code?
Let me show a result neglecting time and just focusing on the intersection. The code is available [here as a gist](https://gist.github.com/mforets/a87d294b41e850ee8929144ad2b0c1b9).
The region in magenta encloses the intersection between the manifolds (in the x-y plane). For example, the intersection area is obtained as
```julia
julia> sum(area, Xint)
0.0005889003496004155
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 1, 2021, 10:46pm UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/9 "2021-09-01T22:46:58Z")
> [@cadojo](#):
>
> Understood about how `saveat` works. I’m still not understanding how saving at time intervals would help to identify times when two separate solutions intersect. If I don’t know the times of intersection ahead of time, how can I use `saveat` to find the intersection between two solutions?
Wait, intersection of `sol1.t` and `sol2.t`?
---
### Author:  [@cadojo](https://discourse.julialang.org/u/cadojo)
#### Post date: [September 1, 2021, 11:05pm UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/10 "2021-09-01T23:05:45Z")
> Wait, intersection of `sol1.t` and `sol2.t` ?
No, I’m looking for when the manifolds physically intersect. There’s a unique time for every physical point on the manifold, so I was saying “when do these manifolds intersect” informally. Sorry about that confusion
---
### Author:  [@cadojo](https://discourse.julialang.org/u/cadojo)
#### Post date: [September 1, 2021, 11:08pm UTC](https://discourse.julialang.org/t/intersection-of-scimlbase-ensemblesolution-s/67390/11 "2021-09-01T23:08:45Z")
Ah thanks so much for such a detailed example! The dynamics are available at [`AstrodynamicalModels.jl](https://github.com/cadojo/AstrodynamicalModels.jl). Also, I can definitely set the time bins for the two manifolds to be equivalent.
I’m on mobile at the moment, but I’ll be able to try out this intersection method this weekend. I’ll update here when I do!