add coupledSDEs type

.gitignore

@@ -6,4 +6,5 @@ deps/build.log
 Manifest.toml
 *style.jl
 *.scss
-*.css
+*.css
+.vscode

Project.toml

@@ -23,5 +23,12 @@ Roots = "1, 2"
 SciMLBase = "1.19.5, 2"
 StateSpaceSets = "1"
 Statistics = "1"
+StochasticDiffEq = "6.66.0"
 SymbolicIndexingInterface = "0.3.4"
 julia = "1.9"
+
+[weakdeps]
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+
+[extensions]
+StochasticSystemsBase = "StochasticDiffEq"

docs/Project.toml

@@ -5,3 +5,4 @@ DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
 DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StateSpaceSets = "40b095a5-5852-4c12-98c7-d43bf788e795"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"

docs/make.jl

@@ -1,9 +1,12 @@
 cd(@__DIR__)
 
 using DynamicalSystemsBase
+using StochasticDiffEq # for extention
+StochasticSystemsBase = Base.get_extension(DynamicalSystemsBase, :StochasticSystemsBase)
 
 pages = [
     "index.md",
+    "CoupledSDEs.md",
 ]
 using DynamicalSystemsBase.SciMLBase
 
@@ -15,5 +18,5 @@ Downloads.download(
 include("build_docs_with_style.jl")
 
 build_docs_with_style(pages,
-    DynamicalSystemsBase, SciMLBase, StateSpaceSets;
+    DynamicalSystemsBase, SciMLBase, StateSpaceSets, StochasticSystemsBase;
 )

docs/src/CoupledSDEs.md

@@ -0,0 +1,59 @@
+# Define a CoupledSDEs system
+
+A `CoupledSDEs` defines a stochastic dynamical system of the form
+
+```math
+\text{d}\vec x = f(\vec x(t); \ p)  \text{d}t + \sigma g(\vec x(t);  \ p) \, \text{d}\mathcal{W} \ ,
+```
+where $\vec x \in \mathbb{R}^\text{D}$, $\sigma > 0$ is the noise strength, $\text{d}\mathcal{W}=\Gamma \cdot \text{d}\mathcal{N}$, and $\mathcal N$ denotes a stochastic process. The (positive definite) noise covariance matrix is $\Sigma = \Gamma \Gamma^\top \in \mathbb R^{N\times N}$.
+
+The function $f$ is the deterministic part of the system and follows the syntax of a `ContinuousTimeDynamicalSystem` in [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSystems.jl/latest/tutorial/), i.e., `f(u, p, t)` for out-of-place (oop) and `f!(du, u, p, t)` for in-place (iip). The function $g$ allows to specify the stochastic dynamics of the system along with the [noise process](#noise-process) $\mathcal{W}$. It should be of the same type (iip or oop) as $f$.
+
+By combining $\sigma$, $g$ and $\mathcal{W}$, you can define different type of stochastic systems. Examples of different types of stochastic systems can be found on the [StochasticDiffEq.jl tutorial page](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/sde_example/). A quick overview of common types of stochastic systems can be found [below](#Defining-stochastic-dynamics).
+
+!!! info
+    Note that nonlinear mixings of the Noise Process $\mathcal{W}$ fall into the class of random ordinary differential equations (RODEs) and is not supported in DynamicalSystems.jl at the moment. See [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/rode_example/) for an implementation.
+
+```@docs
+CoupledSDEs
+```
+
+## Converting between `CoupledSDEs` and `CoupledODEs`
+
+!!! tip "Analyzing deterministic dynamics with DynamicalSystems.jl"
+    The deterministic part of a [`CoupledSDEs`](@ref) system can easily be extracted as a 
+    [`CoupledODEs`](https://juliadynamics.github.io/DynamicalSystems.jl/dev/tutorial/#DynamicalSystemsBase.CoupledODEs), a common subtype of a `ContinuousTimeDynamicalSystem` in DynamicalSystems.jl.
+
+- `CoupledODEs(sde::CoupledSDEs)` extracts the deterministic part of `sde` as a `CoupledODEs`
+- `CoupledSDEs(ode::CoupledODEs, g)`, with `g` the noise function, turns `ode` into a `CoupledSDEs`
+
+```@docs
+CoupledODEs(::CoupledSDEs; kwargs...)
+```
+For example, the
+Lyapunov spectrum of a `CoupledSDEs` in the absence of noise, here exemplified by the
+FitzHugh-Nagumo model, can be computed by typing:
+
+```julia
+using DynamicalSystemsBase, StaticArrays
+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)
+tr = trajectory(sde, 10)
+```
+
+## Available noise processes
+We provide the noise processes $\mathcal{W}$ that can be used in the stochastic simulations through the [DiffEqNoiseProcess.jl](https://docs.sciml.ai/DiffEqNoiseProcess/stable) package. A complete list of the available processes can be found [here](https://docs.sciml.ai/DiffEqNoiseProcess/stable/noise_processes/).
+```
+
+

ext/StochasticSystemsBase.jl

@@ -0,0 +1,17 @@
+module StochasticSystemsBase
+
+using DynamicalSystemsBase: DynamicalSystemsBase, SciMLBase, correct_state, CoupledODEs,
+    CoupledSDEs, StateSpaceSets, isinplace, _delete, set_parameter!,
+    set_state!, dynamic_rule, isdeterministic, current_state,
+    DynamicalSystemsBase, _set_parameter!, u_modified!,
+    additional_details, referrenced_sciml_prob, DEFAULT_DIFFEQ, SVector
+using SciMLBase: SDEProblem, AbstractSDEIntegrator, __init, SDEFunction, step!
+using StochasticDiffEq: SOSRI
+using LinearAlgebra
+
+include("src/CoupledSDEs.jl")
+include("src/classification.jl")
+
+export CoupledSDEs, CoupledODEs
+
+end

ext/src/CoupledSDEs.jl

@@ -0,0 +1,206 @@
+###########################################################################################
+# DiffEq options
+###########################################################################################
+
+const DEFAULT_SDE_SOLVER = SOSRI() # default sciml solver
+const DEFAULT_STOCH_DIFFEQ_KWARGS = (abstol = 1e-2, reltol = 1e-2) # default sciml tol
+const DEFAULT_STOCH_DIFFEQ = (alg = DEFAULT_SDE_SOLVER, DEFAULT_STOCH_DIFFEQ_KWARGS...)
+
+# Function from user `@xlxs4`, see
+# https://github.com/JuliaDynamics/jl/pull/153
+function _decompose_into_sde_solver_and_remaining(diffeq)
+    if haskey(diffeq, :alg)
+        return (diffeq[:alg], _delete(diffeq, :alg))
+    else
+        return (DEFAULT_SDE_SOLVER, diffeq)
+    end
+end
+
+###########################################################################################
+# Constructor functions
+###########################################################################################
+
+function DynamicalSystemsBase.CoupledSDEs(
+        f,
+        g,
+        u0,
+        p = SciMLBase.NullParameters(),
+        noise_strength = 1.0;
+        t0 = 0.0,
+        diffeq = DEFAULT_STOCH_DIFFEQ,
+        noise_rate_prototype = nothing,
+        noise = nothing,
+        seed = UInt64(0)
+)
+    IIP = isinplace(f, 4) # from SciMLBase
+    @assert IIP==isinplace(g, 4) "f and g must both be in-place or out-of-place"
+
+    s = correct_state(Val{IIP}(), u0)
+    T = eltype(s)
+    prob = SDEProblem{IIP}(
+        f,
+        g,
+        s,
+        (T(t0), T(Inf)),
+        p;
+        noise_rate_prototype = noise_rate_prototype,
+        noise = noise,
+        seed = seed
+    )
+    return CoupledSDEs(prob, diffeq, noise_strength)
+end
+
+function DynamicalSystemsBase.CoupledSDEs(
+        prob::SDEProblem, diffeq = DEFAULT_STOCH_DIFFEQ, noise_strength = 1.0; special_kwargs...
+)
+    if haskey(special_kwargs, :diffeq)
+        throw(
+            ArgumentError(
+            "`diffeq` is given as positional argument when an ODEProblem is provided."
+        ),
+        )
+    end
+    IIP = isinplace(prob) # from SciMLBase
+    D = length(prob.u0)
+    P = typeof(prob.p)
+    if prob.tspan === (nothing, nothing)
+        # If the problem was made via MTK, it is possible to not have a default timespan.
+        U = eltype(prob.u0)
+        prob = SciMLBase.remake(prob; tspan = (U(0), U(Inf)))
+    end
+    sde_function = SDEFunction(prob.f, add_noise_strength(noise_strength, prob.g, IIP))
+    prob = SciMLBase.remake(prob; f = sde_function)
+    noise_type = find_noise_type(prob, IIP)
+
+    solver, remaining = _decompose_into_sde_solver_and_remaining(diffeq)
+    integ = __init(
+        prob,
+        solver;
+        remaining...,
+        # Integrators are used exclusively iteratively. There is no reason to save anything.
+        save_start = false,
+        save_end = false,
+        save_everystep = false,
+        # DynamicalSystems.jl operates on integrators and `step!` exclusively,
+        # so there is no reason to limit the maximum time evolution
+        maxiters = Inf
+    )
+    return CoupledSDEs{IIP, D, typeof(integ), P, true}(
+        integ, deepcopy(prob.p), noise_strength, diffeq, noise_type
+    )
+end
+
+"""
+    CoupledSDEs(ds::CoupledODEs, g, p [, σ]; kwargs...)
+
+Converts a [`CoupledODEs`
+](https://juliadynamics.github.io/DynamicalSystems.jl/stable/tutorial/#CoupledODEs)
+system into a [`CoupledSDEs`](@ref).
+"""
+function DynamicalSystemsBase.CoupledSDEs(
+        ds::CoupledODEs,
+        g,
+        p, # the parameter is likely changed as the diffusion function g is added.
+        noise_strength = 1.0;
+        diffeq = DEFAULT_STOCH_DIFFEQ,
+        noise_rate_prototype = nothing,
+        noise = nothing,
+        seed = UInt64(0)
+)
+    return CoupledSDEs(
+        dynamic_rule(ds),
+        g,
+        current_state(ds),
+        p,
+        noise_strength;
+        diffeq = diffeq,
+        noise_rate_prototype = noise_rate_prototype,
+        noise = noise,
+        seed = seed
+    )
+end
+
+"""
+    CoupledODEs(ds::CoupledSDEs; kwargs...)
+
+Converts a [`CoupledSDEs`](@ref) into [`CoupledODEs`](@ref).
+"""
+function DynamicalSystemsBase.CoupledODEs(
+        sys::CoupledSDEs; diffeq = DEFAULT_DIFFEQ, t0 = 0.0)
+    return CoupledODEs(
+        dynamic_rule(sys), SVector{length(sys.integ.u)}(sys.integ.u), sys.p0;
+        diffeq = diffeq, t0 = t0
+    )
+end
+
+# Pretty print
+function DynamicalSystemsBase.additional_details(ds::CoupledSDEs)
+    solver, remaining = _decompose_into_sde_solver_and_remaining(ds.diffeq)
+    return [
+        "Noise strength" => ds.noise_strength,
+        "SDE solver" => string(nameof(typeof(solver))),
+        "SDE kwargs" => remaining,
+        "Noise type" => ds.noise_type
+    ]
+end
+
+###########################################################################################
+# API - obtaining information from the system
+###########################################################################################
+
+SciMLBase.isinplace(::CoupledSDEs{IIP}) where {IIP} = IIP
+StateSpaceSets.dimension(::CoupledSDEs{IIP, D}) where {IIP, D} = D
+DynamicalSystemsBase.current_state(ds::CoupledSDEs) = current_state(ds.integ)
+DynamicalSystemsBase.isdeterministic(ds::CoupledSDEs) = false
+
+function DynamicalSystemsBase.dynamic_rule(ds::CoupledSDEs)
+    # with remake it can happen that we have nested SDEFunctions
+    # sciml integrator deals with this internally well
+    f = ds.integ.f
+    while hasfield(typeof(f), :f)
+        f = f.f
+    end
+    return f
+end
+
+# so that `ds` is printed
+function DynamicalSystemsBase.set_state!(ds::CoupledSDEs, u::AbstractArray)
+    (set_state!(ds.integ, u); ds)
+end
+SciMLBase.step!(ds::CoupledSDEs, args...) = (step!(ds.integ, args...); ds)
+
+# TODO We have to check if for SDEIntegrators a different step interruption is possible.
+function DynamicalSystemsBase.successful_step(integ::AbstractSDEIntegrator)
+    rcode = integ.sol.retcode
+    return rcode == SciMLBase.ReturnCode.Default || rcode == SciMLBase.ReturnCode.Success
+end
+
+# This is here to ensure that `u_modified!` is called
+function DynamicalSystemsBase.set_parameter!(ds::CoupledSDEs, args...)
+    _set_parameter!(ds, args...)
+    u_modified!(ds.integ, true)
+    return nothing
+end
+
+DynamicalSystemsBase.referrenced_sciml_prob(ds::CoupledSDEs) = ds.integ.sol.prob
+
+###########################################################################################
+# Utilities
+###########################################################################################
+
+function σg(σ, g)
+    return (u, p, t) -> σ .* g(u, p, t)
+end
+
+function σg!(σ, g!)
+    function (du, u, p, t)
+        g!(du, u, p, t)
+        du .*= σ
+        return nothing
+    end
+end
+
+function add_noise_strength(σ, g, IIP)
+    newg = IIP ? σg!(σ, g) : σg(σ, g)
+    return newg
+end