revised docs after v0.3 release

Author: reykboerner

CHANGELOG.md

@@ -1,5 +1,14 @@
 # Changelog for `CriticalTransitions.jl`
 
+## v0.3.0
+Major overhaul introducing `CoupledSDEs`
+
+This release replaces the `StochSystem` struct with the new `CoupledSDEs` struct to define a stochastic dynamical system. To see how the new version works, check out the [documentation](https://juliadynamics.github.io/CriticalTransitions.jl/dev/).
+
+The update is a breaking change for almost all functions, because the interface is now built around `CoupledSDEs`. The benefit is that the package now integrates much more seamlessly with DynamicalSystems.jl and DifferentialEquations.jl. To use the package with the old `StochSystem` struct, choose version v0.2.1 or lower.
+
+Full changelog [here](https://github.com/JuliaDynamics/CriticalTransitions.jl/compare/v0.2.1...v0.3.0)
+
 ## v0.2.1
 Freeze before major overhaul `v0.3.0`
 

docs/pages.jl

@@ -4,7 +4,7 @@ pages = [
     "Quickstart" => "quickstart.md",
     "Tutorial" => "tutorial.md",
     "Manual" => Any[
-        "Define a CoupledSDE" => "man/CoupledSDEs.md",
+        "Define a CoupledSDEs system" => "man/CoupledSDEs.md",
         "Stability analysis" => "man/systemanalysis.md",
         "Simulating the system" => "man/simulation.md",
         "Sampling transitions" => "man/sampling.md",

docs/src/figs/CTjl_structure_v0.3.jpg

@@ -4,7 +4,7 @@
 
 Building on [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSystems.jl/stable/) and [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/), this package aims to provide a toolbox for dynamical systems under time-dependent forcing, with a focus on tipping phenomena and metastability.
 
-![CT.jl infographic](./figs/CTjl_structure.png)
+![CT.jl infographic](./figs/CTjl_structure_v0.3_small.jpeg)
 
 !!! info "Current features"
     * **Stability analysis**: Fixed points, linear stability, basins of attraction, edge tracking
@@ -16,6 +16,7 @@ Building on [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSyste
     * **Rare event simulation**: importance sampling, AMS
     * **Quasipotentials**: Ordered line integral method (OLIM)
     * **Rate-induced tipping** tools
+    * **Symbolic** differentiation of action functionals
     * ...?
 
 

docs/src/figs/CTjl_structure_v0.3_small.jpeg

@@ -1,40 +1,39 @@
-# Define a CoupledSDE
+# 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 + g(\vec x(t);  \ p) \text{d}\mathcal{W} \ ,
+\text{d}\vec x = f(\vec x(t); \ p)  \text{d}t + \sigma (\vec x(t);  \ p) \, \text{d}\mathcal{W} \ ,
 ```
-where $\text{d}\mathcal{W}=\Gamma \cdot \text{d}\mathcal{N}$, $\vec x \in \mathbb{R}^\text{dim}$ and $\mathcal N$ denotes a stochastic process. The (positive definite) noise covariance matrix is $\Sigma = \Gamma \Gamma^\top \in \mathbb R^{N\times N}$.
+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 is assumed to be of similar form as is accepted 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 $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$.
 
-The function $g$ represent the stochastics dynamics of the system and should be the of the same type (iip or oop) as $f$.
-
- The keyword `noise` defines the system [noise process](#noise-process). In combination with `g` one can define different type of stochastic systems. Examples of different type of stochastics systems can be found on the [StochasticDiffEq.jl tutorial page](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/sde_example/). A quick overview of the different types of stochastic systems can be found [here](#Type-of-stochastic-system).
+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}$ are not Stochasitic Differential Equations but are a different class of differential equations 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.
-
-
+    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
 ```
-## Type of stochastic system
-Let us make some examples of the different types of stochastic systems that can be defined.
+
+## 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
 import Random # hide
 Random.seed!(10) # hide
-f!(du, u, p, t) = du .= 1.01u
-σ = 0.25
+f!(du, u, p, t) = du .= 1.01u # deterministic part
+σ = 0.25 # noise strength
 ```
 ### Additive noise
-When `g` is independent of the state variables `u`, the noise is called additive.
+When `g \, \text{d}\mathcal{W}` is independent of the state variables `u`, the noise is called additive.
 
 #### Diagonal noise
-A system of diagional 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`.
+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`.
 ```@example type
 t0 = 0.0; W0 = zeros(2);
 W = WienerProcess(t0, W0, 0.0)
@@ -100,9 +99,9 @@ plot(sol)
 ```
 
 !!! warning
-    Non-diagonal problem need specific type of solvers. See the [SciML recommendations](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/#sde_solve).
+    Non-diagonal problems need specific type of solvers. See the [SciML recommendations](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/#sde_solve).
 
-### Corelated noise
+### Correlated noise
 ```@example type
 ρ = 0.3
 Σ = [1 ρ; ρ 1]
@@ -113,11 +112,53 @@ sol = simulate(sde, 1.0, dt=0.01, alg=SOSRA())
 plot(sol)
 ```
 
-## Noise process
-We provide the noise processes $\text{d}\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:
+## 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
+```
+
+## Interface to `DynamicalSystems.jl`
+
+!!! 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.
+    
+
+    - 
+    type of `DynamicalSystems.jl` using the function [`CoupledODEs`](@ref). Vice vera,
+    a `CoupledODEs` system can be converted into a `CoupledSDEs` via `CoupledSDEs(ds::CoupledODEs, g)` with g the noise function.
+
+#### Converting between `CoupledSDEs` and `CoupledODEs`
+
+- `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`
+
+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:
+
+```@example type
+using CriticalTransitions
+using DynamicalSystems: lyapunovspectrum
+
+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
+
+sys = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), [1.,3.,1.,1.,1.,0.], 0.1)
+ls = lyapunovspectrum(CoupledODEs(sys), 10000)
+```
+
+```@docs
+CoupledODEs
 ```

docs/src/index.md

@@ -1,39 +1,4 @@
-# Convenience functions and types
-
-## Interface to `DynamicalSystems.jl`
-
-!!! tip "Using the functionality of DynamicalSystems.jl"
-    A [`CoupledSDEs`](@ref) can easily be turned into a 
-    [`CoupledODEs`](https://juliadynamics.github.io/DynamicalSystems.jl/dev/tutorial/#DynamicalSystemsBase.CoupledODEs)
-    instance of `DynamicalSystems.jl` using the function [`CoupledODEs`](@ref). Vice vera,
-    a `CoupledODEs` system can be converted into a `CoupledSDEs` via `CoupledSDEs(ds::CoupledODEs, g)` with g the noise function.
-
-This way, many of the methods in `DynamicalSystems.jl` can be used directly, even if we have
-not written an analogue method that takes `CoupledSDEs` as input. For example, the
-Lyapunov spectrum of a `CoupledSDEs`, here exemplified by the FitzHugh-Nagumo model, can be
-computed by typing:
-
-```julia
-using CriticalTransitions
-using DynamicalSystems: lyapunovspectrum
-
-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
-
-sys = CoupledSDEs(fitzhugh_nagumo, id_func, zeros(2), σ, [1.,3.,1.,1.,1.,0.])
-ls = lyapunovspectrum(CoupledODEs(sys), 10000)
-```
-
-```@docs
-CoupledODEs(sys::CoupledSDEs; diffeq, t0=0.0)
-```
+# Utility functions
 
 ## `CoupledSDEs` utility functions
 
@@ -56,4 +21,4 @@ make_h5(text::String, relpath::String="")
 
 ```@docs
 intervals_to_box(bmin::Vector, bmax::Vector)
-```
+```

docs/src/man/CoupledSDEs.md

@@ -14,11 +14,11 @@ As this module is not published yet, there are two ways to access it:
 ## Basic usage
 The general workflow of `CriticalTransitions` essentially follows two steps:
 
-1. Define your system (see [Define a CoupledSDE](@ref))
+1. Define your system (see [Define a CoupledSDEs system](@ref))
 2. Investigate the system by calling methods (see [Methods](@ref))
 
-!!! info "Extension to RateSystem and TippingSystem"
-    We are currently working on extending the types of dynamical systems that can be studied with CriticalTransitions.jl. Particularly, we are planning to introduce the overarching structure `TippingSystem`, which has two subtypes: `CoupledSDEs` (as it already exists) and `RateSystem`, a new dynamical system type in which the system parameters may evolve in time.
+!!! info "New system type: RateSystem"
+    We are planning to introduce the the struct `RateSystem` along `CoupledSDEs`. In a `RateSystem`, the time dependence of parameters can conveniently be specified, laying the foundation for a toolbox to study rate-induced tipping, or R-tipping.
 
 ## Methods
 

docs/src/man/utils.md

@@ -3,7 +3,7 @@
 To give you an idea of how our package works, this tutorial provides some example code with explanations.
 
 ## Example: FitzHugh-Nagumo model
-Consider the FitzHugh-Nagumo model,
+Let's consider a simple 2-dimensional dynamical system - the *FitzHugh-Nagumo* model:
 
 ```math
 \begin{aligned}
@@ -19,7 +19,7 @@ Let's investigate this system under stochastic forcing.
 ### System definition
 First, we need to translate the system equations above into Julia code.
 
-This works by defining a function `f(u,p,t)` which takes as input a vector `u` of state variables (``u``,``v``), a vector `p` of parameters, and time `t`. The function must return an array of flow increments ($\text{d}u$, $\text{d}v$). For performance reasons, it is advisable to return a StaticArray `SA[du, dv]` rather than just a Vector `[du, dv]`. This is why we need the `StaticArrays` package.
+This works exactly as in [DynamicalSystems.jl](https://juliadynamics.github.io/DynamicalSystemsBase.jl/dev/) by defining a function `f(u,p,t)` which takes as input a vector `u` of state variables (``u``,``v``), a vector `p` of parameters, and time `t`. The function must return an array of flow increments ($\text{d}u$, $\text{d}v$). For performance reasons, it is advisable to return a StaticArray `SA[du, dv]` rather than just a Vector `[du, dv]`.
 
 ```@example MAIN
 using CriticalTransitions
@@ -37,14 +37,12 @@ function fitzhugh_nagumo(u,p,t)
 end
 ```
 
-Note that the system parameters `ϵ, β, α, γ, κ, I = p[1]` are unpacked as the first component of `p`. This is necessary because in CriticalTransitions.jl one can also define a separate set of parameters for the stochastic component of the system, which would then make up the second component `p[2]` ( see [Define a CoupledSDE](@ref)).
-
 !!! tip "In-place vs. out-of-place"
     The function `fitzhugh_nagumo(u,p,t)` is defined *out-of-place*. It is also possible to define the system *in-place* as `fitzhugh_nagumo!(du,u,p,t)`. For more info, see [here](https://diffeq.sciml.ai/stable/types/ode_types/).
 
 ### CoupledSDE
 
-Next, we turn the `fitzhugh_nagumo` system into a stochastic dynamical system. Suppose we would like to force both state variables ``u`` and ``v`` with additive, uncorrelated Gaussian noise of intensity ``\sigma``. This is the default case. We simply write
+Next, we construct a stochastic system with the `fitzhugh_nagumo` equation as the deterministic part. Suppose we would like to force both state variables ``u`` and ``v`` with additive, uncorrelated Gaussian noise of intensity ``\sigma``. This is the default case. We simply write
 
 ```@example MAIN
 p = [1., 3., 1., 1., 1., 0.] # Parameters (ϵ, β, α, γ, κ, I)
@@ -53,15 +51,15 @@ p = [1., 3., 1., 1., 1., 0.] # Parameters (ϵ, β, α, γ, κ, I)
 # CoupledSDE
 sys = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ)
 ```
-Here we have chosen `zeros(2)` as the initial state of the system. The length of this vector must correspond to the system's dimensionality, but for now the state is just a placeholder that aligns our syntax with that of DifferentialEquations.jl and DynamicalSystems.jl.
+Here the first field `fitzhugh_nagumo` specifies the deterministic dynamics `f` (see [Define a CoupledSDEs system](@ref)), and the second field `idfunc` specifies the noise function `g`. The `idfunc` identity function is predefined for convenience. We have chosen `zeros(2)` as the initial state of the system, which is the third field. The length of this vector must match the system's dimensionality. In the fourth field, we specify the parameter vector, which includes the parameters of `f` followed by the parameters of `g` (in this case, there are no parameters for `g`). Lastly, `σ` sets the noise strength. Since we have not specified a noise process, the default case of an uncorrelated Wiener process is used.
 
 !!! note "Multiplicative and/or correlated noise"
-    Of course, it is also possible to define more complicated noise processes than simple additive white noise. This is done by specifying a custom *noise function* and *covariance matrix* in the `CoupledSDEs` definition. For more info, see [Define a CoupledSDE](@ref).
+    Of course, it is also possible to define more complicated noise processes than simple additive white noise. This is done by specifying a custom *noise function* and *covariance matrix* in the `CoupledSDEs` definition. For more info, see [Define a CoupledSDEs system](@ref).
 
-That's it! Now we can throw the toolbox of `CriticalTransitions` at our stochastic FitzHugh-Nagumo system `sys`.
+That's it! Now we can apply the toolbox of `CriticalTransitions` to our stochastic FitzHugh-Nagumo system `sys`.
 
 ### Find stable equilibria
-For the parameters chosen above, the FitzHugh-Nagumo system is bistable. Let's compute the fixed points using the [`fixedpoints`](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/chaostools/stable/periodicity/#ChaosTools.fixedpoints) function from ChaosTools.jl. As this function is from the `DynamicalSystems` ecosystem, it takes a system of type `CoupledODEs` as input. We can simply convert the CoupledSDEs `sys` via the [`CoupledODEs`](@ref) function:
+For the parameters chosen above, the FitzHugh-Nagumo system is bistable. Let's compute the fixed points using the [`fixedpoints`](https://juliadynamics.github.io/CriticalTransitions.jl/dev/man/systemanalysis/#ChaosTools.fixedpoints) function. This function is borrowed from ChaosTools.jl and is loaded as an extension when we write `using ChaosTools`.
 
 ```@example MAIN
 using ChaosTools
@@ -74,21 +72,21 @@ fp1, fp2 = eqs[stab]
 ```
 
 ### Stochastic simulation
-Using the `simulate` function, we now run a simulation of our system starting out from the fixed point `fp1`:
+Using the [`simulate`](@ref) function, we now run a simulation of our system for `1e3` time units starting out from the fixed point `fp1`:
 
 ```@example MAIN
-sim = simulate(sys, 1e3, fp1, saveat=0.1)
+sim = simulate(sys, 1e3, fp1; saveat=0.1)
 ```
 
-In the keyword arguments, we have specified the time step `dt` and total duration `tmax` of the numerical time integration.
+In the keyword arguments, we have specified at which interval the solution is saved. Further keyword arguments can be used to change the solver (the default is `SOSRA()` for stochastic integration) and other settings.
 
-The simulated trajectory is stored in `sim` as a matrix with 2 rows corresponding to the state variables ``u``, ``v``, and 10,000 columns corresponding to the time steps.
+The simulated trajectory is stored in `sim` in the usual output format of the [`solve`](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/#CommonSolve.solve-Tuple%7BSciMLBase.AbstractDEProblem,%20Vararg%7BAny%7D%7D) method of DifferentialEquations.jl, including the solution `sim.u` and the vector of time points `sim.t`. The solution can also be accessed as a matrix `sim[i, t]`, where `i` is the `i`-th component of `u` and `t` the time index.
 
 Let's plot the result. Did the trajectory transition to the other attractor?
 
 ```@example MAIN
 using Plots
-plt = plot(sim[1, :], sim[2, :], xlabel="u", ylabel="v", legend=false)
+plt = plot(sim[1, :], sim[2, :]; xlabel="u", ylabel="v", legend=false)
 scatter!([fp1[1], fp2[1]], [fp1[2], fp2[2]], color=:red, markersize=4)
 xlims!(-1.2, 1.2)
 ylims!(-0.6, 0.6)