Using new CoupledSDEs

CHANGELOG.md

@@ -1,5 +1,27 @@
 # Changelog for `CriticalTransitions.jl`
 
+## v0.4.0
+New `CoupledSDEs` design
+
+This release upgrades CriticalTransitions.jl to be compatible with the re-design of `CoupledSDEs`, which has now been integrated in [`DynamicalSystemsBase.jl v3.11`](https://juliadynamics.github.io/DynamicalSystemsBase.jl/stable/CoupledSDEs/).
+
+Since we have updated the syntax of defining a `CoupledSDEs` system, this is a breaking change.
+
+#### Changed functions
+- `CoupledSDEs` is now constructed with different args and kwargs
+- `fw_action`, `geometric_action` and `om_action` now normalize the covariance matrix when computing the action functional
+- `noise_strength` is not a function anymore but a kwarg for `CoupledSDEs` and other functions
+- `om_action` now requires `noise_strength` as input argument
+- `relax` was renamed to `deterministic_orbit`
+- `trajectory` has been added to replace `simulate`
+
+#### Deprecations
+- `add_noise_strength`, `idfunc` and `idfunc!` are no longer needed and have been removed
+- the function `noise_strength(::CoupledSDEs)` has been removed
+- `simulate` has been removed
+
+Full changelog [here](https://github.com/JuliaDynamics/CriticalTransitions.jl/compare/v0.3.0...v0.4.0).
+
 ## v0.3.0
 Major overhaul introducing `CoupledSDEs`
 

Project.toml

@@ -2,7 +2,7 @@ name = "CriticalTransitions"
 uuid = "251e6cd3-3112-48a5-99dd-66efcfd18334"
 authors = ["Reyk Börner, Ryan Deeley, Raphael Römer, Orjan Ameye"]
 repo = "https://github.com/juliadynamics/CriticalTransitions.jl.git"
-version = "0.3.0"
+version = "0.4.0"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -29,6 +29,7 @@ StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [weakdeps]
+StateSpaceSets = "40b095a5-5852-4c12-98c7-d43bf788e795"
 Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
 
@@ -38,31 +39,32 @@ CoupledSDEsBaisin = ["ChaosTools", "Attractors"]
 
 [compat]
 Aqua = "0.8.7"
-Attractors = "1.18"
-ChaosTools = "3.1"
+Attractors = "1.19.12"
+ChaosTools = "3.2.1"
 Dates = ">=1.9.0"
 Dierckx = "0.5.3"
 DiffEqNoiseProcess = "5.22"
 DocStringExtensions = "0.9.3"
 Documenter = "^1.4.1"
-DynamicalSystemsBase = "3.7.1"
+DynamicalSystemsBase = "3.11.1"
 ExplicitImports = "1.9"
 Format = "1"
 ForwardDiff = "0.10.36"
 HDF5 = "0.17.1"
 IntervalArithmetic = "0.20"
-JET = "0.9"
+JET = "0.9.9"
 JLD2 = "0.4.46"
 Markdown = ">=1.9.0"
 ModelingToolkit = "9.25"
 Optim = "1.9.3"
-OrdinaryDiffEq = "6.82"
+OrdinaryDiffEq = "6.89"
 ProgressBars = "1.5.1"
 ProgressMeter = "1.10.0"
 Reexport = "1.2.2"
+StateSpaceSets = "2.1.2"
 StaticArrays = "1.9.3"
 Statistics = ">=1.9"
-StochasticDiffEq = "6.65.1"
+StochasticDiffEq = "6.69"
 Symbolics = "5.32"
 julia = "1.9"
 

README.md

@@ -26,22 +26,22 @@ function fitzhugh_nagumo(u, p, t)
     dx = (-α * x^3 + γ * x - κ * y + I) / ϵ
     dy = -β * y + x
 
-    return SA[dx, dy]
+    return SVector{2}([dx, dy])
 end
 
 # System parameters
-p = [1., 3., 1., 1., 1., 0.]
+params = [1., 3., 1., 1., 1., 0.]
 noise_strength = 0.02
+initial_state = zeros(2)
 
 # Define stochastic system
-sys = CoupledSDEs(fitzhugh_nagumo, id_func, zeros(2), p, noise_strength)
+sys = CoupledSDEs(fitzhugh_nagumo, initial_state, params; noise_strength)
 
-# Get stable fixed points
-fps, eigs, stab = fixedpoints(sys, [-2,-2], [2,2])
-fp1, fp2 = fps[stab]
+# Run a sample trajectory
+traj = trajectory(sys, 10.0)
 
-# Generate noise-induced transition from one fixed point to the other
-path, times, success = transition(sys, fp1, fp2)
+# Compute minimum action path using gMAM algorithm
+instanton = geometric_min_action_method(sys, initial_state, current_state(sys))
 
 # ... and more, check out the documentation!
 ```

docs/Project.toml

@@ -1,14 +1,16 @@
 [deps]
 Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
 CriticalTransitions = "251e6cd3-3112-48a5-99dd-66efcfd18334"
+DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
+DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 
 [compat]
 Documenter = "1"

docs/make.jl

@@ -4,6 +4,7 @@ using DocumenterInterLinks
 using Pkg
 
 using CriticalTransitions, ChaosTools, Attractors
+using DynamicalSystemsBase
 
 project_toml = Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))
 package_version = project_toml["version"]
@@ -45,12 +46,14 @@ makedocs(;
         CriticalTransitions.DiffEqNoiseProcess,
         Base.get_extension(CriticalTransitions, :ChaosToolsExt),
         Base.get_extension(CriticalTransitions, :CoupledSDEsBaisin),
+        DynamicalSystemsBase,
+        Base.get_extension(DynamicalSystemsBase, :StochasticSystemsBase)
     ],
     doctest=false,
     format=Documenter.HTML(; html_options...),
-    warnonly=[:missing_docs],
+    warnonly=[:missing_docs, :linkcheck, :cross_references],
     pages=pages,
     plugins=[bib, links],
 )
 
-deploydocs(; repo="github.com/JuliaDynamics/CriticalTransitions.jl.git", push_preview=false)
+deploydocs(; repo="github.com/JuliaDynamics/CriticalTransitions.jl.git", push_preview=false)

docs/pages.jl

@@ -5,7 +5,7 @@ pages = [
     "Tutorial" => "tutorial.md",
     "Manual" => Any[
         "Define a CoupledSDEs system" => "man/CoupledSDEs.md",
-        "Stability analysis" => "man/systemanalysis.md",
+        #"Stability analysis" => "man/systemanalysis.md",
         "Simulating the system" => "man/simulation.md",
         "Sampling transitions" => "man/sampling.md",
         "Large deviation theory" => "man/largedeviations.md",

docs/src/index.md

@@ -7,10 +7,9 @@ Building on [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSyste
 ![CT.jl infographic](./figs/CTjl_structure_v0.3_small.jpeg)
 
 !!! info "Current features"
-    * **Stability analysis**: Fixed points, linear stability, basins of attraction, edge tracking
-    * **Stochastic simulation**: Gaussian noise, uncorrelated and correlated, additive and multiplicative
-    * **Transition path sampling**: Ensemble sampling by direct simulation and Pathspace Langevin MCMC
-    * **Large deviation theory**: Action functionals and minimization algorithms (MAM, gMAM)
+    * **Stochastic simulation** made easy: Gaussian noise, uncorrelated and correlated, additive and multiplicative
+    * **Transition path sampling**: Parallelized ensemble rejection sampling
+    * **Large deviation theory** tools: Action functionals and minimization algorithms (MAM, gMAM)
 
 !!! ukw "Planned features"
     * **Rare event simulation**: importance sampling, AMS

docs/src/man/CoupledSDEs.md

@@ -1,125 +1,118 @@
 # 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 (\vec x(t);  \ p) \, \text{d}\mathcal{W} \ ,
+```@docs
+CoupledSDEs
 ```
-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](#Type-of-stochastic-system).
 
 !!! info
     Note that nonlinear mixings of the Noise Process $\mathcal{W}$ fall into the class of random ordinary differential equations (RODEs) which have a separate set of solvers. See [this example](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/rode_example/) of DifferentialEquations.jl.
 
-```@docs
-CoupledSDEs
-```
 
-## Defining stochastic dynamics
+## [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 CriticalTransitions, Plots
+using DynamicalSystemsBase, StochasticDiffEq, DiffEqNoiseProcess
+using CairoMakie
 import Random # hide
 Random.seed!(10) # hide
 f!(du, u, p, t) = du .= 1.01u # deterministic part
-σ = 0.25 # noise strength
+
+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 \, \text{d}\mathcal{W}` is independent of the state variables `u`, the noise is called additive.
 
-#### Diagonal noise
-A system of diagonal noise is the most common type of noise. It is defined by a vector of random numbers `dW` whose size matches the output of `g` where the noise is applied element-wise, i.e. `g.*dW`.
+### 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!, idfunc!, zeros(2), nothing, σ; noise=W)
+sde = CoupledSDEs(f!, zeros(2);
+    noise_process=W, covariance=[1 0; 0 1], noise_strength=1.0
+    );
 ```
-or equivalently
+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
-sde = CoupledSDEs(f!, idfunc!, zeros(2), nothing, σ)
+tr = trajectory(sde, 1.0)
+plot_trajectory(tr...)
 ```
-where we used the helper function
-```@docs
-idfunc!
-idfunc
+
+#### 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)
 ```
-The vector `dW` is by default zero mean white gaussian noise $\mathcal{N}(0, \text{d}t)$ where the variance is the timestep $\text{d}t$ unit variance (Wiener Process).
+Alternatively, we can parametrise the covariance matrix by defining the diffusion function $g$ ourselves:
 ```@example type
-sol = simulate(sde, 1.0, dt=0.01, alg=SOSRA())
-plot(sol)
+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
-Scalar noise is where a single random variable is applied to all dependent variables. To do this, one has to give the noise process to the `noise` keyword of the `CoupledSDEs` constructor. A common example is the Wiener process starting at `W0=0.0` at time `t0=0.0`.
-
+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!, idfunc!, rand(2)./10, nothing, σ; noise=noise)
-sol = simulate(sde, 1.0, dt=0.01, alg=SOSRA())
-plot(sol)
+sde = CoupledSDEs(f!, rand(2)/10; noise_process=noise)
+
+tr = trajectory(sde, 1.0)
+plot_trajectory(tr...)
 ```
-### Multiplicative noise
-Multiplicative Noise is when $g_i(t, u)=a_i u$
+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[1] = σ*u[1]
-    du[2] = σ*u[2]
+function g!(du, u, p, t)
+    du .= u ./ (1+t)
     return nothing
 end
-sde = CoupledSDEs(f!, g, rand(2)./10)
-sol = simulate(sde, 1.0, dt=0.01, alg=SOSRI())
-plot(sol)
+sde = CoupledSDEs(f!, rand(2)./10; g=g!)
 ```
 
 #### Non-diagonal noise
-Non-diagonal noise allows for the terms to linearly mixed via g being a matrix. Suppose we have two Wiener processes and two dependent random variables such that the output of `g` is a 2x2 matrix. Therefore, we have
+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 can use the `noise_rate_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 utilise the structure of commutative noise.
+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
-function g(du, u, p, t)
+σ = 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
-sde = CoupledSDEs(f!, g, rand(2)./10, noise_rate_prototype = zeros(2, 2))
-sol = simulate(sde, 1.0, dt=0.01, alg=RKMilCommute())
-plot(sol)
+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 type of solvers. See the [SciML recommendations](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/#sde_solve).
-
-### Correlated noise
-```@example type
-ρ = 0.3
-Σ = [1 ρ; ρ 1]
-t0 = 0.0; W0 = zeros(2); Z0 = zeros(2);
-W = CorrelatedWienerProcess(Σ, t0, W0, Z0)
-sde = CoupledSDEs(f!, idfunc!, zeros(2), nothing, σ; noise=W)
-sol = simulate(sde, 1.0, dt=0.01, alg=SOSRA())
-plot(sol)
-```
-
-## 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/). We list some of the most common ones below:
-```@docs
-WienerProcess
-SimpleWienerProcess
-OrnsteinUhlenbeckProcess
-CorrelatedWienerProcess
-```
+    Non-diagonal problems need specific solvers. See the [SciML recommendations](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/#sde_solve).
 
 ## Interface to `DynamicalSystems.jl`
 #### Converting between `CoupledSDEs` and `CoupledODEs`
@@ -128,8 +121,8 @@ CorrelatedWienerProcess
     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`
+- `CoupledODEs(sde::CoupledSDEs)` extracts the deterministic part of `sde` as a `CoupledODEs`.
+- `CoupledSDEs(ode::CoupledODEs; kwargs)` turns `ode` into a `CoupledSDEs`.
 
 ```@docs
 CoupledODEs
@@ -152,6 +145,8 @@ function fitzhugh_nagumo(u, p, t)
     return SA[dx, dy]
 end
 
-sys = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), [1.,3.,1.,1.,1.,0.], 0.1)
+p = [1.,3.,1.,1.,1.,0.]
+
+sys = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=0.1)
 ls = lyapunovspectrum(CoupledODEs(sys), 10000)
 ```