add coupledSDEs type

[oameye]

resolves #202

This PR adds a new CoupledSDEs struct. Here, we follow closely the design of CoupledODEs:

struct CoupledSDEs{IIP,D,I,P,S} <: ContinuousTimeDynamicalSystem
    # D parametrised by length of u0
    # S indicated if the noise strength has been added to the diffusion function
    integ::I
    # things we can't recover from `integ`
    p0::P
    noise_strength
    diffeq # isn't parameterized because it is only used for display
end

Feel free to look at the documentation of CT to understand the API. The type should capture equation of the form $dx = f(x, p) dt +g(x, p) dW$. The function g both captures the noise strength $\sigma$ and the covariance matrix (correlation). In general, g can be quite an exotic function, making it that we lose both the noise strength and covariance matrix. However, we need this information to use algorithms in, e.g., large deviation theory. To solve this, I have opted for two design choses:

• During the construction of the CoupledSDEs the user can opt to give a noise strength as an argument. This noise strength will be saved inside the CoupledSDEs struct. However, during the creation the g wil also be rescaled with the noise strength. The parameter S in CoupledSDEs{IIP,D,I,P,S} indicated if this has happened or not.
• I have merged a PR in DiffEqNoiseProcess.jl that adds a covariance field to the noise process dW. This allows the user to add a correlated noise process like CorrelatedWienerProcess where the covariance matrix is automatically saved in the Integrator.

This solved the problem of not having acces to this information after the construction. Nevertheless, there is still some ambiguity as the user can always choose the g function freely. We can put the responsibility with the user and make tests that try to guess the type of noise. Or we will have to make subtypes such that we can easily distinguish scalar, multiplicative, additive, etc.

---

TODO

• Make add covariance kwarg and giving g function optional
• remove noise strength
• We have to check if for SDEIntegrators a different step ReturnCode is possible.
• add benchmark with CoupledSDEs and SDEProblem
• Think about naming kwargs and fields

[Datseris]

This is looking like a really good base, but there are two key things missing:

1. You need to port over here in this PR the docs from CrtiicalTransitions.jl. Just copy the CoupledSDEs.md here and add it as one more doc page. Given that this is a rather complex new type, it is fine for it to have its own doc page.
2. We need to try this out also with another downstream package, just to make sure that there isn't some assumption we haven't thought of but we "bake in" this new type. What I would recommend is to find attractors using featurizing: if you can apply this simple code snippet https://juliadynamics.github.io/DynamicalSystemsDocs.jl/attractors/stable/examples/#Trivial-featurizing-and-grouping-for-basins-fractions to any sort of stochastic bistable system, and see that there are two attractors, tjhat would be absolutely fantastic. You can just provide the final script jere as a comment, and when this PR is merged I'll add it to Attractors.jl as an example. cc @reykboerner

A bit off topic but https://github.com/JuliaDynamics/CriticalTransitions.jl/blob/ce6d1b19e95eab772b71b158f35f64ee976484d9/src/utils.jl#L12-L14 needs to be reworked to return a static vector. If this is added here I will show you how with a comment.

[Datseris]

We can put the responsibility with the user and make tests that try to guess the type of noise. Or we will have to make subtypes such that we can easily distinguish scalar, multiplicative, additive, etc.

Why do we have to care about this? Is there any algorithm that cares about this information and would need to access it?

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 92.85714% with 11 lines in your changes missing coverage. Please review.

Project coverage is 87.59%. Comparing base (6c98239) to head (de8258e).
Report is 30 commits behind head on main.

Files with missing lines	Patch %	Lines
ext/src/CoupledSDEs.jl	86.41%	11 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #212      +/-   ##
==========================================
+ Coverage   82.00%   87.59%   +5.58%     
==========================================
  Files          15       18       +3     
  Lines         717      927     +210     
==========================================
+ Hits          588      812     +224     
+ Misses        129      115      -14     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[oameye]

