add coupledSDEs type (#212)

* add coupledSDEs type

* add CoupledSDEs docs

* apply suggested format changes

* add idfunc

* Fix docs?

* change coupled sde example to use trajectory

* specify which CoupledODEs docstring

* change to SciML default tolerances

* classify noise

* correct spelling of invertible

* make noisetype field a namedtuple

* make  coupledSDEs an extention

* export pretty print

add extentention to docs

clean up test file

* export CoupledODEs

* add CoupledODEs <-> CoupledSDEs tests

* add extention to modules in makedocs

* define dynamics rule sde specific

* port docs into the docstring

* delete ported documentation

* fix small typos in docstring

* add CoupledSDEs examples to docs

* remove noise strength

* update docs with no noise strength

* add covariance kwarg and change g to kwarg

* update docs

* diffusion <-> covariance matrix

* add noise_strength kwarg

* estimate invertiability with a tolerance

* move to SOSRA with non-adaptive timesteps

update docs to use LambdaEM when we have covariance

remove accidently comitted JuliaFormatter

fix tests

Set default solver to EM

* Reyk's comments

* implemented review comments - part 1

* revised CoupledSDEs docstring

* remove broken docs link

* implemented review comments - part 2

* proofread and revised CoupledSDEs docs

* fixed bug that made tests fail

* Update CHANGELOG.md

* Update Project.toml

---------

Co-authored-by: Datseris <datseris.george@gmail.com>
Co-authored-by: reykboerner <reyk.boerner@reading.ac.uk>

.gitignore

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

CHANGELOG.md

@@ -1,3 +1,7 @@
+# v3.11.0
+
+Brand new dynamical system `CoupledSDEs` that represents stochastic differential equations. It also comes with a dedicated documentation page.
+
 # v3.10.0
 
 - OrdinaryDiffEq.jl dependency reduced to OrdinaryDiffEqTsit5.jl.

Project.toml

@@ -1,7 +1,7 @@
 name = "DynamicalSystemsBase"
 uuid = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 repo = "https://github.com/JuliaDynamics/DynamicalSystemsBase.jl.git"
-version = "3.10.1"
+version = "3.11.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -23,5 +23,12 @@ Roots = "1, 2"
 SciMLBase = "1.19.5, 2"
 StateSpaceSets = "2"
 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

@@ -1,7 +1,9 @@
 [deps]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 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/juliadynamics-dark.scss

@@ -70,6 +70,7 @@ $code-background: $codebg; // for inline code
 $pre-background: $codebg;  // for code blocks
 $documenter-docstring-header-background: lighten-color($body-background-color, 0.5);
 
+
 // Sidebar
 $documenter-sidebar-background: darken-color($maincolor, 1.2); //background color for sidebar
 $documenter-sidebar-color: $text; //font color for sidebar
@@ -86,22 +87,8 @@ $input-focus-border-color: $mainwhite;
 
 $documenter-docstring-shadow: 3px 3px 4px invert($shadow-color);
 
-
-
-// Admonition stuff
-$admbg: lighten-color($body-background-color, 0.5);
-$admonition-background: (
-  'default': $admbg, 'info': $admbg, 'success': $admbg, 'warning': $admbg,
-  'danger': $admbg, 'compat': $admbg
-);
-$admonition-header-background: (
-  'default': #ba3f1f, 'warning': #a88b17, 'danger': #c7524c,
-  'success': #42ac68, 'info': #28c);
-
-// Deprecations
-$admonition-body-color: null;
-$admonition-body-color: null;
-$admonition-header-color: null;
+// This line is required after https://github.com/JuliaDocs/Documenter.jl/pull/2499
+$admonition-default: $darkbg;
 
 // All secondary themes have to be nested in a theme--$(themename) class. When Documenter
 // switches themes, it applies this class to <html> and then disables the primary
@@ -157,4 +144,36 @@ html.theme--#{$themename} {
   .hljs-subst {
     color: #f8f8f2;
   }
+
+
+
+  // fix for the settings/modal panel
+  .modal-card-head {
+    background-color: darken($darkbg,5);
+  }
+  .modal-card-body {
+    background-color: darken($darkbg,5);
+  }
+  .modal-card-foot {
+    background-color: darken($darkbg,5);
+  }
+  .select {
+    > select {
+    background-color: darken($darkbg,5); 
+    }
+  } 
+
+  // fix for the search panel
+  .input {
+    background-color: darken($darkbg,5);
+  }
+  .search-result-title {
+    color: white
+  }
+
+  // match the size of the light theme (additional rule 
+  // that is dark-theme specific and messing around 
+  // with the size of elements, making the size be 0,75em
+  // only for the dark theme, no idea why)
+  font-size: unset !important;
 }

docs/juliadynamics-darkdefs.scss

@@ -24,6 +24,7 @@ $code-background: $codebg; // for inline code
 $pre-background: $codebg;  // for code blocks
 $documenter-docstring-header-background: lighten-color($body-background-color, 0.5);
 
+
 // Sidebar
 $documenter-sidebar-background: darken-color($maincolor, 1.2); //background color for sidebar
 $documenter-sidebar-color: $text; //font color for sidebar
@@ -40,22 +41,8 @@ $input-focus-border-color: $mainwhite;
 
 $documenter-docstring-shadow: 3px 3px 4px invert($shadow-color);
 
-
-
-// Admonition stuff
-$admbg: lighten-color($body-background-color, 0.5);
-$admonition-background: (
-  'default': $admbg, 'info': $admbg, 'success': $admbg, 'warning': $admbg,
-  'danger': $admbg, 'compat': $admbg
-);
-$admonition-header-background: (
-  'default': #ba3f1f, 'warning': #a88b17, 'danger': #c7524c,
-  'success': #42ac68, 'info': #28c);
-
-// Deprecations
-$admonition-body-color: null;
-$admonition-body-color: null;
-$admonition-header-color: null;
+// This line is required after https://github.com/JuliaDocs/Documenter.jl/pull/2499
+$admonition-default: $darkbg;
 
 // All secondary themes have to be nested in a theme--$(themename) class. When Documenter
 // switches themes, it applies this class to <html> and then disables the primary
@@ -111,4 +98,36 @@ html.theme--#{$themename} {
   .hljs-subst {
     color: #f8f8f2;
   }
+
+
+
+  // fix for the settings/modal panel
+  .modal-card-head {
+    background-color: darken($darkbg,5);
+  }
+  .modal-card-body {
+    background-color: darken($darkbg,5);
+  }
+  .modal-card-foot {
+    background-color: darken($darkbg,5);
+  }
+  .select {
+    > select {
+    background-color: darken($darkbg,5); 
+    }
+  } 
+
+  // fix for the search panel
+  .input {
+    background-color: darken($darkbg,5);
+  }
+  .search-result-title {
+    color: white
+  }
+
+  // match the size of the light theme (additional rule 
+  // that is dark-theme specific and messing around 
+  // with the size of elements, making the size be 0,75em
+  // only for the dark theme, no idea why)
+  font-size: unset !important;
 }

docs/make.jl

@@ -1,9 +1,13 @@
 cd(@__DIR__)
 
 using DynamicalSystemsBase
+using StochasticDiffEq, DiffEqNoiseProcess # to enable extention
+# We need this because Documenter doesn't know where to get the docstring from otherwise
+StochasticSystemsBase = Base.get_extension(DynamicalSystemsBase, :StochasticSystemsBase)
 
 pages = [
     "index.md",
+    "CoupledSDEs.md",
 ]
 using DynamicalSystemsBase.SciMLBase
 
@@ -15,5 +19,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,111 @@
+# `CoupledSDEs`
+
+```@docs
+CoupledSDEs
+```
+
+## [Examples defining stochastic dynamics](@id defining-stochastic-dynamics)
+
+Let's look at some examples of the different types of stochastic systems that can be defined.
+
+For simplicity, we choose a slow exponential growth in 2 dimensions as the deterministic dynamics `f`:
+```@example type
+using DynamicalSystemsBase, StochasticDiffEq, DiffEqNoiseProcess
+using CairoMakie
+import Random # hide
+Random.seed!(10) # hide
+f!(du, u, p, t) = du .= 1.01u # deterministic part
+
+function plot_trajectory(Y, t)
+    fig = Figure()
+    ax = Axis(fig[1,1]; xlabel = "time", ylabel = "variable")
+    for var in columns(Y)
+        lines!(ax, t, var)
+    end
+    fig
+end;
+```
+
+### Additive noise
+When $g(u, p, t)$ is independent of the state $u$, the noise is called additive; otherwise, it is multiplicative.
+We can define a simple additive noise system as follows:
+```@example type
+sde = CoupledSDEs(f!, zeros(2));
+```
+which is equivalent to
+```@example type
+t0 = 0.0; W0 = zeros(2);
+W = WienerProcess(t0, W0, 0.0)
+sde = CoupledSDEs(f!, zeros(2);
+    noise_process=W, covariance=[1 0; 0 1], noise_strength=1.0
+    );
+```
+We defined a Wiener process `W`, whose increments are vectors of normally distributed random numbers of length matching the output of `g`. The noise is applied element-wise, i.e., `g.*dW`. Since the noise processes are uncorrelated, meaning the covariance matrix is diagonal, this type of noise is referred to as diagonal.
+
+We can sample a trajectory from this system using the `trajectory` function also used for the deterministic systems:
+```@example type
+tr = trajectory(sde, 1.0)
+plot_trajectory(tr...)
+```
+
+#### Correlated noise
+In the case of correlated noise, the random numbers in a vector increment `dW` are correlated. This can be achieved by specifying the covariance matrix $\Sigma$ via the `covariance` keyword:
+```@example type
+ρ = 0.3
+Σ = [1 ρ; ρ 1]
+diffeq = (alg = LambaEM(), dt=0.1)
+sde = CoupledSDEs(f!, zeros(2); covariance=Σ, diffeq=diffeq)
+```
+Alternatively, we can parametrise the covariance matrix by defining the diffusion function $g$ ourselves:
+```@example type
+g!(du, u, p, t) = (du .= [1 p[1]; p[1] 1]; return nothing) 
+sde = CoupledSDEs(f!, zeros(2), (ρ); g=g!, noise_prototype=zeros(2, 2))
+```
+Here, we had to provide `noise_prototype` to indicate that the diffusion function `g` will output a 2x2 matrix.
+
+#### Scalar noise
+If all state variables are forced by the same single random variable, we have scalar noise.
+To define scalar noise, one has to give an one-dimensional noise process to the `noise_process` keyword of the `CoupledSDEs` constructor. 
+```@example type
+t0 = 0.0; W0 = 0.0;
+noise = WienerProcess(t0, W0, 0.0)
+sde = CoupledSDEs(f!, rand(2)/10; noise_process=noise)
+
+tr = trajectory(sde, 1.0)
+plot_trajectory(tr...)
+```
+We can see that noise applied to each variable is the same.
+
+### Multiplicative and time-dependent noise
+In the SciML ecosystem, multiplicative noise is defined through the condition $g_i(t, u)=a_i u$. However, in the literature the name is more broadly used for any situation where the noise is non-additive and depends on the state $u$, possibly also in a non-linear way. When defining a `CoupledSDEs`, we can make the noise term time- and state-dependent by specifying an explicit time- or state-dependence in the noise function `g`, just like we would define `f`. For example, we can define a system with temporally decreasing multiplicative noise as follows:
+```@example type
+function g!(du, u, p, t)
+    du .= u ./ (1+t)
+    return nothing
+end
+sde = CoupledSDEs(f!, rand(2)./10; g=g!)
+```
+
+#### Non-diagonal noise
+Non-diagonal noise allows for the terms to be linearly mixed (correlated) via `g` being a matrix. Suppose we have two Wiener processes and two state variables such that the output of `g` is a 2x2 matrix. Therefore, we have
+```math
+du_1 = f_1(u,p,t)dt + g_{11}(u,p,t)dW_1 + g_{12}(u,p,t)dW_2 \\
+du_2 = f_2(u,p,t)dt + g_{21}(u,p,t)dW_1 + g_{22}(u,p,t)dW_2
+```
+To indicate the structure that `g` should have, we must use the `noise_prototype` keyword. Let us define a special type of non-diagonal noise called commutative noise. For this we can utilize the `RKMilCommute` algorithm which is designed to utilize the structure of commutative noise.
+
+```@example type
+σ = 0.25 # noise strength
+function g!(du, u, p, t)
+  du[1,1] = σ*u[1]
+  du[2,1] = σ*u[2]
+  du[1,2] = σ*u[1]
+  du[2,2] = σ*u[2]
+    return nothing
+end
+diffeq = (alg = RKMilCommute(), reltol = 1e-3, abstol = 1e-3, dt=0.1)
+sde = CoupledSDEs(f!, rand(2)./10; g=g!, noise_prototype = zeros(2, 2), diffeq = diffeq)
+```
+
+!!! warning
+    Non-diagonal problems need specific solvers. See the [SciML recommendations](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/#sde_solve).

ext/StochasticSystemsBase.jl

@@ -0,0 +1,18 @@
+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, SMatrix, current_parameters
+using SciMLBase: SDEProblem, AbstractSDEIntegrator, __init, SDEFunction, step!
+using StochasticDiffEq: SOSRA
+using LinearAlgebra
+
+include("src/CoupledSDEs.jl")
+include("src/classification.jl")
+
+export CoupledSDEs, CoupledODEs
+
+end