# Ensemble simulations with "paired" trajectories

**URL:** <https://discourse.julialang.org/t/ensemble-simulations-with-paired-trajectories/45194>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [August 19, 2020, 6:06am UTC](https://discourse.julialang.org/t/ensemble-simulations-with-paired-trajectories/45194 "2020-08-19T06:06:45Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Lachlan-Bn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lachlan-bn/32/17008_2.png) [@Lachlan-Bn](https://discourse.julialang.org/u/Lachlan-Bn)\
**Post date:** [August 19, 2020, 6:06am UTC](https://discourse.julialang.org/t/ensemble-simulations-with-paired-trajectories/45194/1 "2020-08-19T06:06:45Z")

</div>

I’m interested in doing something roughly like the following:

1. Choose a (possibly random) initial condition u0
2. From u0, take a (possibly random) “perturbation” of u0, say, x0 = u0 + s
3. Integrate both trajectories, and evaluate some condition at the endpoint of the trajectories (e.g whether they approached the same attractor.)
4. Repeat 1-3 until my endpoint condition is satisfied N times, or maximum iterations is reached.
5. Terminate

I have a vague idea about how I might accomplish this serially (or possibly even with something like `Threads.@threads`), e.g:

```julia
using OrdinaryDiffEq

function eom(dx,x,p,t)
    foo()
end

tspan = (0.0,100.0)
    
N_conditions = 0
N_total = 0
N_maxiters = 100

while N_conditions < N_target && N_total < N_maxiters
    u0 = rand()
    x0 = u0 + rand()
    
    prob = ODEProblem(eom,u0,tspan)
    sol_u0 = solve(prob)

    remake(prob,u0 = x0)
    sol_x0 = solve(prob)

    if check_condition(sol_u0,sol_x0)
        N_conditions += 1
        N_total += 1
    else
        N_total += 1
    end
end

```

But this seems inefficient. I am wondering if this is possible to achieve as some kind of `EnsembleProblem` setup, with a `reduction` that will end the simulation once the target condition is reached. It’s not clear to me what the `prob_func` or `reduction`s should look like in order to achieve this kind of result, so I thought I would ask here.

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [August 19, 2020, 8:24am UTC](https://discourse.julialang.org/t/ensemble-simulations-with-paired-trajectories/45194/2 "2020-08-19T08:24:45Z")

</div>

The parallel integrator of DynamicalSystems.jl offers a framework to evolve multiple initial conditions in parallel at _exactly_ same times

[https://juliadynamics.github.io/DynamicalSystems.jl/dev/advanced/#DynamicalSystemsBase.parallel\_integrator](https://juliadynamics.github.io/DynamicalSystems.jl/dev/advanced/#DynamicalSystemsBase.parallel_integrator)

(it uses DiffEq under the hood)

If you do not care about having both solutions at exactly the same time at all times (e.g. only care that end points match), then `EnsembleProblem` is the way to go.

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [August 19, 2020, 8:59am UTC](https://discourse.julialang.org/t/ensemble-simulations-with-paired-trajectories/45194/3 "2020-08-19T08:59:09Z")

</div>

> [@Lachlan-Bn](#):
>
> Integrate both trajectories, and evaluate some condition at the endpoint of the trajectories (e.g whether they approached the same attractor.)

Just out of curiosity, are you doing this to calculate return times to an attractor, e.g. Poincare return times, to get the distribution of return times?

---

<div class="post-metadata">

**Author:** ![Lachlan-Bn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lachlan-bn/32/17008_2.png) [@Lachlan-Bn](https://discourse.julialang.org/u/Lachlan-Bn)\
**Post date:** [August 19, 2020, 9:18am UTC](https://discourse.julialang.org/t/ensemble-simulations-with-paired-trajectories/45194/4 "2020-08-19T09:18:36Z")

</div>

Thanks for your replies! I have previously played around with DynamicalSystems.jl, but the way I’ve defined my system (and the information I typically want out of the solution) probably lends itself a bit better to a DifferentialEquations setup.

> If you do not care about having both solutions at exactly the same time at all times (e.g. only care that end points match), then `EnsembleProblem` is the way to go.

The only part of the solution I need for this problem is…well, none of it: the quantities I want to look at here are `N_conditions` and `N_total` (i.e, what are they when the simulation stops?) Hence, I’ll only need the endpoints as far as checking the condition.

> Just out of curiosity, are you doing this to calculate return times to an attractor, e.g. Poincare return times, to get the distribution of return times?

No, I’m interested in a fractal dimension related to final-state uncertainty, as in [Redirecting](https://doi.org/10.1016/0375-9601(83)90945-3)

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [August 19, 2020, 9:46am UTC](https://discourse.julialang.org/t/ensemble-simulations-with-paired-trajectories/45194/5 "2020-08-19T09:46:18Z")

</div>

Oh very nice paper, thanks for sharing, I just went through this and I guess it would be nice to have a method that computes the fraction `f` in ChaosTools.jl, if you are interested! c.f. [https://github.com/JuliaDynamics/ChaosTools.jl/issues/121](https://github.com/JuliaDynamics/ChaosTools.jl/issues/121)

The way I imagine it is with either by initializing several integrators in parallel (what EnsembleProblem does) or by doing `reinit!` for every new initial condition (what I typically do in DynamicalSystems.jl, because you typically parallelize at a higher level, e.g. for different parameters).

The user can provide a `terminate` function that checks the state of the integrator and terminates integration on request (e.g. `terminate(integ) = norm(integ.u[3:4]) < 1e-3` a typical condition that terminates when velocity becomes too small). Once terminated, you check the final state of the integrator (always accessible, even if you don’t save any states) and check whether it matches the requested fixed point or not.

Regarding optimization, the function can take a series of sizes `ε` as input, because for `ε2 > ε1` you don’t have to recompute what initial conditions starting in ball `ε1` do; you have already done that in the `ε1` step.

---

<div class="post-metadata">

**Author:** ![Lachlan-Bn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lachlan-bn/32/17008_2.png) [@Lachlan-Bn](https://discourse.julialang.org/u/Lachlan-Bn)\
**Post date:** [August 19, 2020, 10:07am UTC](https://discourse.julialang.org/t/ensemble-simulations-with-paired-trajectories/45194/6 "2020-08-19T10:07:49Z")

</div>

That would be very nice indeed! I’m happy to try to contribute once I get something working for my own case.