e can put the responsibility with the user and make tests that try to guess the type of noise. Or we will have to make subtypes such that we can easily distinguish scalar, multiplicative, additive, etc.

Why do we have to care about this? Is there any algorithm that cares about this information and would need to access it?

Yes, for example, the large deviation algorithms implemented in CT.jl are not applicable to all types of stochastic differentials. For further details, I would rather have @reykboerner, who possesses much more expertise on this topic than I do.

[oameye]

The documentation does not build because of this error:

┌ Error: duplicate docs found for 'CoupledODEs' in `@docs` block in src/index.md:54-56
│ ```@docs
│ CoupledODEs
│ ```

Not sure why?

[reykboerner]

Hi both, yes, as @oameye said, the information on the noise term is important for determining which methods can be applied to a CoupledSDEs system. In DiffEq, the noise function g(u,p,t) can be defined very generally: it can be nonlinear, be state-dependent and time-dependent. However, several methods from large deviation theory require a linear combination of noise processes given by a covariance matrix, which often needs to be invertible and is often assumed constant (i.e. state-independent). There are generalizations to the state-dependent (multiplicative noise) case but they require some care when coding them.

I think the most important criterion to classify a CoupledSDEs would be whether the noise term can be expressed in terms of an invertible covariance matrix or not. We could also consider making a subtype of CoupledSDEs, e.g. StandardCoupledSDEs that fulfills this condition, so that the relevant algorithms only work on that subtype.

Alternatively, we could add warnings and leave it to the user to ensure that they are using an appropriate method. @raphael-roemer and @ryandeeley might have input, too

[Datseris]

There are two options I see: option 1: request the user to specify the noise type by in addition providing a symbol that indicates the noise type such as :white or :correlated or whatever you guys think makes sense. This would be done via a keyword whose defualt value would be :unspecified. 2: attempt to deduce it with the user provided input. I don't know how hard option 2 is. We can have some heuristics that trigger if the default value :unspecified is given.

The noise type would be one more field of the struct that can be examined by other functions.

[Datseris]

The documentation does not build because of this error:

┌ Error: duplicate docs found for 'CoupledODEs' in `@docs` block in src/index.md:54-56
│ ```@docs
│ CoupledODEs
│ ```

Not sure why?

You can't expand a docstring twice @oameye , because then when you refer to it documenter.jl doesn't know where to hyperlink. Sinmply remove the @docs block with the ODEs, as this docstring is already expanded in the index page.

Alternatively, if you really want this docstring to be tehre as well, put ; canonical = false next to the ```@docs line

[oameye]

Here is the example with attractors.jl.

using Attractors, StaticArrays
using DynamicalSystemsBase

using DynamicalSystemsBase: idfunc

function fitzhugh_nagumo(u, p, t)
    x, y = u
    ϵ, β, α, γ, κ, I = p

    dx = (-α * x^3 + γ * x - κ * y + I) / ϵ
    dy = -β * y + x

    return SA[dx, dy]
end

sde = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), [1.,3.,1.,1.,1.,0.], 0.1)
ode = CoupledODEs(sde)

function featurizer(X, t)
    return X[end]
end

mapper = AttractorsViaFeaturizing(ode, featurizer; Ttr = 200, T = 1)

xg = yg = range(-4, 4; length = 101)

region = HRectangle([-4, 4], [1, 1])
sampler, _ = statespace_sampler(region)

fs = basins_fractions(mapper, sampler)

However, should the AttractorsViaFeaturizing work when the sde has much noise? Now I convert the sde to an ode. The trajectory solver is much faster so in general people will want convert the SDEIntegrator to an ODEIntegrator`.

[Datseris]

No I want to see the SDE system being given to AttractorsViaFeaturizing directly to ensure the interfaces work. Just put a very small amount of noise to ensure that the multistability remains intact.

[oameye]

This code snippet works, but takes a long time to complete (~1 min). Lowering the default tolerances to 1e-3 helps.

using Attractors, StaticArrays
using DynamicalSystemsBase

