From 1d610aa5daa61af5ec331dbe418c601e610e2ab5 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Thu, 10 Oct 2024 13:58:11 +0100 Subject: [PATCH 1/3] less abrupt format CI --- .github/workflows/Format.yml | 42 +++++++++++++++++++++++++++++++----- 1 file changed, 37 insertions(+), 5 deletions(-) diff --git a/.github/workflows/Format.yml b/.github/workflows/Format.yml index 76aebbb2..c42031de 100644 --- a/.github/workflows/Format.yml +++ b/.github/workflows/Format.yml @@ -1,9 +1,41 @@ -name: Format suggestions +name: format-check -on: [pull_request] +on: + push: + branches: + - 'main' + tags: '*' + pull_request: jobs: - code-style: - runs-on: ubuntu-latest + build: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1] + julia-arch: [x86] + os: [ubuntu-latest] steps: - - uses: julia-actions/julia-format@v3 \ No newline at end of file + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + + - uses: actions/checkout@v4 + - name: Install JuliaFormatter and format + # This will use the latest version by default but you can set the version like so: + # + # julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.13.0"))' + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' + julia -e 'using JuliaFormatter; format(".", verbose=true)' + - name: Format check + run: | + julia -e ' + out = Cmd(`git diff`) |> read |> String + if out == "" + exit(0) + else + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) + end' From 9151cae941ff6519e499a3d8acdd617ec560332f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Reyk=20B=C3=B6rner?= Date: Fri, 11 Oct 2024 23:30:18 +0100 Subject: [PATCH 2/3] Using new `CoupledSDEs` from DynamicalSystemsBase (#116) * load CoupledSDEs from DynamicalSystemsBase v3.11 * worked on compatibility with new CoupledSDEs * updated LDT functions, made om_action require noise_strength input * export functions from StochasticSystemsBase * worked on updating tests * CoupledSDEs test's * try out new format action * turn on all tests * format * fix downgrade compat up JET compat fix downgrade compat? fix downgrade compat? fix downgrade compat? up StochasticDiffEq compat * better format CI try another format github action turn on formatter in PR * update code example in README * fixed Large Deviations tests by normalizing covariance matrix * updated simulation functions * added simulation tests * removed 'simulate', renamed 'relax' to 'deterministic_orbit' * tried to fix docs, still not finding docstrings from DynamicalSystemsBase * add stochasticSystemsBase do docs modules * fixed typo in makedocs * exported missing functions, removed unneeded dep * worked on fixing docs errors * added dependency to docs * added another dep to docs * updated changelog to v0.4.0 * fixed errors in docstrings --------- Co-authored-by: Orjan Ameye --- CHANGELOG.md | 22 ++ Project.toml | 16 +- README.md | 16 +- docs/Project.toml | 6 +- docs/make.jl | 7 +- docs/pages.jl | 2 +- docs/src/index.md | 7 +- docs/src/man/CoupledSDEs.md | 139 +++++------ docs/src/man/largedeviations.md | 2 +- docs/src/man/sampling.md | 4 +- docs/src/man/simulation.md | 13 +- docs/src/man/systemanalysis.md | 7 - docs/src/man/utils.md | 9 +- docs/src/quickstart.md | 2 +- docs/src/tutorial.md | 10 +- ext/ChaosToolsExt.jl | 3 +- ext/basin/basinsofattraction.jl | 4 +- src/CoupledSDEs.jl | 230 ------------------ src/CoupledSDEs_utils.jl | 62 ----- src/CriticalTransitions.jl | 33 +-- src/largedeviations/action.jl | 51 ++-- src/system_utils.jl | 47 ++++ src/trajectories/equib.jl | 11 +- src/trajectories/simulation.jl | 41 +--- src/trajectories/transition.jl | 24 +- src/utils.jl | 13 + test/CoupledSDEs.jl | 183 +++----------- test/basin/basin_boundary.jl | 18 ++ test/ext/ChaosToolsExt.jl | 2 +- test/ext/CoupledSDEsBaisin.jl | 2 +- test/issue14.jl | 22 ++ test/largedeviations/MAM.jl | 4 +- test/largedeviations/action_fhn.jl | 11 +- test/largedeviations/gMAM.jl | 4 +- test/runtests.jl | 8 +- test/trajactories/simulate.jl | 18 -- test/trajectories/simulate.jl | 13 + .../transition.jl | 7 +- test/utils.jl | 20 -- 39 files changed, 377 insertions(+), 716 deletions(-) delete mode 100644 src/CoupledSDEs.jl delete mode 100644 src/CoupledSDEs_utils.jl create mode 100644 src/system_utils.jl create mode 100644 test/basin/basin_boundary.jl create mode 100644 test/issue14.jl delete mode 100644 test/trajactories/simulate.jl create mode 100644 test/trajectories/simulate.jl rename test/{trajactories => trajectories}/transition.jl (78%) diff --git a/CHANGELOG.md b/CHANGELOG.md index 92e742b2..a8442750 100644 --- a/CHANGELOG.md +++ b/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` diff --git a/Project.toml b/Project.toml index 4d21ceb8..7f744cd8 100644 --- a/Project.toml +++ b/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" diff --git a/README.md b/README.md index 119b67cb..17da4c13 100644 --- a/README.md +++ b/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! ``` diff --git a/docs/Project.toml b/docs/Project.toml index 4c4e5b37..c28f6c33 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/make.jl b/docs/make.jl index c9c227a1..9bda4a59 100644 --- a/docs/make.jl +++ b/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) \ No newline at end of file diff --git a/docs/pages.jl b/docs/pages.jl index 9d5c1bce..879e5aea 100644 --- a/docs/pages.jl +++ b/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", diff --git a/docs/src/index.md b/docs/src/index.md index 2baf6d2a..0bf28ae5 100644 --- a/docs/src/index.md +++ b/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 diff --git a/docs/src/man/CoupledSDEs.md b/docs/src/man/CoupledSDEs.md index 4bc9e786..97a1ade2 100644 --- a/docs/src/man/CoupledSDEs.md +++ b/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) ``` \ No newline at end of file diff --git a/docs/src/man/largedeviations.md b/docs/src/man/largedeviations.md index 543c69fc..99ea8d0b 100644 --- a/docs/src/man/largedeviations.md +++ b/docs/src/man/largedeviations.md @@ -25,7 +25,7 @@ action ## Minimum action paths We provide the following two methods to calculate *instantons*, or minimum action paths, -between two states of a `CoupledSDEs`. +between two states of a `CoupledSDEs` system. ### Minimum action method (MAM) Minimization of the Freidlin-Wentzell action using the L-BFGS algorithm of `Optim.jl`. diff --git a/docs/src/man/sampling.md b/docs/src/man/sampling.md index 07ba8c09..f27b04f3 100644 --- a/docs/src/man/sampling.md +++ b/docs/src/man/sampling.md @@ -8,4 +8,6 @@ transition(sys::CoupledSDEs, x_i, x_f; kwargs...) transitions(sys::CoupledSDEs, x_i, x_f, N=1; kwargs...) ``` -## ... in pathspace \ No newline at end of file +## ... in pathspace + +Coming soon... \ No newline at end of file diff --git a/docs/src/man/simulation.md b/docs/src/man/simulation.md index 43802e32..a1e871c6 100644 --- a/docs/src/man/simulation.md +++ b/docs/src/man/simulation.md @@ -1,16 +1,17 @@ # Simulating the system We provide two main functions to simulate a [`CoupledSDEs`](@ref) forward in time: -* [`relax`](@ref) for the system's deterministic evolution (in the absence of noise) based on the [`ODEProblem`](https://diffeq.sciml.ai/stable/types/ode_types/#SciMLBase.ODEProblem) of `DifferentialEquations.jl` -* [`simulate`](@ref) for stochastic dynamics based on the [`SDEProblem`](https://diffeq.sciml.ai/stable/types/sde_types/#SciMLBase.SDEProblem) of `DifferentialEquations.jl` +* [`trajectory`](@ref), which integrates the stochastic `CoupledSDEs` system forward in time +* [`deterministic_orbit`](@ref), which integrates only the deterministic part of the `CoupledSDEs` system + +## Stochastic dynamics -## Deterministic dynamics ```@docs -relax(sys::CoupledSDEs, T, init; kwargs...) +trajectory ``` -## Stochastic dynamics +## Deterministic dynamics ```@docs -simulate(sys::CoupledSDEs, T, init; kwargs...) +deterministic_orbit ``` \ No newline at end of file diff --git a/docs/src/man/systemanalysis.md b/docs/src/man/systemanalysis.md index 369c6d34..30a3f080 100644 --- a/docs/src/man/systemanalysis.md +++ b/docs/src/man/systemanalysis.md @@ -8,13 +8,6 @@ equilib(sys::CoupledSDEs, state; kwargs...) fixedpoints ``` -## Basins of attraction -```@docs -basins -basinboundary -basboundary -``` - ## Edge tracking The edge tracking algorithm is a simple numerical method to find the *edge state* or (possibly chaotic) saddle on the boundary between two basins of attraction. It is first diff --git a/docs/src/man/utils.md b/docs/src/man/utils.md index c060f175..130683de 100644 --- a/docs/src/man/utils.md +++ b/docs/src/man/utils.md @@ -3,11 +3,14 @@ ## `CoupledSDEs` utility functions ```@docs -noise_strength covariance_matrix +diffusion_matrix +dynamic_rule noise_process -CriticalTransitions.drift -CriticalTransitions.div_drift +current_state +set_state! +drift +div_drift ``` ## Saving data diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md index e5be4931..9e23255c 100644 --- a/docs/src/quickstart.md +++ b/docs/src/quickstart.md @@ -26,5 +26,5 @@ Currently the following functions are implemented to analyze a [`CoupledSDEs`](@ corresponding sample transition paths. ```@index -Pages = ["man/systemanalysis.md", "man/simulation.md", "man/sampling.md", "man/largedeviations.md"] +Pages = ["man/simulation.md", "man/sampling.md", "man/largedeviations.md"] ``` \ No newline at end of file diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 3d374784..d491c70c 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -49,9 +49,9 @@ p = [1., 3., 1., 1., 1., 0.] # Parameters (ϵ, β, α, γ, κ, I) σ = 0.2 # noise strength # CoupledSDE -sys = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ) +sys = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=σ) ``` -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. +Here the first field `fitzhugh_nagumo` specifies the deterministic dynamics `f` (see [Define a CoupledSDEs system](@ref)). We have chosen `zeros(2)` as the initial state of the system, which is the second field. The length of this vector must match the system's dimensionality. In the (optional) third field, we specify the parameter vector `p`, which includes the parameters of `f` followed by the parameters of `g` (in this case, there are no parameters for `g`). Lastly, `noise_strength` 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 CoupledSDEs system](@ref). @@ -72,10 +72,10 @@ fp1, fp2 = eqs[stab] ``` ### Stochastic simulation -Using the [`simulate`](@ref) function, we now run a simulation of our system for `1e3` time units starting out from the fixed point `fp1`: +Using the [`trajectory`](@ref) function, we now run a simulation of our system for `100` time units starting out from the fixed point `fp1`: ```@example MAIN -sim = simulate(sys, 1e3, fp1; saveat=0.1) +sim = trajectory(sys, 100, fp1) ``` 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. @@ -86,7 +86,7 @@ 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][:,1], sim[1][:,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) diff --git a/ext/ChaosToolsExt.jl b/ext/ChaosToolsExt.jl index 2f02fbbd..d2311b2d 100644 --- a/ext/ChaosToolsExt.jl +++ b/ext/ChaosToolsExt.jl @@ -4,6 +4,7 @@ using CriticalTransitions using DocStringExtensions using ForwardDiff using ChaosTools: ChaosTools, fixedpoints +using DynamicalSystemsBase: CoupledODEs, StateSpaceSet export fixedpoints @@ -23,7 +24,7 @@ Returns fixed points, their eigenvalues and stability of the system `sys` within ## Output `[fp, eigs, stable]` -* `fp`: `Dataset` of fixed points +* `fp`: `StateSpaceSet` of fixed points * `eigs`: vector of Jacobian eigenvalues of each fixed point * `stable`: vector of booleans indicating the stability of each fixed point (`true`=stable, `false`=unstable) diff --git a/ext/basin/basinsofattraction.jl b/ext/basin/basinsofattraction.jl index fd569bf2..3fcd8585 100644 --- a/ext/basin/basinsofattraction.jl +++ b/ext/basin/basinsofattraction.jl @@ -1,5 +1,3 @@ -toattractors(V::Dataset) = Dict(i => StateSpaceSet([V[i]]) for i in 1:length(V)) - """ $(TYPEDSIGNATURES) @@ -54,3 +52,5 @@ function CriticalTransitions.basins( return [X, Y, attractors, h] end + +toattractors(V::StateSpaceSet) = Dict(i => StateSpaceSet([V[i]]) for i in 1:length(V)) diff --git a/src/CoupledSDEs.jl b/src/CoupledSDEs.jl deleted file mode 100644 index fdc63067..00000000 --- a/src/CoupledSDEs.jl +++ /dev/null @@ -1,230 +0,0 @@ -using DynamicalSystemsBase.SciMLBase: __init -using DynamicalSystemsBase: - DynamicalSystemsBase, - CoupledODEs, - isinplace, - SciMLBase, - correct_state, - _set_parameter!, - current_state -using StochasticDiffEq: SDEProblem - -########################################################################################### -# DiffEq options -########################################################################################### -const DEFAULT_SOLVER = SOSRA() -const DEFAULT_DIFFEQ_KWARGS = (abstol=1e-6, reltol=1e-6) -const DEFAULT_DIFFEQ = (alg=DEFAULT_SOLVER, DEFAULT_DIFFEQ_KWARGS...) - -# Function from user `@xlxs4`, see -# https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/pull/153 -_delete(a::NamedTuple, s::Symbol) = NamedTuple{filter(≠(s), keys(a))}(a) -function _decompose_into_solver_and_remaining(diffeq) - if haskey(diffeq, :alg) - return (diffeq[:alg], _delete(diffeq, :alg)) - else - return (DEFAULT_SOLVER, diffeq) - end -end - -########################################################################################### -# Type -########################################################################################### - -# Define CoupledSDEs -""" - CoupledSDEs <: ContinuousTimeDynamicalSystem - CoupledSDEs(f, g, u0 [, p, σ]; kwargs...) - -A stochastic continuous time dynamical system defined by a set of -coupled ordinary differential equations as follows: -```math -d\\vec{u} = \\vec{f}(\\vec{u}, p, t) dt + \\vec{g}(\\vec{u}, p, t) dW_t -``` - -Optionally provide the overall noise strength `σ`, the parameter container `p` and initial time as keyword `t0`. If `σ` is provided, the diffusion function `g` is multiplied by `σ`. - -For construction instructions regarding `f, u0` see the [DynamicalSystems.jl tutorial](https://juliadynamics.github.io/DynamicalSystems.jl/latest/tutorial/#DynamicalSystemsBase.CoupledODEs). - -The stochastic part of the differential equation is defined by the function `g` and the keyword arguments `noise_rate_prototype` and `noise`. `noise` indicates the noise process applied during generation and defaults to Gaussian white noise. For details on defining various noise processes, refer to the noise process documentation page. `noise_rate_prototype` indicates the prototype type instance for the noise rates, i.e., the output of `g`. It can be any type which overloads `A_mul_B!` with itself being the middle argument. Commonly, this is a matrix or sparse matrix. If this is not given, it defaults to `nothing`, which means the problem should be interpreted as having diagonal noise. - -## DifferentialEquations.jl interfacing - -The ODEs are evolved via the solvers of DifferentialEquations.jl. -When initializing a `CoupledODEs`, you can specify the solver that will integrate -`f` in time, along with any other integration options, using the `diffeq` keyword. -For example you could use `diffeq = (abstol = 1e-9, reltol = 1e-9)`. -If you want to specify a solver, do so by using the keyword `alg`, e.g.: -`diffeq = (alg = Tsit5(), reltol = 1e-6)`. This requires you to have been first -`using OrdinaryDiffEq` to access the solvers. The default `diffeq` is: - -```julia -$(DEFAULT_DIFFEQ) -``` - -`diffeq` keywords can also include `callback` for [event handling -](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/). - -Dev note: `CoupledSDEs` is a light wrapper of `StochasticDiffEq.SDEIntegrator` from DifferentialEquations.jl. -The integrator is available as the field `integ`, and the `SDEProblem` is `integ.sol.prob`. -The convenience syntax `SDEProblem(ds::CoupledSDEs, tspan = (t0, Inf))` is available -to extract the problem. -""" -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 - -""" - StochSystem - -An alias to [`CoupledSDEs`](@ref). -This was the name these systems had in CriticalTransitions.jl. -""" -const StochSystem = CoupledSDEs - -function CoupledSDEs( - f, - g, - u0, - p=SciMLBase.NullParameters(), - noise_strength=1.0; - t0=0.0, - diffeq=DEFAULT_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 CoupledSDEs(prob::SDEProblem, diffeq=DEFAULT_DIFFEQ, noise_strength=1.0) - 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) - - solver, remaining = _decompose_into_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 - ) -end - -""" -$(TYPEDSIGNATURES) - -Converts a [`CoupledODEs`](https://juliadynamics.github.io/DynamicalSystems.jl/stable/tutorial/#DynamicalSystemsBase.CoupledODEs) -system into a [`CoupledSDEs`](@ref). -""" -function CoupledSDEs( - ds::DynamicalSystemsBase.CoupledODEs, - g, - p, # the parameter contained likely is changed as the diffusion function g is added. - noise_strength=1.0; - diffeq=DEFAULT_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 - -""" -$(TYPEDSIGNATURES) - -Converts a [`CoupledSDEs`](@ref) into [`CoupledODEs`](https://juliadynamics.github.io/DynamicalSystems.jl/stable/tutorial/#DynamicalSystemsBase.CoupledODEs) -from DynamicalSystems.jl. -""" -function DynamicalSystemsBase.CoupledODEs( - sys::CoupledSDEs; diffeq=DynamicalSystemsBase.DEFAULT_DIFFEQ, t0=0.0 -) - return DynamicalSystemsBase.CoupledODEs( - sys.integ.f, SVector{length(sys.integ.u)}(sys.integ.u), sys.p0; diffeq=diffeq, t0=t0 - ) -end - -# Pretty print -function additional_details(ds::CoupledSDEs) - solver, remaining = _decompose_into_solver_and_remaining(ds.diffeq) - return ["SDE solver" => string(nameof(typeof(solver))), "SDE kwargs" => remaining] -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) - -function DynamicalSystemsBase.set_parameter!(ds::CoupledSDEs, args...) - _set_parameter!(ds, args...) - u_modified!(ds.integ, true) - return nothing -end - -referrenced_sciml_prob(ds::CoupledSDEs) = ds.integ.sol.prob - -# so that `ds` is printed -set_state!(ds::CoupledSDEs, u::AbstractArray) = (set_state!(ds.integ, u); ds) -SciMLBase.step!(ds::CoupledSDEs, args...) = (step!(ds.integ, args...); ds) - -# For checking successful step, the `SciMLBase.step!` function checks -# `integ.sol.retcode in (ReturnCode.Default, ReturnCode.Success) || break`. -# But the actual API call would be `successful_retcode(check_error(integ))`. -# The latter however is already used in `step!(integ)` so there is no reason to re-do it. -# Besides, within DynamicalSystems.jl the integration is never expected to terminate. -# Nevertheless here we extend explicitly only for ODE stuff because it may be that for -# other type of DEIntegrators a different step interruption is possible. -function DynamicalSystemsBase.successful_step(integ::SciMLBase.AbstractSDEIntegrator) - rcode = integ.sol.retcode - return rcode == SciMLBase.ReturnCode.Default || rcode == SciMLBase.ReturnCode.Success -end diff --git a/src/CoupledSDEs_utils.jl b/src/CoupledSDEs_utils.jl deleted file mode 100644 index 0c393c67..00000000 --- a/src/CoupledSDEs_utils.jl +++ /dev/null @@ -1,62 +0,0 @@ -diagonal_cov(l) = Diagonal(ones(l)) - -""" -$(TYPEDSIGNATURES) - -Translates the stochastic process specified in `sys` into the language required by the -`SDEProblem` of `DynamicalSystems.jl`. -""" -function noise_process(sys::CoupledSDEs) - return sys.integ.W -end - -""" -$(TYPEDSIGNATURES) - -Gives the covariance matrix specified in `sys`. -""" -function covariance_matrix(sys::CoupledSDEs) - noise = noise_process(sys) - covariance = isnothing(noise) ? nothing : noise.covariance - if isnothing(noise) || isnothing(covariance) - return diagonal_cov(dimension(sys)) - else - return covariance - end -end - -""" -$(TYPEDSIGNATURES) - -Gives the noise strength specified in `sys`. -""" -function noise_strength(sys::CoupledSDEs) - return sys.noise_strength -end - -""" -$(TYPEDSIGNATURES) - -Returns the drift field ``b(x)`` of the CoupledSDEs `sys` at the state vector `x`. -""" -function drift(sys::CoupledSDEs{IIP}, x) where {IIP} - # assumes the drift is time independent - f = dynamic_rule(sys) - if IIP - dx = similar(x) - f(dx, x, sys.p0, 0) - return dx - else - return f(x, sys.p0, 0) - end -end - -""" -$(TYPEDSIGNATURES) - -Computes the divergence of the drift field `sys.f` at the given point `x`. -""" -function div_drift(sys::CoupledSDEs, x) - b(x) = drift(sys, x) - return tr(ForwardDiff.jacobian(b, x)) -end; diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index e7aca763..9b4d70ff 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -2,16 +2,10 @@ module CriticalTransitions # Base using Statistics: Statistics, mean -using LinearAlgebra: LinearAlgebra, Diagonal, I, norm, tr, dot +using LinearAlgebra: LinearAlgebra, I, norm, dot, tr using StaticArrays: StaticArrays, SVector -# core -using DynamicalSystemsBase: - DynamicalSystemsBase, - ContinuousTimeDynamicalSystem, - StateSpaceSets, - dimension, - dynamic_rule +# Core using DiffEqNoiseProcess: DiffEqNoiseProcess using OrdinaryDiffEq: OrdinaryDiffEq, Tsit5 using StochasticDiffEq: @@ -25,6 +19,9 @@ using StochasticDiffEq: step!, terminate!, u_modified! +using DynamicalSystemsBase: + DynamicalSystemsBase, CoupledSDEs, CoupledODEs, + dynamic_rule, current_state, set_state!, trajectory using ForwardDiff: ForwardDiff using IntervalArithmetic: IntervalArithmetic, interval @@ -51,8 +48,7 @@ using Reexport: @reexport include("extention_functions.jl") include("utils.jl") -include("CoupledSDEs.jl") -include("CoupledSDEs_utils.jl") +include("system_utils.jl") include("io.jl") include("trajectories/simulation.jl") include("trajectories/transition.jl") @@ -66,23 +62,18 @@ include("../systems/CTLibrary.jl") using .CTLibrary # Core types -export CoupledSDEs, - idfunc!, - idfunc, - add_noise_strength, - noise_process, - covariance_matrix, - noise_strength, - CoupledODEs +export CoupledSDEs, CoupledODEs, noise_process, covariance_matrix, diffusion_matrix +export dynamic_rule, current_state, set_state!, trajectory # Methods -export equilib, basins, basinboundary, basboundary -export simulate, relax +export drift, div_drift +export equilib, deterministic_orbit export transition, transitions +export basins, basinboundary, basboundary export fw_action, om_action, action, geometric_action export min_action_method, geometric_min_action_method export make_jld2, make_h5, intervals_to_box -# export basins, basinboundary +export covariance_matrix, diffusion_matrix # export edgetracking, bisect_to_edge, AttractorsViaProximity # export fixedpoints # ^ extention tests needed diff --git a/src/largedeviations/action.jl b/src/largedeviations/action.jl index b3cda6cf..bf8b7c63 100644 --- a/src/largedeviations/action.jl +++ b/src/largedeviations/action.jl @@ -2,25 +2,26 @@ $(TYPEDSIGNATURES) Calculates the Freidlin-Wentzell action of a given `path` with time points `time` in a -drift field specified by the deterministic dynamics of `sys`. +drift field specified by the deterministic dynamics `f = dynamic_rule(sys)` and +(normalized) noise covariance matrix `covariance_matrix(sys)`. The path must be a `(D x N)` matrix, where `D` is the dimensionality of the system `sys` and `N` is the number of path points. The `time` array must have length `N`. Returns a single number, which is the value of the action functional -``S_T[\\phi_t] = \\frac{1}{2} \\int_0^T || \\dot \\phi_t - b(\\phi_t) ||^2_Q dt`` +``S_T[\\phi_t] = \\frac{1}{2} \\int_0^T || \\dot \\phi_t - f(\\phi_t) ||^2_Q \\text{d}t`` where ``\\phi_t`` denotes the path in state space, ``b`` the drift field, and ``T`` the total time of the path. The subscript ``Q`` refers to the -generalized norm ``||a||_Q^2 := \\langle a, Q^{-1} b \\rangle`` (see `anorm``). Here -``Q`` is the noise covariance matrix `sys.Σ`. - +generalized norm ``||a||_Q^2 := \\langle a, Q^{-1} b \\rangle`` (see `anorm`). Here +``Q`` is the noise covariance matrix normalized by ``D/L_1(Q)``, with ``L_1`` being the +L1 matrix norm. """ function fw_action(sys::CoupledSDEs, path, time) @assert all(diff(time) .≈ diff(time[1:2])) "Freidlin-Wentzell action is only defined for equispaced time" # Inverse of covariance matrix - A = inv(covariance_matrix(sys)) # Do we have to take the inverse here? + A = inv(normalize_covariance!(covariance_matrix(sys))) # Compute action integral integrand = fw_integrand(sys, path, time, A) @@ -33,27 +34,28 @@ function fw_action(sys::CoupledSDEs, path, time) end; """ -$(TYPEDSIGNATURES) + om_action(sys::CoupledSDEs, path, time, noise_strength) -Calculates the Onsager-Machlup action of a given `path` with time points `time` for the drift field `sys.f` at given `sys.noise_strength`. +Calculates the Onsager-Machlup action of a given `path` with time points `time` for the drift field `f = dynamic_rule(sys)` at given `noise_strength`. The path must be a `(D x N)` matrix, where `D` is the dimensionality of the system `sys` and `N` is the number of path points. The `time` array must have length `N`. Returns a single number, which is the value of the action functional -``I^{\\sigma}_T[\\phi_t] = \\frac{1}{2} \\int_0^T \\left( || \\dot \\phi_t - b(\\phi_t) ||^2_Q + -\\frac{\\sigma^2}{2} \\div(b) \\right) \\, dt`` +``S^{\\sigma}_T[\\phi_t] = \\frac{1}{2} \\int_0^T \\left( || \\dot \\phi_t - f(\\phi_t) ||^2_Q + +\\sigma^2 \\nabla \\cdot f \\right) \\, \\text{d} t`` where ``\\phi_t`` denotes the path in state space, ``b`` the drift field, ``T`` the total time of the path, and ``\\sigma`` the noise strength. The subscript ``Q`` refers to the -generalized norm ``||a||_Q^2 := \\langle a, Q^{-1} b \\rangle`` (see `anorm``). Here -``Q`` is the noise covariance matrix. +generalized norm ``||a||_Q^2 := \\langle a, Q^{-1} b \\rangle`` (see `anorm`). Here +``Q`` is the noise covariance matrix normalized by ``D/L_1(Q)``, with ``L_1`` being the +L1 matrix norm. """ -function om_action(sys::CoupledSDEs, path, time) +function om_action(sys::CoupledSDEs, path, time, noise_strength) @assert all(diff(time) .≈ diff(time[1:2])) "Fw_action is only defined for equispaced time" - σ = noise_strength(sys) + σ = noise_strength # Compute action integral S = 0 for i in 1:(size(path, 2) - 1) @@ -75,11 +77,11 @@ Computes the action functional specified by `functional` for a given CoupledSDEs * `functional = "FW"`: Returns the Freidlin-Wentzell action ([`fw_action`](@ref)) * `functional = "OM"`: Returns the Onsager-Machlup action ([`om_action`](@ref)) """ -function action(sys::CoupledSDEs, path::Matrix, time, functional) +function action(sys::CoupledSDEs, path::Matrix, time, functional; noise_strength=nothing) if functional == "FW" action = fw_action(sys, path, time) elseif functional == "OM" - action = om_action(sys, path, time) + action = om_action(sys, path, time, noise_strength) end return action end; @@ -88,26 +90,29 @@ end; $(TYPEDSIGNATURES) Calculates the geometric action of a given `path` with specified `arclength` for the drift -field `sys.f`. +field specified by the deterministic dynamics `f = dynamic_rule(sys)` and +(normalized) noise covariance matrix `covariance_matrix(sys)`. For a given path ``\\varphi``, the geometric action ``\\bar S`` corresponds to the minimum -of the Freidlin-Wentzell action ``S_T(\\phi)`` over all travel times ``T>0``, where ``\\phi`` +of the Freidlin-Wentzell action ``S_T(\\varphi)`` over all travel times ``T>0``, where ``\\varphi`` denotes the path's parameterization in physical time (see [`fw_action`](@ref)). It is given by the integral -``\\bar S[\\varphi] = \\int_0^L \\left( ||\\varphi'||_Q \\, ||b(\\varphi)||_Q - \\langle \\varphi', \\, - b(\\varphi) \\rangle_Q \\right) \\, ds`` +``\\bar S[\\varphi] = \\int_0^L \\left( ||\\varphi'||_Q \\, ||f(\\varphi)||_Q - \\langle \\varphi', \\, + f(\\varphi) \\rangle_Q \\right) \\, \\text{d}s`` -where ``s`` is the arclength coordinate, ``L`` the arclength, ``b`` the drift field, and the +where ``s`` is the arclength coordinate, ``L`` the arclength, ``f`` the drift field, and the subscript ``Q`` refers to the generalized dot product ``\\langle a, b \\rangle_Q := a^{\\top} -\\cdot Q^{-1} b`` (see `anorm``). Here ``Q`` is the noise covariance matrix `sys.Σ`. +\\cdot Q^{-1} b`` (see `anorm`). Here +``Q`` is the noise covariance matrix normalized by ``D/L_1(Q)``, with ``L_1`` being the +L1 matrix norm. Returns the value of the geometric action ``\\bar S``. """ function geometric_action(sys::CoupledSDEs, path, arclength=1.0) N = size(path, 2) v = path_velocity(path, range(0, arclength; length=N); order=4) - A = inv(covariance_matrix(sys)) + A = inv(normalize_covariance!(covariance_matrix(sys))) b(x) = drift(sys, x) diff --git a/src/system_utils.jl b/src/system_utils.jl new file mode 100644 index 00000000..5917d222 --- /dev/null +++ b/src/system_utils.jl @@ -0,0 +1,47 @@ +StochasticSystemsBase = Base.get_extension(DynamicalSystemsBase, :StochasticSystemsBase) +covariance_matrix = StochasticSystemsBase.covariance_matrix +diffusion_matrix = StochasticSystemsBase.diffusion_matrix + +""" + StochSystem + +An alias to [`CoupledSDEs`](@ref). +This was the name these systems had in CriticalTransitions.jl before v0.3 +""" +const StochSystem = CoupledSDEs + +""" +$(TYPEDSIGNATURES) + +Fetches the stochastic process ``\\mathcal{N}`` specified in the intergrator of `sys`. Returns the type `DiffEqNoiseProcess.NoiseProcess`. +""" +function noise_process(sys::CoupledSDEs) + return sys.integ.W +end + +""" +$(TYPEDSIGNATURES) + +Returns the deterministic drift ``f(x)`` of the CoupledSDEs `sys` at state `x`. For time-dependent systems, the time can be specified as a keyword argument `t` (by default `t=0`). +""" +function drift(sys::CoupledSDEs{IIP}, x; t=0) where {IIP} + f = dynamic_rule(sys) + if IIP + dx = similar(x) + f(dx, x, sys.p0, t) + return dx + else + return f(x, sys.p0, t) + end +end + +""" +$(TYPEDSIGNATURES) + +Computes the divergence of the drift field ``f(x)`` at state `x`. For time- +dependent systems, the time can be specified as a keyword argument `t` (by default `t=0`). +""" +function div_drift(sys::CoupledSDEs, x; t=0) + b(x) = drift(sys, x; t) + return tr(ForwardDiff.jacobian(b, x)) +end; diff --git a/src/trajectories/equib.jl b/src/trajectories/equib.jl index 5e275463..a89a1cab 100644 --- a/src/trajectories/equib.jl +++ b/src/trajectories/equib.jl @@ -1,5 +1,5 @@ """ - equilib(sys::CoupledSDEs, state; kwargs...) + equilib(sys::CoupledSDEs [, state]; kwargs...) Returns the equilibrium solution of the system `sys` for given initial condition `state`. > Warning: This algorithm simply evolves the deterministic system forward in time until a steady-state condition is satisfied. @@ -12,13 +12,14 @@ Returns the equilibrium solution of the system `sys` for given initial condition * `dt = 0.01`: time step of the ODE solver. * `solver = Euler()`: ODE solver used for evolving the state. """ -function equilib(sys::CoupledSDEs, state; dt=0.01, tmax=1e5, abstol=1e-5, solver=Tsit5()) +function equilib( + sys::CoupledSDEs, state=nothing; dt=0.01, tmax=1e5, abstol=1e-5, solver=Tsit5() +) condition(u, t, integrator) = norm(integrator.uprev - u) < abstol affect!(integrator) = terminate!(integrator) equilib_cond = DiscreteCallback(condition, affect!) - - prob = ODEProblem(sys.integ.f, state, (0, tmax), sys.p0) + isnothing(state) ? nothing : set_state!(sys, state) + prob = sys.integ.sol.prob sol = solve(prob, solver; dt=dt, callback=equilib_cond, save_on=false, save_start=false) - return sol.u[1] end; diff --git a/src/trajectories/simulation.jl b/src/trajectories/simulation.jl index e8f03cdf..57c31195 100644 --- a/src/trajectories/simulation.jl +++ b/src/trajectories/simulation.jl @@ -1,41 +1,22 @@ """ $(TYPEDSIGNATURES) -Simulates the CoupledSDEs `sys` forward in time for a duration `T`, starting at initial condition `init`. +Simulates the deterministic (noise-free) dynamics of CoupledSDEs `sys` in time for a +duration `T`, starting at initial condition `init`. -This function uses the [`SDEProblem`](https://diffeq.sciml.ai/stable/types/sde_types/#SciMLBase.SDEProblem) functionality of `DifferentialEquations.jl`. +This function is equivalent to calling [`trajectory`](@ref) on the deterministic part of the +`CoupledSDEs` (with `noise_strength=0`). It works with the ODE solvers used for `CoupledODEs`. ## Keyword arguments -* `alg=SOSRA()`: [SDE solver](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/#sde_solve). Defaults to an adaptive strong order 1.5 method -* `kwargs...`: keyword arguments for [`solve(SDEProblem)`](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/#solver_options) - -For more info, see [`SDEProblem`](https://diffeq.sciml.ai/stable/types/sde_types/#SciMLBase.SDEProblem). - -""" -function simulate( - sys::CoupledSDEs, T, init=current_state(sys); alg=sys.integ.alg, kwargs... -) - prob = remake(referrenced_sciml_prob(sys); u0=init, tspan=(0, T)) - return solve(prob, alg; kwargs...) -end; - -""" -$(TYPEDSIGNATURES) - -Simulates the deterministic dynamics of CoupledSDEs `sys` in time for a duration `T`, starting at initial condition `init`. - -This function integrates `sys.f` forward in time, using the [`ODEProblem`](https://diffeq.sciml.ai/stable/types/ode_types/#SciMLBase.ODEProblem) functionality of `DifferentialEquations.jl`. Thus, `relax` is identical to [`simulate`](@ref) when setting the noise strength `sys.σ = 0`. - -## Keyword arguments -* `solver=Tsit5()`: [ODE solver](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#ode_solve). Defaults to Tsitouras 5/4 Runge-Kutta method -* `kwargs...`: keyword arguments for [`solve(ODEProblem)`](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/#solver_options) +* `diffeq=(alg=Tsit5(), abstol = 1e-6, reltol = 1e-6)`: ODE solver settings (see [`CoupledODEs`](@ref)) +* `kwargs...`: keyword arguments passed to [`trajectory`](@ref) For more info, see [`ODEProblem`](https://diffeq.sciml.ai/stable/types/ode_types/#SciMLBase.ODEProblem). -For stochastic integration, see [`simulate`](@ref). +For stochastic integration, see [`trajectory`](@ref). """ -function relax(sys::CoupledSDEs, T, init=current_state(sys); alg=Tsit5(), kwargs...) - sde_prob = referrenced_sciml_prob(sys) - prob = ODEProblem{isinplace(sde_prob)}(dynamic_rule(sys), init, (0, T), sys.p0) - return solve(prob, alg; kwargs...) +function deterministic_orbit( + sys::CoupledSDEs, T, init=current_state(sys); diffeq=CoupledODEs(sys).diffeq, kwargs... +) + return trajectory(CoupledODEs(sys; diffeq), T, init; kwargs...) end; diff --git a/src/trajectories/transition.jl b/src/trajectories/transition.jl index 546a8f97..4a490fff 100644 --- a/src/trajectories/transition.jl +++ b/src/trajectories/transition.jl @@ -26,23 +26,20 @@ This function simulates `sys` in time, starting from initial condition `x_i`, un ## Keyword arguments * `rad_i=0.1`: radius of ball around `x_i` * `rad_f=0.1`: radius of ball around `x_f` -* `cut_start=true`: if `false`, returns the whole trajectory up to the transition -* `dt=0.01`: time step of integration * `tmax=1e3`: maximum time when the simulation stops even `x_f` has not been reached * `rad_dims=1:length(sys.u)`: the directions in phase space to consider when calculating the radii `rad_i` and `rad_f`. Defaults to all directions. To consider only a subspace of state space, insert a vector of indices of the dimensions to be included. -* `solver=EM()`: numerical solver. Defaults to Euler-Mayurama. -* `progress`: shows a progress bar with respect to `tmax` +* `cut_start=true`: if `false`, returns the whole trajectory up to the transition +* `kwargs...`: keyword arguments passed to [`CommonSolve.solve`](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/#CommonSolve.solve-Tuple{SciMLBase.AbstractDEProblem,%20Vararg{Any}}) ## Output `[path, times, success]` * `path` (Matrix): transition path (size [dim × N], where N is the number of time points) * `times` (Vector): time values (since start of simulation) of the path points (size N) * `success` (bool): if `true`, a transition occured (i.e. the ball around `x_f` has been reached), else `false` -* `kwargs...`: keyword arguments passed to [`simulate`](@ref) -See also [`transitions`](@ref), [`simulate`](@ref). +See also [`transitions`](@ref), [`trajectory`](@ref). """ function transition( sys::CoupledSDEs, @@ -51,15 +48,16 @@ function transition( rad_i=0.1, rad_f=0.1, tmax=1e3, - cut_start=true, rad_dims=1:length(current_state(sys)), + cut_start=true, kwargs..., ) condition(u, t, integrator) = subnorm(u - x_f; directions=rad_dims) < rad_f affect!(integrator) = terminate!(integrator) cb_ball = DiscreteCallback(condition, affect!) - - sim = simulate(sys, tmax, x_i; callback=cb_ball, kwargs...) + + prob = remake(sys.integ.sol.prob; u0=x_i, tspan=(0, tmax)) + sim = solve(prob, sys.integ.alg; callback=cb_ball, kwargs...) success = sim.retcode == SciMLBase.ReturnCode.Terminated simt = sim.t @@ -81,7 +79,7 @@ function transition( end; """ -$(TYPEDSIGNATURES) + function transitions(sys::CoupledSDEs, x_i, x_f, N=1; kwargs...) Generates an ensemble of `N` transition samples of `sys` from point `x_i` to point `x_f`. @@ -90,15 +88,15 @@ This function repeatedly calls the [`transition`](@ref) function to efficiently ## Keyword arguments - `rad_i=0.1`: radius of ball around `x_i` - `rad_f=0.1`: radius of ball around `x_f` - - `cut_start=true`: if `false`, returns the whole trajectory up to the transition - `Nmax`: number of attempts before the algorithm stops even if less than `N` transitions occurred. - - `dt=0.01`: time step of integration - `tmax=1e3`: maximum time when the simulation stops even `x_f` has not been reached - `rad_dims=1:length(sys.u)`: the directions in phase space to consider when calculating the radii `rad_i` and `rad_f`. Defaults to all directions. To consider only a subspace of state space, insert a vector of indices of the dimensions to be included. - - `progress`: shows a progress bar with respect to `Nmax` + - `cut_start=true`: if `false`, returns the whole trajectory up to the transition - `savefile`: if `nothing`, no data is saved to a file. To save to a file, see below. + - `showprogress`: shows a progress bar with respect to `Nmax` + - `kwargs...`: keyword arguments passed to [`CommonSolve.solve`](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/#CommonSolve.solve-Tuple{SciMLBase.AbstractDEProblem,%20Vararg{Any}}) See also [`transition`](@ref). diff --git a/src/utils.jl b/src/utils.jl index e6f2e45e..e62dd15b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -103,6 +103,19 @@ function subnorm(vec; directions=1:length(vec)) return sqrt(sum) end; +""" +$(TYPEDSIGNATURES) + +Normalizes the covariance matrix ``Q`` (in-place) by dividing it by +``L_1(Q)/\\text{dim}(Q)``, where ``L_1(Q)`` is the L1 matrix norm of ``Q`` and +``\\text{dim}(Q)`` is the dimension of ``Q``. +""" +function normalize_covariance!(covariance) + l1norm = norm(covariance, 1) + dim = size(covariance)[1] + return covariance * dim/l1norm +end + # Central finite difference, second derivative function central(f, idx, dz) return (f[idx + 1] - f[idx - 1]) / (2 * dz) diff --git a/test/CoupledSDEs.jl b/test/CoupledSDEs.jl index acab780e..dda8fc2f 100644 --- a/test/CoupledSDEs.jl +++ b/test/CoupledSDEs.jl @@ -1,161 +1,40 @@ -@testset "initialisation" begin - function diagonal_noise!(σ) - function (du, u, p, t) - idfunc!(du, u, p, t) - du .*= σ - return nothing - end - end - diagonal_noise(σ) = (u, p, t) -> σ .* idfunc(u, p, t) - - # Define the system - function meier_stein(u, p, t) # out-of-place +@testset "API" begin + using LinearAlgebra + function fitzhugh_nagumo(u, p, t) x, y = u - dx = x - x^3 - 10 * x * y^2 - dy = -(1 + x^2) * y - return SA[dx, dy] - end - σ = 0.1 - - @testset "diagonal additive noise" begin - # diagonal additive noise: σ*N(0,dt) - # a vector of random numbers dW whose size matches the output of g where the noise is applied element-wise - prob = SDEProblem(meier_stein, diagonal_noise(σ), zeros(2), (0.0, Inf), ()) - sde = CoupledSDEs(meier_stein, idfunc, zeros(2), (), σ) + ϵ, β, α, γ, κ, I = p - @test sde.integ.sol.prob.f isa SDEFunction - @test sde.integ.sol.prob.f.f.f == prob.f.f - # @test sde.integ.sol.prob.g == prob.g - @test sde.integ.sol.prob.u0 == prob.u0 - @test sde.integ.sol.prob.tspan == prob.tspan - @test sde.integ.sol.prob.p == prob.p - @test sde.integ.sol.prob.noise_rate_prototype == prob.noise_rate_prototype - @test sde.integ.sol.prob.noise == prob.noise - @test sde.integ.sol.prob.seed == prob.seed + dx = (-α * x^3 + γ * x - κ * y + I) / ϵ + dy = -β * y + x - noise = sde.integ.sol.prob.noise - @test noise == nothing - end - # @show sde.integ.W.dist - # @show sde.integ.sol.prob - - @testset "Scalar noise Wiener" begin - # Scalar noise Wiener - # a single random variable is applied to all dependent variables - f!(du, u, p, t) = du .= 1.01u # deterministic part - σ = 0.25 # noise strength - W = WienerProcess(0.0, 0.0, 0.0) - prob = SDEProblem(f!, diagonal_noise!(σ), zeros(2), (0.0, Inf), (); noise=W) - sde = CoupledSDEs(f!, idfunc!, zeros(2), (), σ; noise=W) - - @test sde.integ.sol.prob.f isa SDEFunction - @test sde.integ.sol.prob.f.f.f == prob.f.f - # @test sde.integ.sol.prob.g == prob.g - @test sde.integ.sol.prob.u0 == prob.u0 - @test sde.integ.sol.prob.tspan == prob.tspan - @test sde.integ.sol.prob.p == prob.p - @test sde.integ.sol.prob.noise_rate_prototype == prob.noise_rate_prototype - @test sde.integ.sol.prob.noise == prob.noise - @test sde.integ.sol.prob.seed == prob.seed - - noise = sde.integ.sol.prob.noise - @test noise == W - @test length(noise.dW) == 1 - @test W.covariance == nothing - - sol = simulate(sde, 1.0; dt=0.01, alg=SOSRA()) - # every vatiable has the same noise and dynamic_rule - @test all(sol[1, :] .≈ sol[2, :]) - end - - @testset "multiplicative noise Wiener" begin - # multiplicative noise Wiener - # a single random variable is applied to all dependent variables - g_sde(u, p, t) = σ .* u - g_Csde(u, p, t) = σ .* u - prob = SDEProblem(meier_stein, g_sde, zeros(2), (0.0, Inf), ()) - sde = CoupledSDEs(meier_stein, g_Csde, zeros(2), (), σ) - - @test sde.integ.sol.prob.f isa SDEFunction - @test sde.integ.sol.prob.f.f.f == prob.f.f - @test sde.integ.sol.prob.u0 == prob.u0 - @test sde.integ.sol.prob.tspan == prob.tspan - @test sde.integ.sol.prob.p == prob.p - @test sde.integ.sol.prob.noise_rate_prototype == prob.noise_rate_prototype - @test sde.integ.sol.prob.noise == prob.noise - @test sde.integ.sol.prob.seed == prob.seed - - noise = sde.integ.sol.prob.noise - @test noise == nothing + return SA[dx, dy] end - @testset "Non-diagonal noise" begin - # non-diagonal noise allows for the terms to linearly mixed via g being a matrix. - # four Wiener processes and two dependent random variables - # the output of g to be a 2x4 matrix, such that the solution is g(u,p,t)*dW, the matrix multiplication. - f(du, u, p, t) = du .= 1.01u - function g(du, u, p, t) - du[1, 1] = 0.3u[1] - du[1, 2] = 0.6u[1] - du[1, 3] = 0.9u[1] - du[1, 4] = 0.12u[1] - du[2, 1] = 1.2u[2] - du[2, 2] = 0.2u[2] - du[2, 3] = 0.3u[2] - du[2, 4] = 1.8u[2] - return nothing - end - # prob = SDEProblem(f, g, ones(2), (0.0, 1.0), noise_rate_prototype = zeros(2, 4)) - # diffeq =(alg=SRIW1(), noise_rate_prototype = zeros(2, 4)) - sde = CoupledSDEs( - f, g, zeros(2); noise_rate_prototype=zeros(2, 4), diffeq=(alg=RKMilCommute(),) - ) - prob = SDEProblem(f, g, zeros(2), (0.0, Inf); noise_rate_prototype=zeros(2, 4)) - - @test sde.integ.sol.prob.f isa SDEFunction - @test sde.integ.sol.prob.f.f.f == prob.f.f - # @test sde.integ.sol.prob.g == prob.g - @test sde.integ.sol.prob.u0 == prob.u0 - @test sde.integ.sol.prob.tspan == prob.tspan - @test sde.integ.sol.prob.p == prob.p - @test sde.integ.sol.prob.noise_rate_prototype == prob.noise_rate_prototype - @test sde.integ.sol.prob.noise == prob.noise - @test sde.integ.sol.prob.seed == prob.seed + σ = 0.2 + param = [1.0, 3, 1, 1, 1, 0] + u0 = zeros(2) + sys = CoupledSDEs(fitzhugh_nagumo, u0, param; noise_strength=σ) + + @testset "correlation" begin end + + @testset "noise" begin + diff = diffusion_matrix(sys) + cov = covariance_matrix(sys) + @test diff isa AbstractMatrix + @test cov isa AbstractMatrix + @test all(diag(diff) .== σ) + @test diff .^ 2 == cov + + W = noise_process(sys) + int_cov = W.covariance + # The internal covariance in the DiffEqNoiseProcess.NoiseProcess should be nothing + @test int_cov == nothing end - @testset "Correlated noise Wiener" begin - f!(du, u, p, t) = du .= 1.01u - - ρ = 0.3 - Γ = [1 ρ; ρ 1] - t0 = 0.0 - W0 = zeros(2) - Z0 = zeros(2) - W = CorrelatedWienerProcess(Γ, t0, W0, Z0) - prob = SDEProblem(f!, diagonal_noise!(σ), zeros(2), (0.0, Inf), (); noise=W) - sde = CoupledSDEs(f!, idfunc!, zeros(2), (), σ; noise=W) - - @test sde.integ.sol.prob.f isa SDEFunction - @test sde.integ.sol.prob.f.f.f == prob.f.f - @test sde.integ.sol.prob.u0 == prob.u0 - @test sde.integ.sol.prob.tspan == prob.tspan - @test sde.integ.sol.prob.p == prob.p - @test sde.integ.sol.prob.noise_rate_prototype == prob.noise_rate_prototype - @test sde.integ.sol.prob.noise == prob.noise - @test sde.integ.sol.prob.seed == prob.seed - - @test W.covariance == Γ - @testset "covariance" begin - using DiffEqNoiseProcess, Statistics - prob = NoiseProblem(W, (0.0, 1.0)) - output_func = (sol, i) -> (sol.dW, false) - ensemble_prob = EnsembleProblem(prob; output_func=output_func) - - dt = 0.1 - sol = Array(solve(ensemble_prob; dt=dt, trajectories=1_000_000)) - - @test zero(mean(sol; dims=2)[:]) ≈ mean(sol; dims=2)[:] atol = 1e-2 - @test W.covariance ≈ cov(sol; dims=2) / dt rtol = 1e-2 - end + @testset "drift" begin + using CriticalTransitions: drift + drift_vector = drift(sys, [0, 0]) + drift_vector isa SVector{2,Float64} + @test drift_vector == [0, 0] end end diff --git a/test/basin/basin_boundary.jl b/test/basin/basin_boundary.jl new file mode 100644 index 00000000..e2244108 --- /dev/null +++ b/test/basin/basin_boundary.jl @@ -0,0 +1,18 @@ +@testset "fitzhugh_nagumo" begin + # Define systems + g = diag_noise_funtion(0.215) + system = CoupledSDEs(fitzhugh_nagumo, g, [2.0, 0.1], [0.1, 3.0, 1.0, 1.0, 1.0, 0.0]) + + # Set up domain + A, B, C = [0.0, 0.0], [1.0, 0.0], [0.0, 1.0] + box = intervals_to_box([-2.0, -2.0], [2.00, 2.0]) + step = 0.1 + + boas = basins(system, A, B, C, box; bstep=[step, step]) + X, Y, attractors, h = boas + boundary = basinboundary(boas) + @test length(attractors) == 2 + @test length(unique(h)) == 3 + @test h[end ÷ 2 + 1, end ÷ 2 + 1] == -1 + @test boundary[:, end ÷ 2 + 1] == zeros(2) +end diff --git a/test/ext/ChaosToolsExt.jl b/test/ext/ChaosToolsExt.jl index 54991019..9f2ee3c2 100644 --- a/test/ext/ChaosToolsExt.jl +++ b/test/ext/ChaosToolsExt.jl @@ -9,7 +9,7 @@ function meier_stein(u, p, t) # out-of-place end σ = 0.1 -sde = CoupledSDEs(meier_stein, idfunc, zeros(2), (), σ) +sde = CoupledSDEs(meier_stein, zeros(2), (); noise_strength=σ) fps, eigs, stab = fixedpoints(sde, [-3, -3], [3, 3]) diff --git a/test/ext/CoupledSDEsBaisin.jl b/test/ext/CoupledSDEsBaisin.jl index 152288c2..944d35bb 100644 --- a/test/ext/CoupledSDEsBaisin.jl +++ b/test/ext/CoupledSDEsBaisin.jl @@ -2,7 +2,7 @@ using Attractors, ChaosTools # Define systems system = CoupledSDEs( - fitzhugh_nagumo, idfunc, [2.0, 0.1], [0.1, 3.0, 1.0, 1.0, 1.0, 0.0], 0.215 + fitzhugh_nagumo, [2.0, 0.1], [0.1, 3.0, 1.0, 1.0, 1.0, 0.0]; noise_strength=0.215 ) # Set up domain diff --git a/test/issue14.jl b/test/issue14.jl new file mode 100644 index 00000000..a8aac0f0 --- /dev/null +++ b/test/issue14.jl @@ -0,0 +1,22 @@ +using CriticalTransitions, StaticArrays, Test + +function meier_stein(u, p, t) # in-place + x, y = u + dx = x - x^3 - 10 * x * y^2 + dy = -(1 + x^2) * y + return SA[dx, dy] +end +σ = 0.25 +sys = StochSystem(meier_stein, [], zeros(2), σ, idfunc, nothing, I(2), "WhiteGauss") + +# initial path: parabola +xx = range(-1.0, 1.0; length=30) +yy = 0.3 .* (-xx .^ 2 .+ 1) +init = Matrix([xx yy]') + +gm = geometric_min_action_method(sys, init; maxiter=1000) +@test all(isapprox.(path[2, :][(end - 5):end], 0, atol=1e-3)) + +gm[2][end] + +isapprox(gm[2][end], 0.338; atol=1e-3) diff --git a/test/largedeviations/MAM.jl b/test/largedeviations/MAM.jl index 461ec3e0..daef7e39 100644 --- a/test/largedeviations/MAM.jl +++ b/test/largedeviations/MAM.jl @@ -1,7 +1,7 @@ @testset "MAM FitzHugh-Nagumo" begin p = [0.1, 3, 1, 1, 1, 0] σ = 0.1 - fhn = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ) + fhn = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=σ) x_i = SA[sqrt(2 / 3), sqrt(2 / 27)] x_f = SA[0.001, 0.0] N, T = 200, 2.0 @@ -14,7 +14,7 @@ end @testset "MAM Ornstein-Uhlenbeck" begin σ = 0.18 - ou = CoupledSDEs((u, p, t) -> -u, idfunc, SA[1.0], (), σ) + ou = CoupledSDEs((u, p, t) -> -u, SA[1.0]; noise_strength=σ) x0 = -1 xT = 2.0 T = 10.0 diff --git a/test/largedeviations/action_fhn.jl b/test/largedeviations/action_fhn.jl index 333d1477..7f0682ec 100644 --- a/test/largedeviations/action_fhn.jl +++ b/test/largedeviations/action_fhn.jl @@ -4,9 +4,9 @@ Random.seed!(SEED) # System setup - FitzHugh-Nagumo model p = [1.0, 3.0, 1.0, 1.0, 1.0, 0.0] # Parameters (ϵ, β, α, γ, κ, I) σ = 0.2 # noise strength -sys = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ; seed=SEED) +sys = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=σ) -A = inv(CT.covariance_matrix(sys)) +A = inv(covariance_matrix(sys)) T, N = 2.0, 100 x_i = SA[sqrt(2 / 3), sqrt(2 / 27)] @@ -22,7 +22,7 @@ end # Test om_action function @testset "om_action" begin - S = om_action(sys, path, time) + S = om_action(sys, path, time, σ) @test isapprox(S, 0.26, atol=0.01) end @@ -40,11 +40,12 @@ end # Test fw_integrand function @testset "fw_integrand" begin integrand = CriticalTransitions.fw_integrand(sys, path, time, A) - @test minimum(integrand) > 0.18 + @test minimum(integrand) > 0.18 skip = true + #^ This test is does not test something meaningful end # Test div_drift function @testset "div_drift" begin @test CT.div_drift(sys, zeros(2)) == -2.0 @test CT.div_drift(sys, x_i) == -4.0 -end +end \ No newline at end of file diff --git a/test/largedeviations/gMAM.jl b/test/largedeviations/gMAM.jl index 43b9a0ae..fb9a1a3d 100644 --- a/test/largedeviations/gMAM.jl +++ b/test/largedeviations/gMAM.jl @@ -1,7 +1,7 @@ @testset "gMAM FitzHugh-Nagumo" begin p = [0.1, 3, 1, 1, 1, 0] σ = 0.1 - fhn = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ) + fhn = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=σ) x_i = SA[sqrt(2 / 3), sqrt(2 / 27)] x_f = SA[0.001, 0.0] N = 100 @@ -19,7 +19,7 @@ end end σ = 0.25 # sys = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ) - sys = CoupledSDEs(meier_stein, idfunc, zeros(2), (), σ) + sys = CoupledSDEs(meier_stein, zeros(2); noise_strength=σ) # initial path: parabola xx = range(-1.0, 1.0; length=30) diff --git a/test/runtests.jl b/test/runtests.jl index 56e572c9..a6b4e389 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -58,13 +58,13 @@ end include("largedeviations/gMAM.jl") end -@testset "utilities" begin +@testset "Utilities" begin include("utils.jl") end -@testset "Trajactories" begin - include("trajactories/simulate.jl") - include("trajactories/transition.jl") +@testset "Trajectories" begin + include("trajectories/simulate.jl") +# include("trajectories/transition.jl") end @testset "Extentions" begin diff --git a/test/trajactories/simulate.jl b/test/trajactories/simulate.jl deleted file mode 100644 index 0e2dfd41..00000000 --- a/test/trajactories/simulate.jl +++ /dev/null @@ -1,18 +0,0 @@ -@testset "Simulate FitzHugh-Nagumo model" begin - # System setup - p = [0.1, 3.0, 1.0, 1.0, 1.0, 0.0] # Parameters (ϵ, β, α, γ, κ, I) - σ = 0.1 - # Simulate - sys = CoupledSDEs(fitzhugh_nagumo, idfunc, SA[1.0, 0.0], p, σ) - traj = relax(sys, 10) - sys = CoupledSDEs(fitzhugh_nagumo, idfunc, SA[1.0, 0.0], p, σ) - sim = simulate(sys, 10) - @test traj[1][1, 1] == 1.0 - @test sim.u[1][1] == 1.0 - # These tests could be improved - Reyk - - # using DynamicalSystemsBase - # sys = CoupledSDEs(fitzhugh_nagumo, idfunc, SA[1.0, 0.0], p, σ) - # trajectory(sys, 10) # works - # works? -end diff --git a/test/trajectories/simulate.jl b/test/trajectories/simulate.jl new file mode 100644 index 00000000..ac486ce4 --- /dev/null +++ b/test/trajectories/simulate.jl @@ -0,0 +1,13 @@ +@testset "Simulate FitzHugh-Nagumo model" begin + # System setup + p = [0.1, 3.0, 1.0, 1.0, 1.0, 0.0] # Parameters (ϵ, β, α, γ, κ, I) + σ = 0.1 + init = SA[1.0, 0.0] + sys = CoupledSDEs(fitzhugh_nagumo, init, p; noise_strength=σ) + traj = trajectory(sys, 10, init) + relax = deterministic_orbit(sys, 10, init) + + @test traj[1][1,1] == 1.0 + @test isapprox(relax[1][end,1], 0.816; atol=1e-2) + # These tests could be improved - Reyk +end diff --git a/test/trajactories/transition.jl b/test/trajectories/transition.jl similarity index 78% rename from test/trajactories/transition.jl rename to test/trajectories/transition.jl index 15ebd4b9..14c1fe25 100644 --- a/test/trajactories/transition.jl +++ b/test/trajectories/transition.jl @@ -7,7 +7,7 @@ σ = 0.215 # noise strength # CoupledSDEs - sys = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ; seed=SEED) + sys = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=σ, seed=SEED) fp1 = [0.816, 0.272] fp2 = [-0.816, -0.272] @@ -22,7 +22,6 @@ @test ensemble.t_trans ≈ 4.493941793363376 atol = 1e-2 # SEED is different on github # SEED doesn;t work on github - @test length(ensemble.times) == 11 broken=true - @test ensemble.t_res ≈ 5299.98 broken=true - + @test length(ensemble.times) == 11 broken = true + @test ensemble.t_res ≈ 5299.98 broken = true end diff --git a/test/utils.jl b/test/utils.jl index a0e420a8..a5e38593 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,23 +1,3 @@ -# Test for idfunc -@testset "idfunc" begin - u = [1, 2, 3] - p = [0.1, 0.2, 0.3] - t = 0.5 - expected = [1, 1, 1] - @test idfunc(u, p, t) == expected -end - -# Test for idfunc! -@testset "idfunc!" begin - du = zeros(3) - u = [1, 2, 3] - p = [0.1, 0.2, 0.3] - t = 0.5 - expected = [1, 1, 1] - idfunc!(du, u, p, t) - @test du == expected -end - # Test for intervals_to_box @testset "intervals_to_box" begin using IntervalArithmetic From 63b9f82bae8e958d0702767ff5e1139c1abb6976 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sat, 12 Oct 2024 16:26:05 +0200 Subject: [PATCH 3/3] run formatter --- docs/make.jl | 4 ++-- src/CriticalTransitions.jl | 9 +++++++-- src/trajectories/transition.jl | 2 +- src/utils.jl | 2 +- test/largedeviations/action_fhn.jl | 2 +- test/runtests.jl | 2 +- test/trajectories/simulate.jl | 6 +++--- 7 files changed, 16 insertions(+), 11 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 9bda4a59..3057fc92 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -47,7 +47,7 @@ makedocs(; Base.get_extension(CriticalTransitions, :ChaosToolsExt), Base.get_extension(CriticalTransitions, :CoupledSDEsBaisin), DynamicalSystemsBase, - Base.get_extension(DynamicalSystemsBase, :StochasticSystemsBase) + Base.get_extension(DynamicalSystemsBase, :StochasticSystemsBase), ], doctest=false, format=Documenter.HTML(; html_options...), @@ -56,4 +56,4 @@ makedocs(; plugins=[bib, links], ) -deploydocs(; repo="github.com/JuliaDynamics/CriticalTransitions.jl.git", push_preview=false) \ No newline at end of file +deploydocs(; repo="github.com/JuliaDynamics/CriticalTransitions.jl.git", push_preview=false) diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 9b4d70ff..cab1d333 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -20,8 +20,13 @@ using StochasticDiffEq: terminate!, u_modified! using DynamicalSystemsBase: - DynamicalSystemsBase, CoupledSDEs, CoupledODEs, - dynamic_rule, current_state, set_state!, trajectory + DynamicalSystemsBase, + CoupledSDEs, + CoupledODEs, + dynamic_rule, + current_state, + set_state!, + trajectory using ForwardDiff: ForwardDiff using IntervalArithmetic: IntervalArithmetic, interval diff --git a/src/trajectories/transition.jl b/src/trajectories/transition.jl index 4a490fff..d0e95cdd 100644 --- a/src/trajectories/transition.jl +++ b/src/trajectories/transition.jl @@ -55,7 +55,7 @@ function transition( condition(u, t, integrator) = subnorm(u - x_f; directions=rad_dims) < rad_f affect!(integrator) = terminate!(integrator) cb_ball = DiscreteCallback(condition, affect!) - + prob = remake(sys.integ.sol.prob; u0=x_i, tspan=(0, tmax)) sim = solve(prob, sys.integ.alg; callback=cb_ball, kwargs...) success = sim.retcode == SciMLBase.ReturnCode.Terminated diff --git a/src/utils.jl b/src/utils.jl index e62dd15b..502bc922 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -113,7 +113,7 @@ Normalizes the covariance matrix ``Q`` (in-place) by dividing it by function normalize_covariance!(covariance) l1norm = norm(covariance, 1) dim = size(covariance)[1] - return covariance * dim/l1norm + return covariance * dim / l1norm end # Central finite difference, second derivative diff --git a/test/largedeviations/action_fhn.jl b/test/largedeviations/action_fhn.jl index 7f0682ec..65f90791 100644 --- a/test/largedeviations/action_fhn.jl +++ b/test/largedeviations/action_fhn.jl @@ -48,4 +48,4 @@ end @testset "div_drift" begin @test CT.div_drift(sys, zeros(2)) == -2.0 @test CT.div_drift(sys, x_i) == -4.0 -end \ No newline at end of file +end diff --git a/test/runtests.jl b/test/runtests.jl index a6b4e389..42e7572f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -64,7 +64,7 @@ end @testset "Trajectories" begin include("trajectories/simulate.jl") -# include("trajectories/transition.jl") + # include("trajectories/transition.jl") end @testset "Extentions" begin diff --git a/test/trajectories/simulate.jl b/test/trajectories/simulate.jl index ac486ce4..b5daa156 100644 --- a/test/trajectories/simulate.jl +++ b/test/trajectories/simulate.jl @@ -6,8 +6,8 @@ sys = CoupledSDEs(fitzhugh_nagumo, init, p; noise_strength=σ) traj = trajectory(sys, 10, init) relax = deterministic_orbit(sys, 10, init) - - @test traj[1][1,1] == 1.0 - @test isapprox(relax[1][end,1], 0.816; atol=1e-2) + + @test traj[1][1, 1] == 1.0 + @test isapprox(relax[1][end, 1], 0.816; atol=1e-2) # These tests could be improved - Reyk end