diff --git a/README.md b/README.md index 64c6d9d5..815395ca 100644 --- a/README.md +++ b/README.md @@ -14,12 +14,22 @@ See [package documentation](https://juliadynamics.github.io/CriticalTransitions. ```julia using CriticalTransitions +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 + # System parameters p = [1., 3., 1., 1., 1., 0.] noise_strength = 0.02 # Define stochastic system -sys = CoupledSDEs(meier_stein, diagonal_noise(σ), zeros(2), p) +sys = CoupledSDEs(fitzhugh_nagumo, id_func, zeros(2), noise_strength, p) # Get stable fixed points fps, eigs, stab = fixedpoints(sys, [-2,-2], [2,2]) diff --git a/docs/make.jl b/docs/make.jl index 26720205..c9c227a1 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -12,9 +12,9 @@ authors = join(project_toml["authors"], ", ") * " and contributors" github = "https://github.com/juliadynamics/CriticalTransitions.jl" links = InterLinks( - "DiffEqNoiseProcess" => "https://docs.sciml.ai/DiffEqNoiseProcess/stable/", - "DifferentialEquations" => "https://docs.sciml.ai/DiffEqDocs/stable/", - "StochasticDiffEq" => "https://docs.sciml.ai/DiffEqDocs/stable/", +# "DiffEqNoiseProcess" => "https://docs.sciml.ai/DiffEqNoiseProcess/stable/", +# "DifferentialEquations" => "https://docs.sciml.ai/DiffEqDocs/stable/", +# "StochasticDiffEq" => "https://docs.sciml.ai/DiffEqDocs/stable/", ) bib = CitationBibliography(joinpath(@__DIR__, "src", "refs.bib"); style=:numeric) diff --git a/docs/src/man/CoupledSDEs.md b/docs/src/man/CoupledSDEs.md index 6e17e4b1..c6a2b82c 100644 --- a/docs/src/man/CoupledSDEs.md +++ b/docs/src/man/CoupledSDEs.md @@ -44,7 +44,12 @@ or equivalently ```@example type sde = CoupledSDEs(f!, idfunc!, zeros(2), σ) ``` -The `diagonal_noise!` function is a helper function equivalent to `(du, u, p, t) -> σ .* idfunc!(du, u, p, t)` with `idfunc!(du, u, p, t) = (du .= ones(length(u)); return nothing)`. 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). +where we used the helper function +```@docs +idfunc! +idfunc +``` +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). ```@example type sol = simulate(sde, 1.0, dt=0.01, alg=SOSRA()) plot(sol) diff --git a/docs/src/man/largedeviations.md b/docs/src/man/largedeviations.md index 937ca127..543c69fc 100644 --- a/docs/src/man/largedeviations.md +++ b/docs/src/man/largedeviations.md @@ -4,23 +4,23 @@ ### Freidlin-Wentzell action ```@docs -fw_action(sys::CoupledSDEs, path, time; kwargs...) +fw_action ``` ### Geometric Freidlin-Wentzell action ```@docs -geometric_action(sys::CoupledSDEs, path, arclength=1; kwargs...) +geometric_action ``` ### Onsager-Machlup action ```@docs -om_action(sys::CoupledSDEs, path, time; kwargs...) +om_action ``` For convenience, a general [`action`](@ref) function is available where the type of functional is set as an argument: ```@docs -action(sys::CoupledSDEs, path::Matrix, time, functional; kwargs...) +action ``` ## Minimum action paths @@ -31,7 +31,7 @@ between two states of a `CoupledSDEs`. Minimization of the Freidlin-Wentzell action using the L-BFGS algorithm of `Optim.jl`. ```@docs -min_action_method(sys::CoupledSDEs, x_i, x_f, N::Int, T::Float64; kwargs...) +min_action_method ``` ### Geometric minimum action method (gMAM) @@ -39,5 +39,5 @@ Minimization of the geometric action following [Heymann and Vanden-Eijnden, PRL (2008)](https://link.aps.org/doi/10.1103/PhysRevLett.100.140601). ```@docs -geometric_min_action_method(sys::CoupledSDEs, x_i, x_f, arclength=1; kwargs...) +geometric_min_action_method ``` \ No newline at end of file diff --git a/docs/src/man/systemanalysis.md b/docs/src/man/systemanalysis.md index 0d0548e4..369c6d34 100644 --- a/docs/src/man/systemanalysis.md +++ b/docs/src/man/systemanalysis.md @@ -10,9 +10,9 @@ fixedpoints ## Basins of attraction ```@docs -basins(sys::CoupledSDEs, A, B, C, H; kwargs) -basinboundary(X, Y, h; kwargs...) -basboundary(sys::CoupledSDEs, xrange::Vector, yrange::Vector, xspacing::Float64, attractors::Vector; kwargs...) +basins +basinboundary +basboundary ``` ## Edge tracking diff --git a/docs/src/man/utils.md b/docs/src/man/utils.md index 372408ec..b0bb5bf0 100644 --- a/docs/src/man/utils.md +++ b/docs/src/man/utils.md @@ -14,8 +14,20 @@ Lyapunov spectrum of a `CoupledSDEs`, here exemplified by the FitzHugh-Nagumo mo computed by typing: ```julia -using CriticalTransitions, DynamicalSystems: lyapunovspectrum -sys = CoupledSDEs(fitzhugh_nagumo, diagonal_noise(σ), zeros(2), [1.,3.,1.,1.,1.,0.]) +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) ``` @@ -26,7 +38,11 @@ CoupledODEs(sys::CoupledSDEs; diffeq, t0=0.0) ## `CoupledSDEs` utility functions ```@docs +noise_strength +covariance_matrix +noise_process CriticalTransitions.drift +CriticalTransitions.div_drift ``` ## Saving data diff --git a/src/CoupledSDEs.jl b/src/CoupledSDEs.jl index 357d84b8..ab29fc7c 100644 --- a/src/CoupledSDEs.jl +++ b/src/CoupledSDEs.jl @@ -36,7 +36,7 @@ coupled ordinary differential equations as follows: 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 stochastic function `g` is multiplied by `σ`. +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). diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index fe81ebf2..ea40e366 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -14,11 +14,6 @@ using Symbolics using Optim, Dierckx using Printf, DrWatson, Dates, Statistics, Markdown -# Define custom types -# Parameters = Union{Vector{Any}, Nothing}; -# CovMatrix = Union{Matrix, UniformScaling{Bool}, Diagonal{Bool, Vector{Bool}}}; -# State = Union{Vector, SVector} - include("extention_functions.jl") include("utils.jl") include("CoupledSDEs.jl") @@ -41,7 +36,7 @@ export CoupledSDEs, export equilib, fixedpoints, basins, basinboundary, basboundary export simulate, relax export transition, transitions -export fw_integrand, fw_action, om_action, action, geometric_action +export fw_action, om_action, action, geometric_action export min_action_method, geometric_min_action_method export basins, basinboundary export edgetracking, bisect_to_edge, attractor_mapper diff --git a/src/largedeviations/action.jl b/src/largedeviations/action.jl index 79ec27d7..fc6c7dbd 100644 --- a/src/largedeviations/action.jl +++ b/src/largedeviations/action.jl @@ -1,4 +1,4 @@ -@doc doc""" +""" $(TYPEDSIGNATURES) Calculates the Freidlin-Wentzell action of a given `path` with time points `time` in a @@ -32,7 +32,7 @@ function fw_action(sys::CoupledSDEs, path, time) return S / 2 end; -@doc doc""" +""" $(TYPEDSIGNATURES) Calculates the Onsager-Machlup action of a given `path` with time points `time` for the drift field `sys.f` at given `sys.noise_strength`. @@ -66,7 +66,7 @@ function om_action(sys::CoupledSDEs, path, time) return fw_action(sys, path, time) + S / 2 end; -@doc doc""" +""" $(TYPEDSIGNATURES) Computes the action functional specified by `functional` for a given CoupledSDEs `sys` and @@ -84,7 +84,7 @@ function action(sys::CoupledSDEs, path::Matrix, time, functional; kwargs...) return action end; -@doc doc""" +""" $(TYPEDSIGNATURES) Calculates the geometric action of a given `path` with specified `arclength` for the drift @@ -124,12 +124,11 @@ function geometric_action(sys::CoupledSDEs, path, arclength=1.0) return S * arclength / (N - 1) end -@doc doc""" +""" $(TYPEDSIGNATURES) Computes the squared ``A``-norm ``|| \\dot \\phi_t - b ||^2_A`` (see `fw_action` for -details). Returns a vector of length `N` containing the values of the above squared norm for -each time point in the vector `time`. +details). Returns a vector of length `N` containing the values of the above squared norm for each time point in the vector `time`. """ function fw_integrand(sys::CoupledSDEs, path, time, A) v = path_velocity(path, time; order=4) @@ -143,7 +142,7 @@ function fw_integrand(sys::CoupledSDEs, path, time, A) return sqnorm end; -@doc doc""" +""" $(TYPEDSIGNATURES) Returns the velocity along a given `path` with time points given by `time`. diff --git a/src/largedeviations/min_action_method.jl b/src/largedeviations/min_action_method.jl index 33a684b8..20e39a63 100644 --- a/src/largedeviations/min_action_method.jl +++ b/src/largedeviations/min_action_method.jl @@ -73,7 +73,7 @@ initial condition `init`, given a system `sys` and total path time `T`. The initial path `init` must be a matrix of size `(D, N)`, where `D` is the dimension of the system and `N` is the number of path points. The physical time of the path is specified by `T`, such that the time step between consecutive path points is -``\\Deltat = T/N``. +``\\Delta t = T/N``. For more information see the main method, [`min_action_method(sys::CoupledSDEs, x_i, x_f, N::Int, T::Real; kwargs...)`](@ref). diff --git a/src/utils.jl b/src/utils.jl index 4970fdcc..2142f6b2 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -6,33 +6,24 @@ Functions in this file are independent of types/functions defined in CriticalTra """ $(TYPEDSIGNATURES) -Identity function for noise function of `CoupledSDEs` (out-of-place). +Identity function for a diffusion function `g` of `CoupledSDEs` (out-of-place). +Equivalent to `(u, p, t) -> ones(length(u))`, """ -function idfunc(u, p, t) +function idfunc(u, p, t) # TODO: return same type as u return ones(length(u)) end; """ $(TYPEDSIGNATURES) -Identity function for noise function `CoupledSDEs` (in-place). +Identity function for a diffusion function `g` of `CoupledSDEs` (in-place). +Equivalent to `idfunc!(du, u, p, t) = (du .= ones(length(u)); return nothing)` """ function idfunc!(du, u, p, t) du .= ones(length(u)) return nothing end; -# """ -# """ -# 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) - function σg(σ, g) return (u, p, t) -> σ .* g(u, p, t) end @@ -74,50 +65,7 @@ function intervals_to_box(bmin::Vector, bmax::Vector) return box end; -# function additive_idx!(du, u, p, t, idx) -# du[idx] = 1.0 -# return nothing -# end; - -# function additive_idx(u, p, t, idx) -# du = zeros(length(u)) -# du[idx] = 1.0 -# return SVector{length(u)}(du) -# end; - -# function multiplicative_idx!(du, u, p, t, idx) -# return du[idx] = u[idx] -# end; - -# function multiplicative_idx(u, p, t, idx) -# du = zeros(length(u)) -# du[idx] = u[idx] -# return SVector{length(u)}(du) -# end - -# function multiplicative_first!(du, u, p, t) -# return du[1] = u[1] -# end; - -# function multiplicative_first(u, p, t) -# du = zeros(length(u)) -# du[1] = u[1] - -# return SVector{length(u)}(du) -# end; - -# function additive_first!(du, u, p, t) -# return du[1] = 1 -# end; - -# function additive_first(u, p, t) -# du = zeros(length(u)) -# du[1] = 1 - -# return SVector{length(u)}(du) -# end; - -@doc doc""" +""" Calculates the generalized ``A``-norm of the vector `vec`, ``||v||_A := \\sqrt(v^\\top \\cdot A \\cdot v)``, where `A` is a square matrix of dimension `(length(vec) x length(vec))`. @@ -130,7 +78,7 @@ function anorm(vec, A; square=false) return square ? normsquared : sqrt(normsquared) end; -@doc doc""" +""" $(TYPEDSIGNATURES) Returns the Euclidean norm of the vector `vec`; however, if `directions` are specified, the @@ -164,7 +112,7 @@ function central2(f, idx, dz) return (f[idx + 1] - 2f[idx] + f[idx - 1]) / (dz^2) end; -@doc doc""" +""" $(TYPEDSIGNATURES) Smooth approximation of `abs(x)`, ``|x| = x \\tanh(\\xi x)``, where ``xi`` controls the accuracy of the approximation. The exact absolute value function is obtained in the limit diff --git a/systems/CTLibrary.jl b/systems/CTLibrary.jl index 3d86c0bf..2c10145b 100644 --- a/systems/CTLibrary.jl +++ b/systems/CTLibrary.jl @@ -10,6 +10,4 @@ include("rooth.jl") include("stommel.jl") include("rivals.jl") -export fitzhugh_nagumo - end diff --git a/test/largedeviations/action_fhn.jl b/test/largedeviations/action_fhn.jl index e4b6f58f..7e1d963e 100644 --- a/test/largedeviations/action_fhn.jl +++ b/test/largedeviations/action_fhn.jl @@ -39,7 +39,7 @@ end # Test fw_integrand function @testset "fw_integrand" begin - integrand = fw_integrand(sys, path, time, A) + integrand = CriticalTransitions.fw_integrand(sys, path, time, A) @test minimum(integrand) > 0.18 end