using DynamicalSystemsBase: idfunc

function fitzhugh_nagumo(u, p, t)
    x, y = u
    ϵ, β, α, γ, κ, I = p

    dx = (-α * x^3 + γ * x - κ * y + I) / ϵ
    dy = -β * y + x

    return SA[dx, dy]
end

sde = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), [1.,3.,1.,1.,1.,0.], 0.01)

function featurizer(X, t)
    return X[end]
end

mapper = AttractorsViaFeaturizing(sde, featurizer; Ttr = 10, T = 1)

xg = yg = range(-1, 1; length = 51)

region = HRectangle([-1, 1], [1, 1])
sampler, _ = statespace_sampler(region)

fs = basins_fractions(mapper, sampler)

[oameye]

You can't expand a docstring twice @oameye , because then when you refer to it documenter.jl doesn't know where to hyperlink. Sinmply remove the @docs block with the ODEs, as this docstring is already expanded in the index page.

Thanks @Datseris. Specifying which docstring I wanted did the trick.

[oameye]

Regarding the Attractors example. I looked up the default tolerances by StochasticDiffEq. Apperently they only use 1e-2, see here. So that's why it took so long. Will change the default abstol and relto to 1e-2.

[Datseris]

This code snippet works, but takes a long time to complete (~1 min). Lowering the default tolerances to 1e-3 helps.

That's surprising. Which tolerance are you referring to? There is no tolerance in the code snippet you pasted. By default basin fractions will integrate 1,000 initial conditions, and then cluster them. I assume 99.9% of the time here is spent on integrating the initial conditions. Is solving a 2D SDE really so costly that you need about 1 milisecond per initial condition?

[Datseris]

SA[dx, dy]

Can you try changing this to SVector(dx, dy) and see if it makes things faster? I've never seen such a syntax before.

[Datseris]

In any case, it is absolutely fantastic that this just works out of the box with Attractors.jl.

[oameye]

Regarding the Attractors example. I looked up the default tolerances by StochasticDiffEq. Apperently they only use 1e-2, see here. So that's why it took so long. Will change the default abstol and relto to 1e-2.

This is what solved it. SDE solvers require much lower tolerances.
This snippet runs in miliseconds with the new set tolerances.

using Attractors, StaticArrays
using DynamicalSystemsBase

using DynamicalSystemsBase: idfunc

function fitzhugh_nagumo(u, p, t)
    x, y = u
    ϵ, β, α, γ, κ, I = p

    dx = (-α * x^3 + γ * x - κ * y + I) / ϵ
    dy = -β * y + x

    return SA[dx, dy]
end

sde = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), [1.,3.,1.,1.,1.,0.], 0.01)

function featurizer(X, t)
    return X[end]
end

mapper = AttractorsViaFeaturizing(sde, featurizer; Ttr = 200, T = 1)

xg = yg = range(-1, 1; length = 101)

region = HRectangle([-1, 1], [1, 1])
sampler, _ = statespace_sampler(region)

fs = basins_fractions(mapper, sampler)

Can you try changing this to SVector(dx, dy) and see if it makes things faster? I've never seen such a syntax before.

It doesn't make a difference. It is the same.

In any case, it is absolutely fantastic that this just works out of the box with Attractors.jl.

Yeah, I agree. You laid out a really good foundation :)

[Datseris]

Just to be clear as I haven't run the code block: basins_fractions finds two attractors, right?

[oameye]

Just to be clear as I haven't run the code block: basins_fractions finds two attractors, right?

Yes :) Output is:

Dict{Int64, Float64} with 2 entries:
  2 => 0.353
  1 => 0.647

There are two options I see: option 1: request the user to specify the noise type by in addition providing a symbol that indicates the noise type such as :white or :correlated or whatever you guys think makes sense. This would be done via a keyword whose defualt value would be :unspecified. 2: attempt to deduce it with the user provided input. I don't know how hard option 2 is. We can have some heuristics that trigger if the default value :unspecified is given.

I think option 2 is feasible. I will give it a go. Also having option 1, so that the user can override it, is probably favorable.

[Datseris]

Great. I think that's the only thing remaining from the code side. I will have a final review of the docs though. and then we can merge!

[oameye]

I'll ask @reykboerner if he is planning to review through another channel.

The implementation of coupledSDEs is already implemented as an extension 😁

[oameye]

Ha but StochasticDiffEq is still an main dependency in the project.toml. we should change this (on phone atm)

[ChrisRackauckas]

A split is planned, but it will take time.

[reykboerner]

I was waiting on @reykboerner for final commments as I was/am too busy for an additional thorough look. Do you know if he plans to submit a review soon? If not I can marge.

I will finish and submit my review today.

[reykboerner]

Hey @oameye, thanks a lot for all the work so far and sorry for the delay in reviewing!

There are a few minor changes I would implement, plus a question about the robustness of testing the noise function for time-independence, state-independence etc.

Let me know if I can help with making any of these changes!

I also have only skimmed the docs text and the tests so far.

[Datseris]

I assume this is good to go @reykboerner and @oameye ?

[oameye]

Not yet. Let me take a quick look later today

[reykboerner]

I assume this is good to go @reykboerner and @oameye ?

I need to finish the remaining edits tomorrow morning, we'll let you know!

[reykboerner]

From my side this is ready now @oameye @Datseris!

[oameye]

@reykboerner Thanks a lot Reyk for carrying it to finish line!

@Datseris We are also good to go for me! :D

[Datseris]

mucho successo !

[oameye]

So we merge?

[Datseris]

Fantastic job @oameye @reykboerner . I am merging and tagging now.

@oameye can you paste me again please the example you used to show that CoupledODEs works out of the box with Attractors.jl? I will put it into the main tutorial of DynamicalSystems.jl.

[Datseris]

I assume it was this one:

using Attractors, StaticArrays
using DynamicalSystemsBase

using DynamicalSystemsBase: idfunc

function fitzhugh_nagumo(u, p, t)
    x, y = u
    ϵ, β, α, γ, κ, I = p

    dx = (-α * x^3 + γ * x - κ * y + I) / ϵ
    dy = -β * y + x

    return SA[dx, dy]
end

sde = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), [1.,3.,1.,1.,1.,0.], 0.01)

function featurizer(X, t)
    return X[end]
end

mapper = AttractorsViaFeaturizing(sde, featurizer; Ttr = 200, T = 1)

xg = yg = range(-1, 1; length = 101)

region = HRectangle([-1, 1], [1, 1])
sampler, _ = statespace_sampler(region)

fs = basins_fractions(mapper, sampler)

but needs updating (idfunc was removed?)

[oameye]

using Attractors, StaticArrays
using DynamicalSystemsBase, StochasticDiffEq # load extention for CoupledSDEs

function fitzhugh_nagumo(u, p, t)
    x, y = u
    ϵ, β, α, γ, κ, I = p

    dx = (-α * x^3 + γ * x - κ * y + I) / ϵ
    dy = -β * y + x

    return SA[dx, dy]
end

sde = CoupledSDEs(fitzhugh_nagumo, zeros(2), [1.,3.,1.,1.,1.,0.], noise_strength = 0.05)

function featurizer(X, t)
    return X[end]
end

mapper = AttractorsViaFeaturizing(sde, featurizer; Ttr = 200, T = 1)

xg = yg = range(-1, 1; length = 101)

region = HRectangle([-1, 1], [1, 1])
sampler, _ = statespace_sampler(region)

fs = basins_fractions(mapper, sampler)

[oameye]

We probably should write a disclaimer in the example that the attractors are found faster when using CoupledODEs.

[reykboerner]

And we could also note that basins of attraction aren't well defined anymore theoretically in the presence of noise (essentially, noise makes any multistable system monostable); this example should simply demonstrate that existing functions work with CoupledSDEs as one might expect.