From e58aeabbf1decff109576523d2423a3329d40062 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 8 Oct 2024 18:18:10 +0100 Subject: [PATCH 01/76] added first draft for RateSystem --- src/RateSys.jl | 53 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) create mode 100644 src/RateSys.jl diff --git a/src/RateSys.jl b/src/RateSys.jl new file mode 100644 index 00000000..8414e723 --- /dev/null +++ b/src/RateSys.jl @@ -0,0 +1,53 @@ +# | | | | +# t_i autonomous t_ni non-autonomous t_nf autonomous t_f +# | | | | + +using DynamicalSystems +using Plots + +# we write + +function RateSyst(tni,tnf,f,λ,p_λ,initvals) + func(u,p,t) = combined_system(u,t,tni,tnf,f,λ,p_λ); + return ContinuousDynamicalSystem(func, initvals, Float64[], t0=tni) +end + +function combined_system(u,t,tni,tnf,f,λ,p_λ) + lambda = t < tni ? λ(u,p_λ,tni) : tni <= t <= tnf ? λ(u,p_λ,t) : λ(u,p_λ,tnf) + return f(u,lambda,t) +end; + + +############### +# user writes + +function f(u::SVector{1, Float64},p::SVector{1, Float64},t::Float64) + x = u[1] + λ = p[1] + dx = (x+λ)^2 - 1 + return SVector{1}(dx) +end; + +function λ(p::Vector{Float64},t::Float64) + r,λ_max = p + lambda = (λ_max/2)*(tanh(λ_max*r*t/2) +1) + return SVector{1}(lambda) +end; + +tni=-10.; tnf=10.; p_λ = [1.0,3.0]; initvals = [-4.]; +sys = RateSyst(tni,tnf,f,λ,p_λ,initvals) + +traj=trajectory(sys,40.,t0=-20.) + +trajic=[traj[1][i][1] for i in 1:length(traj[1])] +plot(collect(traj[2]),trajic) + + +# further ideas +# function futureSyst(RateSyst) +# Canard Trajectories + +# \dot(x) = f(x(t),λ(t)) + + + From 525c6a1d73f4fae849f03609c9be494b9ea2ed9b Mon Sep 17 00:00:00 2001 From: reykboerner Date: Wed, 16 Oct 2024 18:17:29 +0100 Subject: [PATCH 02/76] incorporated RateSystem file --- src/CriticalTransitions.jl | 3 +++ src/{RateSys.jl => RateSystem.jl} | 26 ++------------------------ test/RateSystem.jl | 22 ++++++++++++++++++++++ 3 files changed, 27 insertions(+), 24 deletions(-) rename src/{RateSys.jl => RateSystem.jl} (52%) create mode 100644 test/RateSystem.jl diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index cab1d333..dffb3562 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -21,6 +21,7 @@ using StochasticDiffEq: u_modified! using DynamicalSystemsBase: DynamicalSystemsBase, + ContinuousTimeDynamicalSystem, CoupledSDEs, CoupledODEs, dynamic_rule, @@ -55,6 +56,7 @@ include("extention_functions.jl") include("utils.jl") include("system_utils.jl") include("io.jl") +include("RateSystem.jl") include("trajectories/simulation.jl") include("trajectories/transition.jl") include("trajectories/equib.jl") @@ -69,6 +71,7 @@ using .CTLibrary # Core types export CoupledSDEs, CoupledODEs, noise_process, covariance_matrix, diffusion_matrix export dynamic_rule, current_state, set_state!, trajectory +export RateSystem # Methods export drift, div_drift diff --git a/src/RateSys.jl b/src/RateSystem.jl similarity index 52% rename from src/RateSys.jl rename to src/RateSystem.jl index 8414e723..dd6c0de0 100644 --- a/src/RateSys.jl +++ b/src/RateSystem.jl @@ -2,12 +2,9 @@ # t_i autonomous t_ni non-autonomous t_nf autonomous t_f # | | | | -using DynamicalSystems -using Plots - # we write -function RateSyst(tni,tnf,f,λ,p_λ,initvals) +function RateSystem(tni,tnf,f,λ,p_λ,initvals) func(u,p,t) = combined_system(u,t,tni,tnf,f,λ,p_λ); return ContinuousDynamicalSystem(func, initvals, Float64[], t0=tni) end @@ -21,26 +18,7 @@ end; ############### # user writes -function f(u::SVector{1, Float64},p::SVector{1, Float64},t::Float64) - x = u[1] - λ = p[1] - dx = (x+λ)^2 - 1 - return SVector{1}(dx) -end; - -function λ(p::Vector{Float64},t::Float64) - r,λ_max = p - lambda = (λ_max/2)*(tanh(λ_max*r*t/2) +1) - return SVector{1}(lambda) -end; - -tni=-10.; tnf=10.; p_λ = [1.0,3.0]; initvals = [-4.]; -sys = RateSyst(tni,tnf,f,λ,p_λ,initvals) - -traj=trajectory(sys,40.,t0=-20.) - -trajic=[traj[1][i][1] for i in 1:length(traj[1])] -plot(collect(traj[2]),trajic) +# ...moved to test/RateSystem.jl (Reyk) # further ideas diff --git a/test/RateSystem.jl b/test/RateSystem.jl new file mode 100644 index 00000000..768b7794 --- /dev/null +++ b/test/RateSystem.jl @@ -0,0 +1,22 @@ +using Plots + +function f(u::SVector{1, Float64},p::SVector{1, Float64},t::Float64) + x = u[1] + λ = p[1] + dx = (x+λ)^2 - 1 + return SVector{1}(dx) +end; + +function λ(p::Vector{Float64},t::Float64) + r,λ_max = p + lambda = (λ_max/2)*(tanh(λ_max*r*t/2) +1) + return SVector{1}(lambda) +end; + +tni=-10.; tnf=10.; p_λ = [1.0,3.0]; initvals = [-4.]; +sys = RateSystem(tni,tnf,f,λ,p_λ,initvals) + +traj=trajectory(sys,40.,t0=-20.) + +trajic=[traj[1][i][1] for i in 1:length(traj[1])] +plot(collect(traj[2]),trajic) \ No newline at end of file From 4e8366c2578017dd70d552f234cfc53023727a8f Mon Sep 17 00:00:00 2001 From: reykboerner Date: Wed, 16 Oct 2024 18:25:50 +0100 Subject: [PATCH 03/76] small typo fix --- src/RateSystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/RateSystem.jl b/src/RateSystem.jl index dd6c0de0..8046b9b0 100644 --- a/src/RateSystem.jl +++ b/src/RateSystem.jl @@ -6,7 +6,7 @@ function RateSystem(tni,tnf,f,λ,p_λ,initvals) func(u,p,t) = combined_system(u,t,tni,tnf,f,λ,p_λ); - return ContinuousDynamicalSystem(func, initvals, Float64[], t0=tni) + return ContinuousTimeDynamicalSystem(func, initvals, Float64[], t0=tni) end function combined_system(u,t,tni,tnf,f,λ,p_λ) From 727740fa962f46851321b8d2bff90ac4e43ae5da Mon Sep 17 00:00:00 2001 From: reykboerner Date: Wed, 16 Oct 2024 18:29:23 +0100 Subject: [PATCH 04/76] undid wrong fix --- src/CriticalTransitions.jl | 1 - src/RateSystem.jl | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index dffb3562..526af9be 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -21,7 +21,6 @@ using StochasticDiffEq: u_modified! using DynamicalSystemsBase: DynamicalSystemsBase, - ContinuousTimeDynamicalSystem, CoupledSDEs, CoupledODEs, dynamic_rule, diff --git a/src/RateSystem.jl b/src/RateSystem.jl index 8046b9b0..b2dd11ed 100644 --- a/src/RateSystem.jl +++ b/src/RateSystem.jl @@ -6,7 +6,7 @@ function RateSystem(tni,tnf,f,λ,p_λ,initvals) func(u,p,t) = combined_system(u,t,tni,tnf,f,λ,p_λ); - return ContinuousTimeDynamicalSystem(func, initvals, Float64[], t0=tni) + return CoupledODEs(func, initvals, Float64[], t0=tni) end function combined_system(u,t,tni,tnf,f,λ,p_λ) From 3acaddbabaa319734f41473df09ae6d64fc12a28 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Wed, 19 Mar 2025 18:41:06 +0000 Subject: [PATCH 05/76] Removed RateSystem Lines --- src/CriticalTransitions.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 526af9be..cab1d333 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -55,7 +55,6 @@ include("extention_functions.jl") include("utils.jl") include("system_utils.jl") include("io.jl") -include("RateSystem.jl") include("trajectories/simulation.jl") include("trajectories/transition.jl") include("trajectories/equib.jl") @@ -70,7 +69,6 @@ using .CTLibrary # Core types export CoupledSDEs, CoupledODEs, noise_process, covariance_matrix, diffusion_matrix export dynamic_rule, current_state, set_state!, trajectory -export RateSystem # Methods export drift, div_drift From d32c1c563c748bc7edfacfdc49d26876ec87d3dd Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Wed, 19 Mar 2025 19:09:12 +0000 Subject: [PATCH 06/76] Return RateSystem changes --- src/CriticalTransitions.jl | 106 +++++++++++++++++-------------------- 1 file changed, 48 insertions(+), 58 deletions(-) diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 4c9f9e12..526af9be 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -2,15 +2,23 @@ module CriticalTransitions # Base using Statistics: Statistics, mean -using LinearAlgebra: LinearAlgebra, I, norm, dot, tr, det +using LinearAlgebra: LinearAlgebra, I, norm, dot, tr using StaticArrays: StaticArrays, SVector -using SparseArrays: spdiagm -using DataStructures: CircularBuffer -using Random: Random # Core -using SciMLBase: EnsembleThreads, DiscreteCallback, remake, terminate! -using StochasticDiffEq: StochasticDiffEq +using DiffEqNoiseProcess: DiffEqNoiseProcess +using OrdinaryDiffEq: OrdinaryDiffEq, Tsit5 +using StochasticDiffEq: + StochasticDiffEq, + DiscreteCallback, + ODEProblem, + SDEFunction, + SOSRA, + remake, + solve, + step!, + terminate!, + u_modified! using DynamicalSystemsBase: DynamicalSystemsBase, CoupledSDEs, @@ -18,80 +26,62 @@ using DynamicalSystemsBase: dynamic_rule, current_state, set_state!, - trajectory, - jacobian, - StateSpaceSet + trajectory -using Interpolations: linear_interpolation -using Optimization -using OptimizationOptimisers: Optimisers -using LinearSolve: LinearProblem, KLUFactorization, solve +using ForwardDiff: ForwardDiff +using IntervalArithmetic: IntervalArithmetic, interval +using Dierckx: Dierckx, ParametricSpline +using Optim: Optim, LBFGS +using Symbolics: Symbolics # io and documentation using Format: Format +using Dates: Dates using Printf: Printf -using DocStringExtensions: TYPEDSIGNATURES, TYPEDEF, TYPEDFIELDS, METHODLIST -using ProgressMeter: Progress, next! +using Markdown: Markdown +using DocStringExtensions: TYPEDSIGNATURES +using HDF5: HDF5, h5open, push! +using JLD2: JLD2, jldopen +using ProgressBars: ProgressBars, tqdm +using ProgressMeter: ProgressMeter # reexport using Reexport: @reexport @reexport using StaticArrays -@reexport using DynamicalSystemsBase +@reexport using StochasticDiffEq +@reexport using DiffEqNoiseProcess -include("extension_functions.jl") +include("extention_functions.jl") include("utils.jl") -include("sde_utils.jl") - -include("trajectories/TransitionEnsemble.jl") +include("system_utils.jl") +include("io.jl") +include("RateSystem.jl") include("trajectories/simulation.jl") include("trajectories/transition.jl") - -include("transition_path_theory/TransitionPathMesh.jl") -include("transition_path_theory/langevin.jl") -include("transition_path_theory/committor.jl") -include("transition_path_theory/invariant_pdf.jl") -include("transition_path_theory/reactive_current.jl") -include("transition_path_theory/probability.jl") - -include("largedeviations/utils.jl") +include("trajectories/equib.jl") +include("noiseprocesses/stochprocess.jl") include("largedeviations/action.jl") -include("largedeviations/MinimumActionPath.jl") include("largedeviations/min_action_method.jl") include("largedeviations/geometric_min_action_method.jl") -include("largedeviations/sgMAM.jl") -include("largedeviations/string_method.jl") - include("../systems/CTLibrary.jl") using .CTLibrary # Core types export CoupledSDEs, CoupledODEs, noise_process, covariance_matrix, diffusion_matrix -export drift, div_drift, solver -export StateSpaceSet - -export sgmam, SgmamSystem -export fw_action, om_action, action, geometric_action -export min_action_method, geometric_min_action_method, string_method -export MinimumActionPath +export dynamic_rule, current_state, set_state!, trajectory +export RateSystem -export deterministic_orbit +# Methods +export drift, div_drift +export equilib, deterministic_orbit export transition, transitions - -export distmesh2D, dellipse, ddiff -export TransitionPathMesh, Committor -export get_ellipse, reparametrization -export find_boundary, huniform, dunion - -export committor, - invariant_pdf, reactive_current, probability_reactive, probability_last_A, Langevin - -# Error hint for extensions stubs -function __init__() - Base.Experimental.register_error_hint( - _basin_error_hinter(intervals_to_box), MethodError - ) - return nothing -end - +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 covariance_matrix, diffusion_matrix +# export edgetracking, bisect_to_edge, AttractorsViaProximity +# export fixedpoints +# ^ extention tests needed end # module CriticalTransitions From 5ee2ed9dafd7be5d479b4a7227c6500c2bb1621b Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Wed, 19 Mar 2025 19:17:00 +0000 Subject: [PATCH 07/76] Again delete Rate stuff --- src/CriticalTransitions.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 526af9be..cab1d333 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -55,7 +55,6 @@ include("extention_functions.jl") include("utils.jl") include("system_utils.jl") include("io.jl") -include("RateSystem.jl") include("trajectories/simulation.jl") include("trajectories/transition.jl") include("trajectories/equib.jl") @@ -70,7 +69,6 @@ using .CTLibrary # Core types export CoupledSDEs, CoupledODEs, noise_process, covariance_matrix, diffusion_matrix export dynamic_rule, current_state, set_state!, trajectory -export RateSystem # Methods export drift, div_drift From 444eb72847582331bbb7b281eefda9cf73f1e446 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Wed, 19 Mar 2025 19:25:51 +0000 Subject: [PATCH 08/76] Include Rate stuff again --- src/CriticalTransitions.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index cab1d333..fcec3409 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -62,6 +62,7 @@ include("noiseprocesses/stochprocess.jl") include("largedeviations/action.jl") include("largedeviations/min_action_method.jl") include("largedeviations/geometric_min_action_method.jl") +include("RateSystem.jl") include("../systems/CTLibrary.jl") using .CTLibrary @@ -79,6 +80,7 @@ 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 covariance_matrix, diffusion_matrix +export RateSystem # export edgetracking, bisect_to_edge, AttractorsViaProximity # export fixedpoints # ^ extention tests needed From b8f145a70e40ea17bcb29906afd42bb889f39103 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Wed, 19 Mar 2025 19:29:18 +0000 Subject: [PATCH 09/76] added new RateSystem Version --- src/RateSystemDraft2.jl | 115 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 115 insertions(+) create mode 100644 src/RateSystemDraft2.jl diff --git a/src/RateSystemDraft2.jl b/src/RateSystemDraft2.jl new file mode 100644 index 00000000..51694e5d --- /dev/null +++ b/src/RateSystemDraft2.jl @@ -0,0 +1,115 @@ +using GLMakie +using CriticalTransitions +using DifferentialEquations + +# we consider the ODE dxₜ/dt = f(xₜ,λ(rt)) +# here λ = λ(t) ∈ Rᵐ is a function containing all the system parameters + +# We ask the user to define: +# 1) a system that should be investigated and +# 2) a protocol for λ(rt) by calling the function rate_protocol + + +############################################################################################################## +# we write behind the scenes +############################################################################################################## + +mutable struct RateProtocol + λ::Function + p_lambda::Vector + r::Float64 + t_start::Float64 + t_end::Float64 +end + +RateProtocol(λ::Function,p_lambda::Vector,r::Float64)=RateProtocol(λ,p_lambda,r,-Inf,Inf) +RateProtocol(λ::Function,r::Float64)=RateProtocol(λ,[],r,-Inf,Inf) +#RateProtocol(λ::Function,p_lambda::Vector,r::Float64,t_start::Float64)=RateProtocol(λ,p_lambda,r,t_start,Inf) +#RateProtocol(λ::Function,r::Float64,t_start::Float64)=RateProtocol(λ,[],r,t_start,Inf) + + +function modified_drift(u,p,t,ds::ContinuousTimeDynamicalSystem,λ::Function,t_start::Float64,t_end::Float64,r::Float64; + kwargs...) + + if t_start > t_end + error("Please ensure that t_start ≤ t_end.") + end + + p̃ = r*t ≤ t_start ? λ(p,t_start;kwargs...) : t_start < r*t < t_end ? λ(p,r*t;kwargs...) : λ(p,t_end;kwargs...); # the value(s) of λ(rt) + return ds.integ.f(u,p̃,t) +end; + +function RateSystem(auto_sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t0::Float64; + kwargs...) + # we wish to return a continuous time dynamical system with modified drift field + + f(u,p,t) = modified_drift(u,p,t,auto_sys,rp.λ,rp.t_start,rp.t_end,rp.r;kwargs...); + prob = remake(auto_sys.integ.sol.prob;f,p=rp.p_lambda,tspan=(t0,Inf)); + nonauto_sys = CoupledODEs(prob,auto_sys.diffeq); + return nonauto_sys +end; + +############################################################################################################## +# user writes (pseudo example - see the bottom of the page for an explicit example) +############################################################################################################## + + +############################################################################################################## +# test (prototypical model used for studying R-tipping, with critical rate r = 4/3) +############################################################################################################## + + +# autonomous system + +function f(u,p,t) + x = u[1] + λ = p[1] + dx = (x+λ)^2 - 1 + return SVector{1}(dx) +end; + +lambda = 0.; +p = [lambda]; + +z₀ = [-1.]; + +Δt = 1.e-3; +diffeq = (alg = AutoVern9(Rodas5(autodiff=false)),abstol=1.e-16,reltol=1.e-16); + +auto_sys = CoupledODEs(f,z₀,p;diffeq) + +# time-dependent parameter function + +function λ(p,t) + λ_max = p[1] + lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) + return SVector{1}(lambda) +end; + +# calling RateSystem + +r = 4/3-0.02; +λ_max = 3.; p_lambda = [λ_max]; + +#t_start = -Inf; t_end = Inf; +t0 = -10.; + +## test v2 + +rp = RateProtocol(λ,p_lambda,r); + +nonauto_sys = RateSystem(auto_sys,rp,t0); + +#nonauto_sys = RateSystem(auto_sys,λ,r,t0;p_lambda,t_start,t_end); + +T = 20.; + +auto_traj = trajectory(auto_sys,T,z₀); +nonauto_traj = trajectory(nonauto_sys,T,z₀); + +fig = Figure(); axs = Axis(fig[1,1]); +lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}"); +lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}"); +axislegend(axs,position=:rc,labelsize=40); +ylims!(-4.5,1.5); +fig \ No newline at end of file From f8230d773b13f9a17d62350edf053bab7d6871a4 Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Wed, 19 Mar 2025 21:09:37 +0100 Subject: [PATCH 10/76] Compressed RateSystemDraft2.jl script, created RateSystemTestsDraft.jl script, created a dev folder, updated src/CriticalTransitions.jl file. --- dev/RateSystemDraft.jl | 42 +++++++++++++ dev/RateSystemTestsDraft.jl | 58 ++++++++++++++++++ src/CriticalTransitions.jl | 6 +- src/RateSystemDraft2.jl | 115 ------------------------------------ 4 files changed, 104 insertions(+), 117 deletions(-) create mode 100644 dev/RateSystemDraft.jl create mode 100644 dev/RateSystemTestsDraft.jl delete mode 100644 src/RateSystemDraft2.jl diff --git a/dev/RateSystemDraft.jl b/dev/RateSystemDraft.jl new file mode 100644 index 00000000..41c8ee44 --- /dev/null +++ b/dev/RateSystemDraft.jl @@ -0,0 +1,42 @@ +# we consider the ODE dxₜ/dt = f(xₜ,λ(rt)) +# here λ = λ(t) ∈ Rᵐ is a function containing all the system parameters + +# We ask the user to define: +# 1) a system that should be investigated and +# 2) a protocol for the time-dependent forcing with the struct RateProtocol + +mutable struct RateProtocol + λ::Function + p_lambda::Vector + r::Float64 + t_start::Float64 + t_end::Float64 +end + +# convenience functions + +RateProtocol(λ::Function,p_lambda::Vector,r::Float64)=RateProtocol(λ,p_lambda,r,-Inf,Inf) +RateProtocol(λ::Function,r::Float64)=RateProtocol(λ,[],r,-Inf,Inf) +#RateProtocol(λ::Function,p_lambda::Vector,r::Float64,t_start::Float64)=RateProtocol(λ,p_lambda,r,t_start,Inf) +#RateProtocol(λ::Function,r::Float64,t_start::Float64)=RateProtocol(λ,[],r,t_start,Inf) + +function modified_drift(u,p,t,ds::ContinuousTimeDynamicalSystem,λ::Function,t_start::Float64,t_end::Float64,r::Float64; + kwargs...) + + if t_start > t_end + error("Please ensure that t_start ≤ t_end.") + end + + p̃ = r*t ≤ t_start ? λ(p,t_start;kwargs...) : t_start < r*t < t_end ? λ(p,r*t;kwargs...) : λ(p,t_end;kwargs...); # the value(s) of λ(rt) + return ds.integ.f(u,p̃,t) +end; + +function RateSystem(auto_sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t0::Float64; + kwargs...) + # we wish to return a continuous time dynamical system with modified drift field + + f(u,p,t) = modified_drift(u,p,t,auto_sys,rp.λ,rp.t_start,rp.t_end,rp.r;kwargs...); + prob = remake(auto_sys.integ.sol.prob;f,p=rp.p_lambda,tspan=(t0,Inf)); + nonauto_sys = CoupledODEs(prob,auto_sys.diffeq); + return nonauto_sys +end; diff --git a/dev/RateSystemTestsDraft.jl b/dev/RateSystemTestsDraft.jl new file mode 100644 index 00000000..2decbdf4 --- /dev/null +++ b/dev/RateSystemTestsDraft.jl @@ -0,0 +1,58 @@ +############################################################################################################## +# test (prototypical model used for studying R-tipping, with critical rate r = 4/3) +############################################################################################################## + +# autonomous system + +function f(u,p,t) + x = u[1] + λ = p[1] + dx = (x+λ)^2 - 1 + return SVector{1}(dx) +end; + +lambda = 0.; +p = [lambda]; + +z₀ = [-1.]; + +Δt = 1.e-3; +diffeq = (alg = AutoVern9(Rodas5(autodiff=false)),abstol=1.e-16,reltol=1.e-16); + +auto_sys = CoupledODEs(f,z₀,p;diffeq) + +# time-dependent parameter function + +function λ(p,t) + λ_max = p[1] + lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) + return SVector{1}(lambda) +end; + +# calling RateSystem + +r = 4/3-0.02; +λ_max = 3.; p_lambda = [λ_max]; + +#t_start = -Inf; t_end = Inf; +t0 = -10.; + +## test v2 + +rp = RateProtocol(λ,p_lambda,r); + +nonauto_sys = RateSystem(auto_sys,rp,t0); + +#nonauto_sys = RateSystem(auto_sys,λ,r,t0;p_lambda,t_start,t_end); + +T = 20.; + +auto_traj = trajectory(auto_sys,T,z₀); +nonauto_traj = trajectory(nonauto_sys,T,z₀); + +# fig = Figure(); axs = Axis(fig[1,1]); +# lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}"); +# lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}"); +# axislegend(axs,position=:rc,labelsize=40); +# ylims!(-4.5,1.5); +# fig \ No newline at end of file diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index fcec3409..53d2c35e 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -62,7 +62,8 @@ include("noiseprocesses/stochprocess.jl") include("largedeviations/action.jl") include("largedeviations/min_action_method.jl") include("largedeviations/geometric_min_action_method.jl") -include("RateSystem.jl") +#include("RateSystem.jl") +include("RateSystemDraft.jl") include("../systems/CTLibrary.jl") using .CTLibrary @@ -80,7 +81,8 @@ 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 covariance_matrix, diffusion_matrix -export RateSystem +#export RateSystem +export RateProtocol,RateSystem # export edgetracking, bisect_to_edge, AttractorsViaProximity # export fixedpoints # ^ extention tests needed diff --git a/src/RateSystemDraft2.jl b/src/RateSystemDraft2.jl deleted file mode 100644 index 51694e5d..00000000 --- a/src/RateSystemDraft2.jl +++ /dev/null @@ -1,115 +0,0 @@ -using GLMakie -using CriticalTransitions -using DifferentialEquations - -# we consider the ODE dxₜ/dt = f(xₜ,λ(rt)) -# here λ = λ(t) ∈ Rᵐ is a function containing all the system parameters - -# We ask the user to define: -# 1) a system that should be investigated and -# 2) a protocol for λ(rt) by calling the function rate_protocol - - -############################################################################################################## -# we write behind the scenes -############################################################################################################## - -mutable struct RateProtocol - λ::Function - p_lambda::Vector - r::Float64 - t_start::Float64 - t_end::Float64 -end - -RateProtocol(λ::Function,p_lambda::Vector,r::Float64)=RateProtocol(λ,p_lambda,r,-Inf,Inf) -RateProtocol(λ::Function,r::Float64)=RateProtocol(λ,[],r,-Inf,Inf) -#RateProtocol(λ::Function,p_lambda::Vector,r::Float64,t_start::Float64)=RateProtocol(λ,p_lambda,r,t_start,Inf) -#RateProtocol(λ::Function,r::Float64,t_start::Float64)=RateProtocol(λ,[],r,t_start,Inf) - - -function modified_drift(u,p,t,ds::ContinuousTimeDynamicalSystem,λ::Function,t_start::Float64,t_end::Float64,r::Float64; - kwargs...) - - if t_start > t_end - error("Please ensure that t_start ≤ t_end.") - end - - p̃ = r*t ≤ t_start ? λ(p,t_start;kwargs...) : t_start < r*t < t_end ? λ(p,r*t;kwargs...) : λ(p,t_end;kwargs...); # the value(s) of λ(rt) - return ds.integ.f(u,p̃,t) -end; - -function RateSystem(auto_sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t0::Float64; - kwargs...) - # we wish to return a continuous time dynamical system with modified drift field - - f(u,p,t) = modified_drift(u,p,t,auto_sys,rp.λ,rp.t_start,rp.t_end,rp.r;kwargs...); - prob = remake(auto_sys.integ.sol.prob;f,p=rp.p_lambda,tspan=(t0,Inf)); - nonauto_sys = CoupledODEs(prob,auto_sys.diffeq); - return nonauto_sys -end; - -############################################################################################################## -# user writes (pseudo example - see the bottom of the page for an explicit example) -############################################################################################################## - - -############################################################################################################## -# test (prototypical model used for studying R-tipping, with critical rate r = 4/3) -############################################################################################################## - - -# autonomous system - -function f(u,p,t) - x = u[1] - λ = p[1] - dx = (x+λ)^2 - 1 - return SVector{1}(dx) -end; - -lambda = 0.; -p = [lambda]; - -z₀ = [-1.]; - -Δt = 1.e-3; -diffeq = (alg = AutoVern9(Rodas5(autodiff=false)),abstol=1.e-16,reltol=1.e-16); - -auto_sys = CoupledODEs(f,z₀,p;diffeq) - -# time-dependent parameter function - -function λ(p,t) - λ_max = p[1] - lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) - return SVector{1}(lambda) -end; - -# calling RateSystem - -r = 4/3-0.02; -λ_max = 3.; p_lambda = [λ_max]; - -#t_start = -Inf; t_end = Inf; -t0 = -10.; - -## test v2 - -rp = RateProtocol(λ,p_lambda,r); - -nonauto_sys = RateSystem(auto_sys,rp,t0); - -#nonauto_sys = RateSystem(auto_sys,λ,r,t0;p_lambda,t_start,t_end); - -T = 20.; - -auto_traj = trajectory(auto_sys,T,z₀); -nonauto_traj = trajectory(nonauto_sys,T,z₀); - -fig = Figure(); axs = Axis(fig[1,1]); -lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}"); -lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}"); -axislegend(axs,position=:rc,labelsize=40); -ylims!(-4.5,1.5); -fig \ No newline at end of file From 320fa82cc1d9d8f766ca2b4d67e6954b3b58b2f0 Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Wed, 19 Mar 2025 21:18:14 +0100 Subject: [PATCH 11/76] Corrected src/CriticalTransitions.jl file. --- src/CriticalTransitions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 53d2c35e..1fa7452b 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -63,7 +63,7 @@ include("largedeviations/action.jl") include("largedeviations/min_action_method.jl") include("largedeviations/geometric_min_action_method.jl") #include("RateSystem.jl") -include("RateSystemDraft.jl") +include("../dev/RateSystemDraft.jl") include("../systems/CTLibrary.jl") using .CTLibrary From 189f62b3f2bfd2065e3313e1ff4714224f025680 Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Wed, 19 Mar 2025 23:14:55 +0100 Subject: [PATCH 12/76] Moved contents of dev/ file to test/ file and reset src/CriticalTransitions.jl file to be in line with that of the main branch. --- src/CriticalTransitions.jl | 112 ++++++++++-------- test/{ => ratesystem}/RateSystem.jl | 0 {dev => test/ratesystem}/RateSystemDraft.jl | 0 .../ratesystem}/RateSystemTestsDraft.jl | 0 4 files changed, 62 insertions(+), 50 deletions(-) rename test/{ => ratesystem}/RateSystem.jl (100%) rename {dev => test/ratesystem}/RateSystemDraft.jl (100%) rename {dev => test/ratesystem}/RateSystemTestsDraft.jl (100%) diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 1fa7452b..6f8c7365 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -2,23 +2,15 @@ module CriticalTransitions # Base using Statistics: Statistics, mean -using LinearAlgebra: LinearAlgebra, I, norm, dot, tr +using LinearAlgebra: LinearAlgebra, I, norm, dot, tr, det using StaticArrays: StaticArrays, SVector +using SparseArrays: spdiagm +using DataStructures: CircularBuffer +using Random: Random # Core -using DiffEqNoiseProcess: DiffEqNoiseProcess -using OrdinaryDiffEq: OrdinaryDiffEq, Tsit5 -using StochasticDiffEq: - StochasticDiffEq, - DiscreteCallback, - ODEProblem, - SDEFunction, - SOSRA, - remake, - solve, - step!, - terminate!, - u_modified! +using SciMLBase: EnsembleThreads, DiscreteCallback, remake, terminate! +using StochasticDiffEq: StochasticDiffEq using DynamicalSystemsBase: DynamicalSystemsBase, CoupledSDEs, @@ -26,64 +18,84 @@ using DynamicalSystemsBase: dynamic_rule, current_state, set_state!, - trajectory + trajectory, + jacobian, + StateSpaceSet -using ForwardDiff: ForwardDiff -using IntervalArithmetic: IntervalArithmetic, interval -using Dierckx: Dierckx, ParametricSpline -using Optim: Optim, LBFGS -using Symbolics: Symbolics +using Interpolations: linear_interpolation +using Optimization +using OptimizationOptimisers: Optimisers +using LinearSolve: LinearProblem, KLUFactorization, solve # io and documentation using Format: Format -using Dates: Dates using Printf: Printf -using Markdown: Markdown -using DocStringExtensions: TYPEDSIGNATURES -using HDF5: HDF5, h5open, push! -using JLD2: JLD2, jldopen -using ProgressBars: ProgressBars, tqdm -using ProgressMeter: ProgressMeter +using DocStringExtensions: TYPEDSIGNATURES, TYPEDEF, TYPEDFIELDS, METHODLIST +using ProgressMeter: Progress, next! # reexport using Reexport: @reexport @reexport using StaticArrays -@reexport using StochasticDiffEq -@reexport using DiffEqNoiseProcess +@reexport using DynamicalSystemsBase -include("extention_functions.jl") +include("extension_functions.jl") include("utils.jl") -include("system_utils.jl") -include("io.jl") +include("sde_utils.jl") + +include("trajectories/TransitionEnsemble.jl") include("trajectories/simulation.jl") include("trajectories/transition.jl") -include("trajectories/equib.jl") -include("noiseprocesses/stochprocess.jl") + +include("transition_path_theory/TransitionPathMesh.jl") +include("transition_path_theory/langevin.jl") +include("transition_path_theory/committor.jl") +include("transition_path_theory/invariant_pdf.jl") +include("transition_path_theory/reactive_current.jl") +include("transition_path_theory/probability.jl") + +include("largedeviations/utils.jl") include("largedeviations/action.jl") +include("largedeviations/MinimumActionPath.jl") include("largedeviations/min_action_method.jl") include("largedeviations/geometric_min_action_method.jl") -#include("RateSystem.jl") -include("../dev/RateSystemDraft.jl") + +include("largedeviations/sgMAM.jl") +include("largedeviations/string_method.jl") include("../systems/CTLibrary.jl") using .CTLibrary +include("../test/ratesystem/RateSystemDraft.jl") + # Core types export CoupledSDEs, CoupledODEs, noise_process, covariance_matrix, diffusion_matrix -export dynamic_rule, current_state, set_state!, trajectory +export drift, div_drift, solver +export StateSpaceSet -# Methods -export drift, div_drift -export equilib, deterministic_orbit -export transition, transitions -export basins, basinboundary, basboundary +export sgmam, SgmamSystem 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 covariance_matrix, diffusion_matrix -#export RateSystem -export RateProtocol,RateSystem -# export edgetracking, bisect_to_edge, AttractorsViaProximity -# export fixedpoints -# ^ extention tests needed +export min_action_method, geometric_min_action_method, string_method +export MinimumActionPath + +export deterministic_orbit +export transition, transitions + +export distmesh2D, dellipse, ddiff +export TransitionPathMesh, Committor +export get_ellipse, reparametrization +export find_boundary, huniform, dunion + +export committor, + invariant_pdf, reactive_current, probability_reactive, probability_last_A, Langevin + +export RateProtocol, RateSystem + +# Error hint for extensions stubs +function __init__() + Base.Experimental.register_error_hint( + _basin_error_hinter(intervals_to_box), MethodError + ) + return nothing +end + end # module CriticalTransitions diff --git a/test/RateSystem.jl b/test/ratesystem/RateSystem.jl similarity index 100% rename from test/RateSystem.jl rename to test/ratesystem/RateSystem.jl diff --git a/dev/RateSystemDraft.jl b/test/ratesystem/RateSystemDraft.jl similarity index 100% rename from dev/RateSystemDraft.jl rename to test/ratesystem/RateSystemDraft.jl diff --git a/dev/RateSystemTestsDraft.jl b/test/ratesystem/RateSystemTestsDraft.jl similarity index 100% rename from dev/RateSystemTestsDraft.jl rename to test/ratesystem/RateSystemTestsDraft.jl From bba545bc9015774bdc636167c496fcfca115617a Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Wed, 19 Mar 2025 23:40:15 +0100 Subject: [PATCH 13/76] Collected all old and new RateSystem scripts into the test/ratesystem folder. --- src/RateSystem.jl => test/ratesystem/RateSystemDraft1.jl | 3 --- test/ratesystem/{RateSystemDraft.jl => RateSystemDraft2.jl} | 0 test/ratesystem/{RateSystem.jl => RateSystemTestDraft1.jl} | 2 +- .../{RateSystemTestsDraft.jl => RateSystemTestDraft2.jl} | 0 4 files changed, 1 insertion(+), 4 deletions(-) rename src/RateSystem.jl => test/ratesystem/RateSystemDraft1.jl (99%) rename test/ratesystem/{RateSystemDraft.jl => RateSystemDraft2.jl} (100%) rename test/ratesystem/{RateSystem.jl => RateSystemTestDraft1.jl} (94%) rename test/ratesystem/{RateSystemTestsDraft.jl => RateSystemTestDraft2.jl} (100%) diff --git a/src/RateSystem.jl b/test/ratesystem/RateSystemDraft1.jl similarity index 99% rename from src/RateSystem.jl rename to test/ratesystem/RateSystemDraft1.jl index b2dd11ed..019952ad 100644 --- a/src/RateSystem.jl +++ b/test/ratesystem/RateSystemDraft1.jl @@ -26,6 +26,3 @@ end; # Canard Trajectories # \dot(x) = f(x(t),λ(t)) - - - diff --git a/test/ratesystem/RateSystemDraft.jl b/test/ratesystem/RateSystemDraft2.jl similarity index 100% rename from test/ratesystem/RateSystemDraft.jl rename to test/ratesystem/RateSystemDraft2.jl diff --git a/test/ratesystem/RateSystem.jl b/test/ratesystem/RateSystemTestDraft1.jl similarity index 94% rename from test/ratesystem/RateSystem.jl rename to test/ratesystem/RateSystemTestDraft1.jl index 768b7794..4fb306f4 100644 --- a/test/ratesystem/RateSystem.jl +++ b/test/ratesystem/RateSystemTestDraft1.jl @@ -19,4 +19,4 @@ sys = RateSystem(tni,tnf,f,λ,p_λ,initvals) traj=trajectory(sys,40.,t0=-20.) trajic=[traj[1][i][1] for i in 1:length(traj[1])] -plot(collect(traj[2]),trajic) \ No newline at end of file +plot(collect(traj[2]),trajic) diff --git a/test/ratesystem/RateSystemTestsDraft.jl b/test/ratesystem/RateSystemTestDraft2.jl similarity index 100% rename from test/ratesystem/RateSystemTestsDraft.jl rename to test/ratesystem/RateSystemTestDraft2.jl From dcaf8aa577cac9150dafbc95e3ac52345fd04170 Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Thu, 20 Mar 2025 15:45:39 +0100 Subject: [PATCH 14/76] yep --- src/CriticalTransitions.jl | 2 +- systems/CTLibrary.jl | 3 +++ systems/truscottbrindley_mod_gen.jl | 22 ++++++++++++++++++++++ 3 files changed, 26 insertions(+), 1 deletion(-) create mode 100644 systems/truscottbrindley_mod_gen.jl diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 6f8c7365..52c4148b 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -65,7 +65,7 @@ include("largedeviations/string_method.jl") include("../systems/CTLibrary.jl") using .CTLibrary -include("../test/ratesystem/RateSystemDraft.jl") +include("../test/ratesystem/RateSystemDraft1.jl") # Core types export CoupledSDEs, CoupledODEs, noise_process, covariance_matrix, diffusion_matrix diff --git a/systems/CTLibrary.jl b/systems/CTLibrary.jl index a7e0ba8e..2c46fa4a 100644 --- a/systems/CTLibrary.jl +++ b/systems/CTLibrary.jl @@ -3,9 +3,11 @@ module CTLibrary using CriticalTransitions: CriticalTransitions using IntervalArithmetic: IntervalArithmetic, interval using StaticArrays: SA, SVector +using PolynomialRoots: roots include("fitzhughnagumo.jl") include("truscottbrindley_mod.jl") +include("truscottbrindley_mod_gen.jl") include("truscottbrindley_orig.jl") include("truscottbrindley_orig1.jl") include("rooth.jl") @@ -16,5 +18,6 @@ export fitzhugh_nagumo!, fitzhugh_nagumo, stommel, rivals!, rivals, cessi, rooth modifiedtruscottbrindleywithdimensions!, modifiedtruscottbrindleywithdimensions originaltruscottbrindley!, originaltruscottbrindley rampedoriginaltruscottbrindley!, rampedoriginaltruscottbrindley +truscottbrindley_mod_gen_det end diff --git a/systems/truscottbrindley_mod_gen.jl b/systems/truscottbrindley_mod_gen.jl new file mode 100644 index 00000000..1d6e328a --- /dev/null +++ b/systems/truscottbrindley_mod_gen.jl @@ -0,0 +1,22 @@ +function modTBgen_det(z,p,t) + x,y = z; # the system variables + α,β,ε = p; # the system parameters + # the drift field + dx = α*x*(1-x)-x^2*y/(β^2+x^2); + dy = ε*(x^2*y/(β^2+x^2)-y^2); + return SVector{2}([dx,dy]) +end; + +function truscottbrindley_mod_gen_det(p::Vector{Float64},diffeq; + t0::Float64 = 0., + kwargs...) + # finding the real equilibrium with greatest x-value + α,β = p[1:2]; + xeq = roots([α*β^4,-α*β^4,2*α*β^2,-1-2*α*β^2,α,-α]); + xeqreal = real.(xeq[abs.(imag.(xeq).<1.e-12)]); + xᵣ = findmax(xeqreal)[1]; + zᵣ = [xᵣ,α*(1-xᵣ)*(β^2+xᵣ^2)/xᵣ]; + # returning the CoupledODEs + return CoupledODEs(modTBgen_det,zᵣ,p;diffeq,t0,kwargs...) +end + From 650b6fa407c466d17258c8f7b8aa917183b3ab98 Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Thu, 20 Mar 2025 15:52:19 +0100 Subject: [PATCH 15/76] Added PolynomialRoots to dependencies. --- src/CriticalTransitions.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 52c4148b..537f7f15 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -26,6 +26,7 @@ using Interpolations: linear_interpolation using Optimization using OptimizationOptimisers: Optimisers using LinearSolve: LinearProblem, KLUFactorization, solve +using PolynomialRoots: roots # io and documentation using Format: Format From f6d97e325b7ec783b49013a06384c0021061834d Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Thu, 20 Mar 2025 16:03:17 +0100 Subject: [PATCH 16/76] See previous. --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 72ffbe03..4a2c4299 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +PolynomialRoots = "3a141323-8675-5d76-9d11-e1df1406c778" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" From 3febe829e0b86367467ae04aceda6e06547494aa Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Thu, 20 Mar 2025 16:05:24 +0100 Subject: [PATCH 17/76] See previous again. --- src/CriticalTransitions.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 537f7f15..52c4148b 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -26,7 +26,6 @@ using Interpolations: linear_interpolation using Optimization using OptimizationOptimisers: Optimisers using LinearSolve: LinearProblem, KLUFactorization, solve -using PolynomialRoots: roots # io and documentation using Format: Format From 8e0a7b3d4b82e0f19c740fb4969402975faeb4f4 Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Thu, 20 Mar 2025 16:15:18 +0100 Subject: [PATCH 18/76] Exported truscottbrindley-mod-gen-det function. --- systems/CTLibrary.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/systems/CTLibrary.jl b/systems/CTLibrary.jl index 2c46fa4a..d3ddb7f1 100644 --- a/systems/CTLibrary.jl +++ b/systems/CTLibrary.jl @@ -15,9 +15,9 @@ include("stommel.jl") include("rivals.jl") export fitzhugh_nagumo!, fitzhugh_nagumo, stommel, rivals!, rivals, cessi, rooth_smooth -modifiedtruscottbrindleywithdimensions!, modifiedtruscottbrindleywithdimensions -originaltruscottbrindley!, originaltruscottbrindley -rampedoriginaltruscottbrindley!, rampedoriginaltruscottbrindley -truscottbrindley_mod_gen_det +export truscottbrindley_mod_gen_det +#modifiedtruscottbrindleywithdimensions!, modifiedtruscottbrindleywithdimensions +#originaltruscottbrindley!, originaltruscottbrindley +#rampedoriginaltruscottbrindley!, rampedoriginaltruscottbrindley end From e2f006ef0868b63a7d1cdfe3af90c0bd886050cd Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Thu, 20 Mar 2025 16:44:51 +0100 Subject: [PATCH 19/76] Removed CoupledODEs functionality from systems/ folder. --- Project.toml | 1 - src/CriticalTransitions.jl | 2 +- systems/CTLibrary.jl | 1 - systems/truscottbrindley_mod_gen.jl | 18 ++---------------- 4 files changed, 3 insertions(+), 19 deletions(-) diff --git a/Project.toml b/Project.toml index 4a2c4299..72ffbe03 100644 --- a/Project.toml +++ b/Project.toml @@ -16,7 +16,6 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" -PolynomialRoots = "3a141323-8675-5d76-9d11-e1df1406c778" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 52c4148b..d55f483c 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -65,7 +65,7 @@ include("largedeviations/string_method.jl") include("../systems/CTLibrary.jl") using .CTLibrary -include("../test/ratesystem/RateSystemDraft1.jl") +include("../test/ratesystem/RateSystemDraft2.jl") # Core types export CoupledSDEs, CoupledODEs, noise_process, covariance_matrix, diffusion_matrix diff --git a/systems/CTLibrary.jl b/systems/CTLibrary.jl index d3ddb7f1..10bdd5af 100644 --- a/systems/CTLibrary.jl +++ b/systems/CTLibrary.jl @@ -3,7 +3,6 @@ module CTLibrary using CriticalTransitions: CriticalTransitions using IntervalArithmetic: IntervalArithmetic, interval using StaticArrays: SA, SVector -using PolynomialRoots: roots include("fitzhughnagumo.jl") include("truscottbrindley_mod.jl") diff --git a/systems/truscottbrindley_mod_gen.jl b/systems/truscottbrindley_mod_gen.jl index 1d6e328a..257fdd53 100644 --- a/systems/truscottbrindley_mod_gen.jl +++ b/systems/truscottbrindley_mod_gen.jl @@ -1,22 +1,8 @@ -function modTBgen_det(z,p,t) +function truscottbrindley_mod_gen_det(z,p,t) x,y = z; # the system variables α,β,ε = p; # the system parameters # the drift field dx = α*x*(1-x)-x^2*y/(β^2+x^2); dy = ε*(x^2*y/(β^2+x^2)-y^2); return SVector{2}([dx,dy]) -end; - -function truscottbrindley_mod_gen_det(p::Vector{Float64},diffeq; - t0::Float64 = 0., - kwargs...) - # finding the real equilibrium with greatest x-value - α,β = p[1:2]; - xeq = roots([α*β^4,-α*β^4,2*α*β^2,-1-2*α*β^2,α,-α]); - xeqreal = real.(xeq[abs.(imag.(xeq).<1.e-12)]); - xᵣ = findmax(xeqreal)[1]; - zᵣ = [xᵣ,α*(1-xᵣ)*(β^2+xᵣ^2)/xᵣ]; - # returning the CoupledODEs - return CoupledODEs(modTBgen_det,zᵣ,p;diffeq,t0,kwargs...) -end - +end; \ No newline at end of file From 74ba1ae76387292fec6c5f2865466cbbf801fce3 Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Thu, 20 Mar 2025 17:24:53 +0100 Subject: [PATCH 20/76] Removed the system I added. --- systems/CTLibrary.jl | 2 -- systems/truscottbrindley_mod_gen.jl | 8 -------- 2 files changed, 10 deletions(-) delete mode 100644 systems/truscottbrindley_mod_gen.jl diff --git a/systems/CTLibrary.jl b/systems/CTLibrary.jl index 10bdd5af..04ddc79a 100644 --- a/systems/CTLibrary.jl +++ b/systems/CTLibrary.jl @@ -6,7 +6,6 @@ using StaticArrays: SA, SVector include("fitzhughnagumo.jl") include("truscottbrindley_mod.jl") -include("truscottbrindley_mod_gen.jl") include("truscottbrindley_orig.jl") include("truscottbrindley_orig1.jl") include("rooth.jl") @@ -14,7 +13,6 @@ include("stommel.jl") include("rivals.jl") export fitzhugh_nagumo!, fitzhugh_nagumo, stommel, rivals!, rivals, cessi, rooth_smooth -export truscottbrindley_mod_gen_det #modifiedtruscottbrindleywithdimensions!, modifiedtruscottbrindleywithdimensions #originaltruscottbrindley!, originaltruscottbrindley #rampedoriginaltruscottbrindley!, rampedoriginaltruscottbrindley diff --git a/systems/truscottbrindley_mod_gen.jl b/systems/truscottbrindley_mod_gen.jl deleted file mode 100644 index 257fdd53..00000000 --- a/systems/truscottbrindley_mod_gen.jl +++ /dev/null @@ -1,8 +0,0 @@ -function truscottbrindley_mod_gen_det(z,p,t) - x,y = z; # the system variables - α,β,ε = p; # the system parameters - # the drift field - dx = α*x*(1-x)-x^2*y/(β^2+x^2); - dy = ε*(x^2*y/(β^2+x^2)-y^2); - return SVector{2}([dx,dy]) -end; \ No newline at end of file From bdd468e51ce140c46a0b9e92a6cd043e0443db26 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Sat, 28 Jun 2025 12:44:09 +0200 Subject: [PATCH 21/76] Starting a RateSystem Example --- Project.toml | 1 + docs/src/examples/RateSystem.md | 247 ++++++++++++++++++++++++++++++++ 2 files changed, 248 insertions(+) create mode 100644 docs/src/examples/RateSystem.md diff --git a/Project.toml b/Project.toml index a05bfb6b..fde25719 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ repo = "https://github.com/juliadynamics/CriticalTransitions.jl.git" version = "0.6.0" [deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4" diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md new file mode 100644 index 00000000..a95534ff --- /dev/null +++ b/docs/src/examples/RateSystem.md @@ -0,0 +1,247 @@ + +# RateSystem + +```@example GMAM +using CriticalTransitions + +using CairoMakie +using CairoMakie.Makie.MathTeXEngine: get_font +font = (; + regular=get_font(:regular), bold=get_font(:bold), + italic=get_font(:italic), bold_italic=get_font(:bolditalic) +); +``` + +Let us explore the features of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl) with Maier-Stein model. + +## Maier-stein model + +The [Maier-Stein model](https://doi.org/10.1007/BF02183736) (J. Stat. Phys 83, 3–4 (1996)) is commonly used in the field of nonlinear dynamics for benchmarking Large Deviation Theory (LDT) techniques, e.g., stoachastic transitions between different stable states. It is a simple model that describes the dynamics of a system with two degrees of freedom ``u`` and ``v``, and is given by the following set of ordinary differential equations: +```math +\begin{aligned} + \dot{u} &= u-u^3 - \beta*u*v^2\\ + \dot{v} &= -\alpha (1+u^2)*v +\end{aligned} +``` +The parameter ``\alpha>0`` controls the strength of the drift field and ``\beta>0`` represents the softening of that drift field. + +```@example GMAM +function meier_stein!(du, u, p, t) # in-place + x, y = u + du[1] = x - x^3 - 10 * x * y^2 + du[2] = -(1 + x^2) * y +end +function meier_stein(u, p, t) # out-of-place + x, y = u + dx = x - x^3 - 10 * x * y^2 + dy = -(1 + x^2) * y + SA[dx, dy] +end +σ = 0.25 +sys = CoupledSDEs(meier_stein, zeros(2), (); noise_strength=σ) +``` + +A good reference to read about the large deviations methods is [this](https://homepages.warwick.ac.uk/staff/T.Grafke/simplified-geometric-minimum-action-method-for-the-computation-of-instantons.html) or [this](https://homepages.warwick.ac.uk/staff/T.Grafke/simplified-geometric-minimum-action-method-for-the-computation-of-instantons.html) blog post by Tobias Grafke. + +## Attractors + +We start by investigating the deterministic dynamics of the Maier-Stein model. + +The function `fixed points` return the attractors, their eigenvalues and stability within the state space volume defined by `bmin` and `bmax`. + +```@example GMAM +using ChaosTools + +u_min = -1.1; +u_max = 1.1; +v_min = -0.4; +v_max = 0.4; +bmin = [u_min, v_min]; +bmax = [u_max, v_max]; +fp, eig, stab = fixedpoints(sys, bmin, bmax) +stable_fp = fp[stab] +``` + +```@example GMAM +using LinearAlgebra: norm +res = 100 +u_range = range(u_min, u_max, length=res) +v_range = range(v_min, v_max, length=res) + +du(u, v) = u - u^3 - 10 * u * v^2 +dv(u, v) = -(1 + u^2) * v +odeSol(u, v) = Point2f(du(u, v), dv(u, v)) + +z = [norm([du(x, y), dv(x, y)]) for x in u_range, y in v_range] +zmin, zmax = minimum(z), maximum(z) + +fig = Figure(size=(600, 400), fontsize=13) +ax = Axis(fig[1, 1], xlabel="u", ylabel="v", aspect=1.4, + xgridcolor=:transparent, ygridcolor=:transparent, + ylabelrotation=0) + +hm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax)) +Colorbar(fig[1, 2], hm; label="", width=15, ticksize=15, tickalign=1) +streamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max); + gridsize=(20, 20), arrow_size=10, stepsize=0.01, + colormap=[:black, :black] +) +colgap!(fig.layout, 7) +limits!(u_min, u_max, v_min, v_max) +fig + +[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue, + markersize=10) for i in eachindex(fp)] +fig +``` + +We can simulate a stochastic trajectory using the function `trajectory`. + +```@example GMAM +tr, ts = trajectory(sys, 1000) + +fig = Figure(size=(1000, 400), fontsize=13) +ax1 = Axis(fig[1, 1], xlabel="t", ylabel="u", aspect=1.2, + xgridcolor=:transparent, ygridcolor=:transparent, + ylabelrotation=0) +ax2 = Axis(fig[1, 2], xlabel="u", ylabel="v", aspect=1.2, + xgridcolor=:transparent, ygridcolor=:transparent, + ylabelrotation=0) + +lines!(ax1, ts, first.(tr); linewidth=2, color=:black) + +hm = heatmap!(ax2, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax)) +Colorbar(fig[1, 3], hm; label="", width=15, ticksize=15, tickalign=1) +streamplot!(ax2, odeSol, (u_min, u_max), (v_min, v_max); + gridsize=(20, 20), arrow_size=10, stepsize=0.01, + colormap=[:white, :white] +) +colgap!(fig.layout, 7) +limits!(u_min, u_max, v_min, v_max) +fig + +[scatter!(ax2, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue, + markersize=10) for i in eachindex(fp)] + +lines!(ax2, reduce(hcat, tr), linewidth=1, color=(:black, 0.2)) +fig +``` + +## Transitions + +We can quickly find a path which computes a transition from one attractor to another using the function `transition. + +```@example GMAM +paths_ends = (fp[stab][1], fp[stab][2]) +path, time, succes = transition(sys, paths_ends...); +``` + +```@example GMAM +fig = Figure(size=(600, 400), fontsize=13) +ax = Axis(fig[1, 1], xlabel="u", ylabel="v", aspect=1.4, + xgridcolor=:transparent, ygridcolor=:transparent, + ylabelrotation=0) + +hm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax)) +Colorbar(fig[1, 2], hm; label="", width=15, ticksize=15, tickalign=1) +streamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max); + gridsize=(20, 20), arrow_size=10, stepsize=0.01, + colormap=[:white, :white] +) +colgap!(fig.layout, 7) +limits!(u_min, u_max, v_min, v_max) +fig + +[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue, + markersize=10) for i in eachindex(fp)] +fig + +lines!(ax, path, color=:black) +fig +``` + +If we want to compute many: `transitions` is the function to use. + +```@example GMAM +tt = transitions(sys, paths_ends..., 3; tmax=1e3); +``` + +```@example GMAM +fig = Figure(size=(600, 400), fontsize=13) +ax = Axis(fig[1, 1], xlabel="u", ylabel="v", aspect=1.4, + xgridcolor=:transparent, ygridcolor=:transparent, + ylabelrotation=0) + +hm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax)) +Colorbar(fig[1, 2], hm; label="", width=15, ticksize=15, tickalign=1) +streamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max); + gridsize=(20, 20), arrow_size=10, stepsize=0.01, + colormap=[:black, :black] +) +colgap!(fig.layout, 7) +limits!(u_min, u_max, v_min, v_max) +fig + +[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue, + markersize=10) for i in eachindex(fp)] + +for i in 1:length(tt.paths) + lines!(ax, tt.paths[i]) +end +fig +``` + +## Large deviation theory + +In the context of nonlinear dynamics, Large Deviation Theory provides tools to quantify the probability of rare events that deviate significantly from the system's typical behavior. These rare events might be extreme values of a system's output, sudden transitions between different states, or other phenomena that occur with very low probability but can have significant implications for the system's overall behavior. + +Large deviation theory applies principles from probability theory and statistical mechanics to develop a rigorous mathematical description of these rare events. It uses the concept of a rate function, which measures the exponential decay rate of the probability of large deviations from the mean or typical behavior. This rate function plays a crucial role in quantifying the likelihood of rare events and understanding their impact on the system. + +For example, in a system exhibiting chaotic behavior, LDT can help quantify the probability of sudden large shifts in the system's trajectory. Similarly, in a system with multiple stable states, it can provide insight into the likelihood and pathways of transitions between these states under fluctuations. In the context of the Minimum Action Method (MAM) and the Geometric Minimum Action Method (gMAM), Large Deviation Theory is used to handle the large deviations action functional on the space of curves. This is a key part of how these methods analyze dynamical systems. + +The Maier-Stein model is a typical benchmark to test such LDT techniques. Let us try to reproduce the following figure from [Tobias Grafke's blog post](https://homepages.warwick.ac.uk/staff/T.Grafke/simplified-geometric-minimum-action-method-for-the-computation-of-instantons.html): + +![maier_stein](./maierstein-dynamics.png) + +Let us first make an initial path: + +```@example GMAM +xx = range(-1.0, 1.0, length=100) +yy = 0.3 .* (-xx .^ 2 .+ 1) +init = Matrix([xx yy]') +``` + +`geometric_min_action_method` computes the minimizer of the Freidlin-Wentzell action using the geometric minimum action method (gMAM), to find the minimum action path (instanton) between an initial state x_i and final state x_f. The Minimum Action Method (MAM) is a more traditional approach, while the Geometric Minimum Action Method (gMAM) is a blend of the original MAM and the [string method](https://doi.org/10.1103/PhysRevB.66.052301). + +```@example GMAM +gm = geometric_min_action_method(sys, init, maxiter=500; show_progress=false) +MLP = gm.path +``` + +```@example GMAM +fig = Figure(size=(600, 400), fontsize=13) +ax = Axis(fig[1, 1], xlabel="u", ylabel="v", aspect=1.4, + xgridcolor=:transparent, ygridcolor=:transparent, + ylabelrotation=0) + +hm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax)) +Colorbar(fig[1, 2], hm; label="", width=15, ticksize=15, tickalign=1) +streamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max); + gridsize=(20, 20), arrow_size=10, stepsize=0.01, + colormap=[:black, :black] +) +colgap!(fig.layout, 7) +limits!(u_min, u_max, v_min, v_max) +fig + +[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue, + markersize=10) for i in eachindex(fp)] + +lines!(ax, init, linewidth=3, color=:black, linestyle=:dash) +lines!(ax, MLP, linewidth=3, color=:orange) +fig +``` + +--- +Author: Orjan Ameye (orjan.ameye@hotmail.com) +Date: 13 Feb 2024 From 020b80f076aaec895e0c2f5e5b224d20df99b375 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Sun, 29 Jun 2025 08:42:43 +0200 Subject: [PATCH 22/76] work on example of RateSystem --- docs/src/examples/RateSystem.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index a95534ff..2d6227ee 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -1,7 +1,7 @@ # RateSystem -```@example GMAM +```@example RateSystem using CriticalTransitions using CairoMakie @@ -14,9 +14,9 @@ font = (; Let us explore the features of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl) with Maier-Stein model. -## Maier-stein model +## Prototypical model for R-tipping, with critical rate r = 4/3 -The [Maier-Stein model](https://doi.org/10.1007/BF02183736) (J. Stat. Phys 83, 3–4 (1996)) is commonly used in the field of nonlinear dynamics for benchmarking Large Deviation Theory (LDT) techniques, e.g., stoachastic transitions between different stable states. It is a simple model that describes the dynamics of a system with two degrees of freedom ``u`` and ``v``, and is given by the following set of ordinary differential equations: +The following simple one-dimensional model with one attractor is given by the following ordinary differential equations: ```math \begin{aligned} \dot{u} &= u-u^3 - \beta*u*v^2\\ From 894945b8c8d0b789aa3d56a949856109e3d4ccb5 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Sun, 29 Jun 2025 12:11:48 +0200 Subject: [PATCH 23/76] work on example of RateSystem --- docs/src/examples/RateSystem.md | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index 2d6227ee..b7c11c73 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -19,11 +19,10 @@ Let us explore the features of [CriticalTransitions.jl](https://github.com/Julia The following simple one-dimensional model with one attractor is given by the following ordinary differential equations: ```math \begin{aligned} - \dot{u} &= u-u^3 - \beta*u*v^2\\ - \dot{v} &= -\alpha (1+u^2)*v + \dot{x} &= (x+\lambda)^2 - 1 \end{aligned} ``` -The parameter ``\alpha>0`` controls the strength of the drift field and ``\beta>0`` represents the softening of that drift field. +The parameter ``\lambda`` shifts the location of the extrema of the drift field. ```@example GMAM function meier_stein!(du, u, p, t) # in-place From af0efc47aca39aad336c21797718c89b4ce70e08 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Sun, 29 Jun 2025 12:20:17 +0200 Subject: [PATCH 24/76] work on example of RateSystem --- docs/src/examples/RateSystem.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index b7c11c73..6cf7a9af 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -132,7 +132,7 @@ We can quickly find a path which computes a transition from one attractor to ano ```@example GMAM paths_ends = (fp[stab][1], fp[stab][2]) -path, time, succes = transition(sys, paths_ends...); +path, time, success = transition(sys, paths_ends...); ``` ```@example GMAM From 9e1f68dcbdff9710f6e28d11a247807cd19b5373 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Sun, 29 Jun 2025 12:31:02 +0200 Subject: [PATCH 25/76] work on example of RateSystem --- docs/src/examples/RateSystem.md | 44 ++++++++++++++++++++++----------- 1 file changed, 30 insertions(+), 14 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index 6cf7a9af..f3788e7e 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -24,23 +24,39 @@ The following simple one-dimensional model with one attractor is given by the fo ``` The parameter ``\lambda`` shifts the location of the extrema of the drift field. -```@example GMAM -function meier_stein!(du, u, p, t) # in-place - x, y = u - du[1] = x - x^3 - 10 * x * y^2 - du[2] = -(1 + x^2) * y -end -function meier_stein(u, p, t) # out-of-place - x, y = u - dx = x - x^3 - 10 * x * y^2 - dy = -(1 + x^2) * y - SA[dx, dy] +```@example RateSystem +function f(u,p,t) # out-of-place + x = u[1] + λ = p[1] + dx = (x+λ)^2 - 1 + return SVector{1}(dx) end -σ = 0.25 -sys = CoupledSDEs(meier_stein, zeros(2), (); noise_strength=σ) +lambda = 0.0 +p = [lambda] +z₀ = [-1.] +auto_sys = CoupledODEs(f,z₀,p) ``` -A good reference to read about the large deviations methods is [this](https://homepages.warwick.ac.uk/staff/T.Grafke/simplified-geometric-minimum-action-method-for-the-computation-of-instantons.html) or [this](https://homepages.warwick.ac.uk/staff/T.Grafke/simplified-geometric-minimum-action-method-for-the-computation-of-instantons.html) blog post by Tobias Grafke. +## Non-autonomous case + +Now, we define the time-dependent parameter function ``\lambda(t)`` to make the system non-autonomous and investigate the system's behaviour under parameter rampings. + +```@example RateSystem +function λ(p,t) + λ_max = p[1] + lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) + return SVector{1}(lambda) +end +``` +We define the following parameters +``` +r = 4/3-0.02 +λ_max = 3. +p_lambda = [λ_max] +t_start = -Inf +t_end = Inf +t0 = -10. +``` ## Attractors From d3c9ecef9de7bcb564fad1a2b88097be540a199c Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Sun, 29 Jun 2025 12:36:55 +0200 Subject: [PATCH 26/76] added documentation in RateSystem file --- test/ratesystem/RateSystemDraft2.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/ratesystem/RateSystemDraft2.jl b/test/ratesystem/RateSystemDraft2.jl index 41c8ee44..04126a4d 100644 --- a/test/ratesystem/RateSystemDraft2.jl +++ b/test/ratesystem/RateSystemDraft2.jl @@ -2,9 +2,11 @@ # here λ = λ(t) ∈ Rᵐ is a function containing all the system parameters # We ask the user to define: -# 1) a system that should be investigated and +# 1) a ContinuousTimeDynamicalSystem that should be investigated and # 2) a protocol for the time-dependent forcing with the struct RateProtocol +# Then we give back the ContinuousTimeDynamicalSystem with the parameter +# changing according to the rate protocol mutable struct RateProtocol λ::Function p_lambda::Vector From e3d2156ed71fca8d4ab1ef62a70f92e7926bb209 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 30 Jun 2025 22:08:21 +0200 Subject: [PATCH 27/76] corrected mistake in RateSystem file --- test/ratesystem/RateSystemDraft2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/ratesystem/RateSystemDraft2.jl b/test/ratesystem/RateSystemDraft2.jl index 04126a4d..838b24f7 100644 --- a/test/ratesystem/RateSystemDraft2.jl +++ b/test/ratesystem/RateSystemDraft2.jl @@ -41,4 +41,4 @@ function RateSystem(auto_sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t prob = remake(auto_sys.integ.sol.prob;f,p=rp.p_lambda,tspan=(t0,Inf)); nonauto_sys = CoupledODEs(prob,auto_sys.diffeq); return nonauto_sys -end; +end From 90f42283a4f1b1acc582ef5faae8fd8212b6fe8c Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 30 Jun 2025 22:23:36 +0200 Subject: [PATCH 28/76] improved documentation of test2 --- test/ratesystem/RateSystemTestDraft2.jl | 72 ++++++++++++------------- 1 file changed, 33 insertions(+), 39 deletions(-) diff --git a/test/ratesystem/RateSystemTestDraft2.jl b/test/ratesystem/RateSystemTestDraft2.jl index 2decbdf4..42921043 100644 --- a/test/ratesystem/RateSystemTestDraft2.jl +++ b/test/ratesystem/RateSystemTestDraft2.jl @@ -3,56 +3,50 @@ ############################################################################################################## # autonomous system - function f(u,p,t) x = u[1] λ = p[1] dx = (x+λ)^2 - 1 return SVector{1}(dx) -end; +end -lambda = 0.; -p = [lambda]; - -z₀ = [-1.]; +lambda = 0. +p = [lambda] +x0 = [-1.] -Δt = 1.e-3; -diffeq = (alg = AutoVern9(Rodas5(autodiff=false)),abstol=1.e-16,reltol=1.e-16); +Δt = 1.e-3 +#diffeq = (alg = AutoVern9(Rodas5(autodiff=false)),abstol=1.e-16,reltol=1.e-16) +#auto_sys = CoupledODEs(f,x0,p;diffeq) +auto_sys = CoupledODEs(f,x0,p) -auto_sys = CoupledODEs(f,z₀,p;diffeq) -# time-dependent parameter function +# now we make the system non-autonomous +# We define a time-dependent parameter function function λ(p,t) λ_max = p[1] lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) return SVector{1}(lambda) -end; - -# calling RateSystem - -r = 4/3-0.02; -λ_max = 3.; p_lambda = [λ_max]; - -#t_start = -Inf; t_end = Inf; -t0 = -10.; - -## test v2 - -rp = RateProtocol(λ,p_lambda,r); - -nonauto_sys = RateSystem(auto_sys,rp,t0); - -#nonauto_sys = RateSystem(auto_sys,λ,r,t0;p_lambda,t_start,t_end); - -T = 20.; - -auto_traj = trajectory(auto_sys,T,z₀); -nonauto_traj = trajectory(nonauto_sys,T,z₀); - -# fig = Figure(); axs = Axis(fig[1,1]); -# lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}"); -# lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}"); -# axislegend(axs,position=:rc,labelsize=40); -# ylims!(-4.5,1.5); -# fig \ No newline at end of file +end + +# we set the parameters for RateProtocol +λ_max = 3. +p_lambda = [λ_max] +r = 4/3-0.02 # r just below critical rate +t_start = -Inf # start time of non-autonomous part +t_end = Inf # end time of non-autonomous part +# And the initial time of the system +t0 = -10. + +rp = RateProtocol(λ,p_lambda,r,t_start,t_end) +nonauto_sys = RateSystem(auto_sys,rp,t0) + +T = 20. # final simulation time +auto_traj = trajectory(auto_sys,T,x0) +nonauto_traj = trajectory(nonauto_sys,T,x0) + +fig = Figure(); axs = Axis(fig[1,1]) +lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}") +lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}") +axislegend(axs,position=:rc,labelsize=10) +fig \ No newline at end of file From 1703378a149f8e9a77a4a4740761877137a25d47 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 30 Jun 2025 22:34:52 +0200 Subject: [PATCH 29/76] work on the example for RateSystem --- docs/src/examples/RateSystem.md | 229 ++++---------------------------- 1 file changed, 29 insertions(+), 200 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index f3788e7e..2c07aa61 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -3,7 +3,8 @@ ```@example RateSystem using CriticalTransitions - +using DifferentialEquations +using ModelingToolkit using CairoMakie using CairoMakie.Makie.MathTeXEngine: get_font font = (; @@ -12,7 +13,7 @@ font = (; ); ``` -Let us explore the features of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl) with Maier-Stein model. +Let us explore an example of the RateSystem setting of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl). ## Prototypical model for R-tipping, with critical rate r = 4/3 @@ -33,8 +34,8 @@ function f(u,p,t) # out-of-place end lambda = 0.0 p = [lambda] -z₀ = [-1.] -auto_sys = CoupledODEs(f,z₀,p) +x0 = [-1.] +auto_sys = CoupledODEs(f,x0,p) ``` ## Non-autonomous case @@ -49,214 +50,42 @@ function λ(p,t) end ``` We define the following parameters -``` -r = 4/3-0.02 +```@example RateSystem λ_max = 3. -p_lambda = [λ_max] -t_start = -Inf -t_end = Inf +p_lambda = [λ_max] # parameter of the function lambda +r = 4/3-0.02 # r just below critical rate +t_start = -Inf # start time of non-autonomous part +t_end = Inf # end time of non-autonomous part +# And the initial time of the system t0 = -10. ``` -## Attractors - -We start by investigating the deterministic dynamics of the Maier-Stein model. - -The function `fixed points` return the attractors, their eigenvalues and stability within the state space volume defined by `bmin` and `bmax`. - -```@example GMAM -using ChaosTools - -u_min = -1.1; -u_max = 1.1; -v_min = -0.4; -v_max = 0.4; -bmin = [u_min, v_min]; -bmax = [u_max, v_max]; -fp, eig, stab = fixedpoints(sys, bmin, bmax) -stable_fp = fp[stab] -``` - -```@example GMAM -using LinearAlgebra: norm -res = 100 -u_range = range(u_min, u_max, length=res) -v_range = range(v_min, v_max, length=res) - -du(u, v) = u - u^3 - 10 * u * v^2 -dv(u, v) = -(1 + u^2) * v -odeSol(u, v) = Point2f(du(u, v), dv(u, v)) - -z = [norm([du(x, y), dv(x, y)]) for x in u_range, y in v_range] -zmin, zmax = minimum(z), maximum(z) - -fig = Figure(size=(600, 400), fontsize=13) -ax = Axis(fig[1, 1], xlabel="u", ylabel="v", aspect=1.4, - xgridcolor=:transparent, ygridcolor=:transparent, - ylabelrotation=0) - -hm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax)) -Colorbar(fig[1, 2], hm; label="", width=15, ticksize=15, tickalign=1) -streamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max); - gridsize=(20, 20), arrow_size=10, stepsize=0.01, - colormap=[:black, :black] -) -colgap!(fig.layout, 7) -limits!(u_min, u_max, v_min, v_max) -fig - -[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue, - markersize=10) for i in eachindex(fp)] -fig -``` - -We can simulate a stochastic trajectory using the function `trajectory`. - -```@example GMAM -tr, ts = trajectory(sys, 1000) - -fig = Figure(size=(1000, 400), fontsize=13) -ax1 = Axis(fig[1, 1], xlabel="t", ylabel="u", aspect=1.2, - xgridcolor=:transparent, ygridcolor=:transparent, - ylabelrotation=0) -ax2 = Axis(fig[1, 2], xlabel="u", ylabel="v", aspect=1.2, - xgridcolor=:transparent, ygridcolor=:transparent, - ylabelrotation=0) - -lines!(ax1, ts, first.(tr); linewidth=2, color=:black) - -hm = heatmap!(ax2, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax)) -Colorbar(fig[1, 3], hm; label="", width=15, ticksize=15, tickalign=1) -streamplot!(ax2, odeSol, (u_min, u_max), (v_min, v_max); - gridsize=(20, 20), arrow_size=10, stepsize=0.01, - colormap=[:white, :white] -) -colgap!(fig.layout, 7) -limits!(u_min, u_max, v_min, v_max) -fig - -[scatter!(ax2, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue, - markersize=10) for i in eachindex(fp)] - -lines!(ax2, reduce(hcat, tr), linewidth=1, color=(:black, 0.2)) -fig -``` - -## Transitions - -We can quickly find a path which computes a transition from one attractor to another using the function `transition. - -```@example GMAM -paths_ends = (fp[stab][1], fp[stab][2]) -path, time, success = transition(sys, paths_ends...); -``` - -```@example GMAM -fig = Figure(size=(600, 400), fontsize=13) -ax = Axis(fig[1, 1], xlabel="u", ylabel="v", aspect=1.4, - xgridcolor=:transparent, ygridcolor=:transparent, - ylabelrotation=0) - -hm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax)) -Colorbar(fig[1, 2], hm; label="", width=15, ticksize=15, tickalign=1) -streamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max); - gridsize=(20, 20), arrow_size=10, stepsize=0.01, - colormap=[:white, :white] -) -colgap!(fig.layout, 7) -limits!(u_min, u_max, v_min, v_max) -fig - -[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue, - markersize=10) for i in eachindex(fp)] -fig - -lines!(ax, path, color=:black) -fig -``` - -If we want to compute many: `transitions` is the function to use. - -```@example GMAM -tt = transitions(sys, paths_ends..., 3; tmax=1e3); -``` - -```@example GMAM -fig = Figure(size=(600, 400), fontsize=13) -ax = Axis(fig[1, 1], xlabel="u", ylabel="v", aspect=1.4, - xgridcolor=:transparent, ygridcolor=:transparent, - ylabelrotation=0) +We define the RateProtocol -hm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax)) -Colorbar(fig[1, 2], hm; label="", width=15, ticksize=15, tickalign=1) -streamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max); - gridsize=(20, 20), arrow_size=10, stepsize=0.01, - colormap=[:black, :black] -) -colgap!(fig.layout, 7) -limits!(u_min, u_max, v_min, v_max) -fig - -[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue, - markersize=10) for i in eachindex(fp)] - -for i in 1:length(tt.paths) - lines!(ax, tt.paths[i]) -end -fig +```@example RateSystem +rp = RateProtocol(λ,p_lambda,r,t_start,t_end) ``` -## Large deviation theory - -In the context of nonlinear dynamics, Large Deviation Theory provides tools to quantify the probability of rare events that deviate significantly from the system's typical behavior. These rare events might be extreme values of a system's output, sudden transitions between different states, or other phenomena that occur with very low probability but can have significant implications for the system's overall behavior. - -Large deviation theory applies principles from probability theory and statistical mechanics to develop a rigorous mathematical description of these rare events. It uses the concept of a rate function, which measures the exponential decay rate of the probability of large deviations from the mean or typical behavior. This rate function plays a crucial role in quantifying the likelihood of rare events and understanding their impact on the system. - -For example, in a system exhibiting chaotic behavior, LDT can help quantify the probability of sudden large shifts in the system's trajectory. Similarly, in a system with multiple stable states, it can provide insight into the likelihood and pathways of transitions between these states under fluctuations. In the context of the Minimum Action Method (MAM) and the Geometric Minimum Action Method (gMAM), Large Deviation Theory is used to handle the large deviations action functional on the space of curves. This is a key part of how these methods analyze dynamical systems. - -The Maier-Stein model is a typical benchmark to test such LDT techniques. Let us try to reproduce the following figure from [Tobias Grafke's blog post](https://homepages.warwick.ac.uk/staff/T.Grafke/simplified-geometric-minimum-action-method-for-the-computation-of-instantons.html): - -![maier_stein](./maierstein-dynamics.png) - -Let us first make an initial path: +And use it to create the system with past and future autonomous systems and non-autonomous ramping in between: -```@example GMAM -xx = range(-1.0, 1.0, length=100) -yy = 0.3 .* (-xx .^ 2 .+ 1) -init = Matrix([xx yy]') -``` - -`geometric_min_action_method` computes the minimizer of the Freidlin-Wentzell action using the geometric minimum action method (gMAM), to find the minimum action path (instanton) between an initial state x_i and final state x_f. The Minimum Action Method (MAM) is a more traditional approach, while the Geometric Minimum Action Method (gMAM) is a blend of the original MAM and the [string method](https://doi.org/10.1103/PhysRevB.66.052301). +```@example RateSystem +nonauto_sys = RateSystem(auto_sys,rp,t0) -```@example GMAM -gm = geometric_min_action_method(sys, init, maxiter=500; show_progress=false) -MLP = gm.path +T = 20. # final simulation time +auto_traj = trajectory(auto_sys,T,x0) +nonauto_traj = trajectory(nonauto_sys,T,x0) ``` -```@example GMAM -fig = Figure(size=(600, 400), fontsize=13) -ax = Axis(fig[1, 1], xlabel="u", ylabel="v", aspect=1.4, - xgridcolor=:transparent, ygridcolor=:transparent, - ylabelrotation=0) +We plot the two trajectories -hm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax)) -Colorbar(fig[1, 2], hm; label="", width=15, ticksize=15, tickalign=1) -streamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max); - gridsize=(20, 20), arrow_size=10, stepsize=0.01, - colormap=[:black, :black] -) -colgap!(fig.layout, 7) -limits!(u_min, u_max, v_min, v_max) -fig - -[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue, - markersize=10) for i in eachindex(fp)] - -lines!(ax, init, linewidth=3, color=:black, linestyle=:dash) -lines!(ax, MLP, linewidth=3, color=:orange) +```@example RateSystem +fig = Figure(); axs = Axis(fig[1,1]) +lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}") +lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}") +axislegend(axs,position=:rc,labelsize=10) fig ``` ---- -Author: Orjan Ameye (orjan.ameye@hotmail.com) -Date: 13 Feb 2024 +----- +Author: Raphael Roemer +Date: 30 Jun 202r From 2a5619e8bb6870c4b3aca6470f7207f88983bfb5 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 30 Jun 2025 22:48:48 +0200 Subject: [PATCH 30/76] added codes for the RateSystem example --- docs/src/examples/RateSystem.md | 4 +- examples/RateSystem.jl | 96 +++++++++++++++++++++++++++++++++ 2 files changed, 99 insertions(+), 1 deletion(-) create mode 100644 examples/RateSystem.jl diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index 2c07aa61..7f06564f 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -1,4 +1,6 @@ - +```@meta +EditURL = "../../../examples/RateSystem.jl" +``` # RateSystem ```@example RateSystem diff --git a/examples/RateSystem.jl b/examples/RateSystem.jl new file mode 100644 index 00000000..4359119f --- /dev/null +++ b/examples/RateSystem.jl @@ -0,0 +1,96 @@ + +# # Rate System Example + +using CriticalTransitions +using DifferentialEquations +using ModelingToolkit +using CairoMakie +using CairoMakie.Makie.MathTeXEngine: get_font +font = (; + regular=get_font(:regular), bold=get_font(:bold), + italic=get_font(:italic), bold_italic=get_font(:bolditalic) +); + +# Let us explore an example of the RateSystem setting of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl). + +# ## Prototypical model for R-tipping, with critical rate r = 4/3 + +# The following simple one-dimensional model with one attractor is given by the following ordinary differential equations: +# ```math +# \begin{aligned} +# \dot{x} &= (x+\lambda)^2 - 1 +# \end{aligned} +# ``` +# The parameter ``\lambda`` shifts the location of the extrema of the drift field. + +function f(u,p,t) # out-of-place + x = u[1] + λ = p[1] + dx = (x+λ)^2 - 1 + return SVector{1}(dx) +end +lambda = 0.0 +p = [lambda] +x0 = [-1.] +auto_sys = CoupledODEs(f,x0,p) + +## Non-autonomous case + +# Now, we define the time-dependent parameter function ``\lambda(t)`` to make the system non-autonomous and investigate the system's behaviour under parameter rampings. + +# ```@example RateSystem +# function λ(p,t) +# λ_max = p[1] +# lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) +# return SVector{1}(lambda) +# end +# ``` + +function λ(p,t) + λ_max = p[1] + lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) + return SVector{1}(lambda) +end + +# We define the following parameters +# ```@example RateSystem +# λ_max = 3. +# p_lambda = [λ_max] # parameter of the function lambda +# r = 4/3-0.02 # r just below critical rate +# t_start = -Inf # start time of non-autonomous part +# t_end = Inf # end time of non-autonomous part +# # And the initial time of the system +# t0 = -10. +# ``` + +λ_max = 3. +p_lambda = [λ_max] # parameter of the function lambda +r = 4/3-0.02 # r just below critical rate +t_start = -Inf # start time of non-autonomous part +t_end = Inf # end time of non-autonomous part +# And the initial time of the system +t0 = -10. + + +# We define the RateProtocol + +# ```@example RateSystem +# rp = RateProtocol(λ,p_lambda,r,t_start,t_end) +# ``` +rp = RateProtocol(λ,p_lambda,r,t_start,t_end) + + +# We plot the two trajectories + +# ```@example RateSystem +# fig = Figure(); axs = Axis(fig[1,1]) +# lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}") +# lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}") +# axislegend(axs,position=:rc,labelsize=10) +# fig +# ``` +fig = Figure(); axs = Axis(fig[1,1]) +lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}") +lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}") +axislegend(axs,position=:rc,labelsize=10) +fig From 260e18f41ae08e2082e427cb35861010cb24998d Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 30 Jun 2025 22:50:41 +0200 Subject: [PATCH 31/76] correction in MaierStein example --- docs/src/examples/gMAM_Maierstein.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/examples/gMAM_Maierstein.md b/docs/src/examples/gMAM_Maierstein.md index f0b6b7c5..705b85d8 100644 --- a/docs/src/examples/gMAM_Maierstein.md +++ b/docs/src/examples/gMAM_Maierstein.md @@ -18,7 +18,7 @@ font = (; nothing #hide ```` -Let us explore the features of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl) with Maier-Stein model. +Let us explore the features of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl) with the Maier-Stein model. ## Maier-stein model From 64d775c5656f7a074b0a1486f3ad813d5533f5e8 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 22 Jul 2025 09:53:17 +0200 Subject: [PATCH 32/76] small corrections in Test and RateSystemDraft --- test/ratesystem/RateSystemDraft2.jl | 6 +++--- test/ratesystem/RateSystemTestDraft2.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/ratesystem/RateSystemDraft2.jl b/test/ratesystem/RateSystemDraft2.jl index 838b24f7..60ad7229 100644 --- a/test/ratesystem/RateSystemDraft2.jl +++ b/test/ratesystem/RateSystemDraft2.jl @@ -37,8 +37,8 @@ function RateSystem(auto_sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t kwargs...) # we wish to return a continuous time dynamical system with modified drift field - f(u,p,t) = modified_drift(u,p,t,auto_sys,rp.λ,rp.t_start,rp.t_end,rp.r;kwargs...); - prob = remake(auto_sys.integ.sol.prob;f,p=rp.p_lambda,tspan=(t0,Inf)); - nonauto_sys = CoupledODEs(prob,auto_sys.diffeq); + f(u,p,t) = modified_drift(u,p,t,auto_sys,rp.λ,rp.t_start,rp.t_end,rp.r;kwargs...) + prob = remake(auto_sys.integ.sol.prob;f,p=rp.p_lambda,tspan=(t0,Inf)) + nonauto_sys = CoupledODEs(prob,auto_sys.diffeq) return nonauto_sys end diff --git a/test/ratesystem/RateSystemTestDraft2.jl b/test/ratesystem/RateSystemTestDraft2.jl index 42921043..d97ce6ab 100644 --- a/test/ratesystem/RateSystemTestDraft2.jl +++ b/test/ratesystem/RateSystemTestDraft2.jl @@ -38,7 +38,7 @@ t_end = Inf # end time of non-autonomous part # And the initial time of the system t0 = -10. -rp = RateProtocol(λ,p_lambda,r,t_start,t_end) +rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end) nonauto_sys = RateSystem(auto_sys,rp,t0) T = 20. # final simulation time From cdd41f4b727daa9077ec29eb33ccbe55a16abe54 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 22 Jul 2025 10:19:29 +0200 Subject: [PATCH 33/76] addition to RateSystem example --- docs/src/examples/RateSystem.md | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index 7f06564f..089bc571 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -1,6 +1,3 @@ -```@meta -EditURL = "../../../examples/RateSystem.jl" -``` # RateSystem ```@example RateSystem @@ -15,17 +12,18 @@ font = (; ); ``` -Let us explore an example of the RateSystem setting of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl). +Let us explore an example of how to use the RateSystem setting of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl). ## Prototypical model for R-tipping, with critical rate r = 4/3 -The following simple one-dimensional model with one attractor is given by the following ordinary differential equations: +We first consider the following simple one-dimensional autonomous system with one attractor, given by the ordinary differential equation: ```math \begin{aligned} \dot{x} &= (x+\lambda)^2 - 1 \end{aligned} ``` -The parameter ``\lambda`` shifts the location of the extrema of the drift field. +The parameter ``\lambda`` shifts the location of the extrema of the drift field. +We implement this system as follows: ```@example RateSystem function f(u,p,t) # out-of-place @@ -34,6 +32,7 @@ function f(u,p,t) # out-of-place dx = (x+λ)^2 - 1 return SVector{1}(dx) end + lambda = 0.0 p = [lambda] x0 = [-1.] @@ -42,7 +41,7 @@ auto_sys = CoupledODEs(f,x0,p) ## Non-autonomous case -Now, we define the time-dependent parameter function ``\lambda(t)`` to make the system non-autonomous and investigate the system's behaviour under parameter rampings. +Now, we want to explore the system in a non-autonomous setting with a time-dependent parameter function ``\lambda(rt)``. Choosing different values of the parameter ``r`` allows us to vary the speed of the parameter ramping. ```@example RateSystem function λ(p,t) @@ -56,21 +55,22 @@ We define the following parameters λ_max = 3. p_lambda = [λ_max] # parameter of the function lambda r = 4/3-0.02 # r just below critical rate -t_start = -Inf # start time of non-autonomous part -t_end = Inf # end time of non-autonomous part -# And the initial time of the system -t0 = -10. ``` We define the RateProtocol ```@example RateSystem + +t_start = -Inf # start time of non-autonomous part +t_end = Inf # end time of non-autonomous part + rp = RateProtocol(λ,p_lambda,r,t_start,t_end) ``` And use it to create the system with past and future autonomous systems and non-autonomous ramping in between: ```@example RateSystem +t0 = -10. # initial time of the system nonauto_sys = RateSystem(auto_sys,rp,t0) T = 20. # final simulation time @@ -90,4 +90,4 @@ fig ----- Author: Raphael Roemer -Date: 30 Jun 202r +Date: 30 Jun 2025 From cbe4abed12397c60084327cd956754fd04bdbbb4 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 22 Jul 2025 10:27:43 +0200 Subject: [PATCH 34/76] additions to RateSystem example --- docs/src/examples/RateSystem.md | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index 089bc571..fbee2555 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -41,26 +41,25 @@ auto_sys = CoupledODEs(f,x0,p) ## Non-autonomous case -Now, we want to explore the system in a non-autonomous setting with a time-dependent parameter function ``\lambda(rt)``. Choosing different values of the parameter ``r`` allows us to vary the speed of the parameter ramping. +Now, we want to explore a non-autonomous version of the system. +We consider a setting where in the past and in the future the system is autnonomous and in between there is a non-autonomous period ``[t_start, t_end]`` with a time-dependent parameter ramping given by the function ``\lambda(rt)``. Choosing different values of the parameter ``r`` allows us to vary the speed of the parameter ramping. +We start by defining the function ``\lambda(t)``: ```@example RateSystem function λ(p,t) λ_max = p[1] lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) return SVector{1}(lambda) end -``` -We define the following parameters -```@example RateSystem + λ_max = 3. p_lambda = [λ_max] # parameter of the function lambda -r = 4/3-0.02 # r just below critical rate ``` -We define the RateProtocol +Now, we define the RateProtocol that describes the non-autonomous period: ```@example RateSystem - +r = 4/3-0.02 # r just below critical rate t_start = -Inf # start time of non-autonomous part t_end = Inf # end time of non-autonomous part @@ -73,7 +72,7 @@ And use it to create the system with past and future autonomous systems and non- t0 = -10. # initial time of the system nonauto_sys = RateSystem(auto_sys,rp,t0) -T = 20. # final simulation time +T = 20. # final simulation time auto_traj = trajectory(auto_sys,T,x0) nonauto_traj = trajectory(nonauto_sys,T,x0) ``` From 6ab08ec89367e6da4f2db1bd71a6fb1528ae7de8 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 22 Jul 2025 10:30:14 +0200 Subject: [PATCH 35/76] additions to RateSystem example --- docs/src/examples/RateSystem.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index fbee2555..b3437083 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -66,12 +66,15 @@ t_end = Inf # end time of non-autonomous part rp = RateProtocol(λ,p_lambda,r,t_start,t_end) ``` -And use it to create the system with past and future autonomous systems and non-autonomous ramping in between: +Now, we set up the combined system with autonomous past and future and non-autonomous ramping in between: ```@example RateSystem t0 = -10. # initial time of the system nonauto_sys = RateSystem(auto_sys,rp,t0) +``` +We can compute trajectories of this new system in the familiar way: +```@example RateSystem T = 20. # final simulation time auto_traj = trajectory(auto_sys,T,x0) nonauto_traj = trajectory(nonauto_sys,T,x0) From c47fdd430ac149c3c4d9f392991d03104a9f7266 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 22 Jul 2025 10:31:43 +0200 Subject: [PATCH 36/76] correction in RateSystem example --- docs/src/examples/RateSystem.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index b3437083..586bc169 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -63,7 +63,7 @@ r = 4/3-0.02 # r just below critical rate t_start = -Inf # start time of non-autonomous part t_end = Inf # end time of non-autonomous part -rp = RateProtocol(λ,p_lambda,r,t_start,t_end) +rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end) ``` Now, we set up the combined system with autonomous past and future and non-autonomous ramping in between: From 5266412abe062072648f7c588bf1d64041af7ffb Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 22 Jul 2025 10:35:41 +0200 Subject: [PATCH 37/76] correction in RateSystem.jl --- examples/RateSystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/RateSystem.jl b/examples/RateSystem.jl index 4359119f..7ecad01b 100644 --- a/examples/RateSystem.jl +++ b/examples/RateSystem.jl @@ -77,7 +77,7 @@ t0 = -10. # ```@example RateSystem # rp = RateProtocol(λ,p_lambda,r,t_start,t_end) # ``` -rp = RateProtocol(λ,p_lambda,r,t_start,t_end) +rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end) # We plot the two trajectories From 12ab98f931941d1b662d9ebbe25dcad1f156a700 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Wed, 23 Jul 2025 18:39:04 +0200 Subject: [PATCH 38/76] added RateSystem.md to pages.jl --- docs/pages.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/pages.jl b/docs/pages.jl index 1b7c3571..01239674 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -8,6 +8,7 @@ pages = [ "sgMAM for the Kerr Parametric Oscillator" => "examples/sgMAM_KPO.md", "Transition Path Theory using finite element method" => "examples/duffing_TPT.md", "Minimal action method as an Optimal Control problem" => "examples/OC_mam.md", + "Setting up a RateSystem" => "examples/RateSystem.md" ], "Manual" => Any[ "Define a CoupledSDEs system" => "man/CoupledSDEs.md", From e9f90c8aa90a3ce7e6559c0743118333dcf84008 Mon Sep 17 00:00:00 2001 From: reykboerner Date: Thu, 24 Jul 2025 17:04:38 +0200 Subject: [PATCH 39/76] moved RateSystem source code to src --- src/CriticalTransitions.jl | 4 ++-- .../RateSystemDraft2.jl => src/r_tipping/RateSystem.jl | 0 2 files changed, 2 insertions(+), 2 deletions(-) rename test/ratesystem/RateSystemDraft2.jl => src/r_tipping/RateSystem.jl (100%) diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 8b8dc52b..f6b66e1b 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -62,11 +62,11 @@ include("largedeviations/geometric_min_action_method.jl") include("largedeviations/sgMAM.jl") include("largedeviations/string_method.jl") +include("r_tipping/RateSystem.jl") + include("../systems/CTLibrary.jl") using .CTLibrary -include("../test/ratesystem/RateSystemDraft2.jl") - # Core types export CoupledSDEs, CoupledODEs, noise_process, covariance_matrix, diffusion_matrix export drift, div_drift, solver diff --git a/test/ratesystem/RateSystemDraft2.jl b/src/r_tipping/RateSystem.jl similarity index 100% rename from test/ratesystem/RateSystemDraft2.jl rename to src/r_tipping/RateSystem.jl From 14a97d606c32c262f9927e1590b531091cbc5479 Mon Sep 17 00:00:00 2001 From: reykboerner Date: Thu, 24 Jul 2025 17:14:30 +0200 Subject: [PATCH 40/76] applied formatter, fixed typo, disabled spell check until PR review --- .github/workflows/SpellCheck.yml | 1 - docs/src/examples/Langevin_MCMC.md | 2 +- examples/RateSystem.jl | 45 +++++++++++++------- src/r_tipping/RateSystem.jl | 55 +++++++++++++++++-------- test/ratesystem/RateSystemDraft1.jl | 20 +++++---- test/ratesystem/RateSystemTestDraft1.jl | 19 +++++---- test/ratesystem/RateSystemTestDraft2.jl | 48 +++++++++++++-------- 7 files changed, 120 insertions(+), 70 deletions(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 94597ca3..17be9084 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -5,7 +5,6 @@ on: types: - opened - reopened - - synchronize - ready_for_review jobs: diff --git a/docs/src/examples/Langevin_MCMC.md b/docs/src/examples/Langevin_MCMC.md index 0f1f61d0..0468b298 100644 --- a/docs/src/examples/Langevin_MCMC.md +++ b/docs/src/examples/Langevin_MCMC.md @@ -13,7 +13,7 @@ The LMCMC ## FitzHugh Nagumo System (FHN) -As an example, we consider a gradient 2D system that is driven by a two dimensional Wiener rpocess, the FitzHugh Nagumo System (FHN). +As an example, we consider a gradient 2D system that is driven by a two dimensional Wiener process, the FitzHugh Nagumo System (FHN). The system's deterministic dynamics are given by: diff --git a/examples/RateSystem.jl b/examples/RateSystem.jl index 7ecad01b..bac867a9 100644 --- a/examples/RateSystem.jl +++ b/examples/RateSystem.jl @@ -7,8 +7,10 @@ using ModelingToolkit using CairoMakie using CairoMakie.Makie.MathTeXEngine: get_font font = (; - regular=get_font(:regular), bold=get_font(:bold), - italic=get_font(:italic), bold_italic=get_font(:bolditalic) + regular=get_font(:regular), + bold=get_font(:bold), + italic=get_font(:italic), + bold_italic=get_font(:bolditalic), ); # Let us explore an example of the RateSystem setting of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl). @@ -23,16 +25,16 @@ font = (; # ``` # The parameter ``\lambda`` shifts the location of the extrema of the drift field. -function f(u,p,t) # out-of-place +function f(u, p, t) # out-of-place x = u[1] λ = p[1] dx = (x+λ)^2 - 1 return SVector{1}(dx) end -lambda = 0.0 +lambda = 0.0 p = [lambda] -x0 = [-1.] -auto_sys = CoupledODEs(f,x0,p) +x0 = [-1.0] +auto_sys = CoupledODEs(f, x0, p) ## Non-autonomous case @@ -46,7 +48,7 @@ auto_sys = CoupledODEs(f,x0,p) # end # ``` -function λ(p,t) +function λ(p, t) λ_max = p[1] lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) return SVector{1}(lambda) @@ -63,22 +65,20 @@ end # t0 = -10. # ``` -λ_max = 3. +λ_max = 3.0 p_lambda = [λ_max] # parameter of the function lambda r = 4/3-0.02 # r just below critical rate t_start = -Inf # start time of non-autonomous part t_end = Inf # end time of non-autonomous part # And the initial time of the system -t0 = -10. - +t0 = -10.0 # We define the RateProtocol # ```@example RateSystem # rp = RateProtocol(λ,p_lambda,r,t_start,t_end) # ``` -rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end) - +rp = CriticalTransitions.RateProtocol(λ, p_lambda, r, t_start, t_end) # We plot the two trajectories @@ -89,8 +89,21 @@ rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end) # axislegend(axs,position=:rc,labelsize=10) # fig # ``` -fig = Figure(); axs = Axis(fig[1,1]) -lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}") -lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}") -axislegend(axs,position=:rc,labelsize=10) +fig = Figure(); +axs = Axis(fig[1, 1]) +lines!( + axs, + t0 .+ auto_traj[2], + auto_traj[1][:, 1]; + linewidth=2, + label=L"\text{Autonomous system}", +) +lines!( + axs, + nonauto_traj[2], + nonauto_traj[1][:, 1]; + linewidth=2, + label=L"\text{Nonautonomous system}", +) +axislegend(axs; position=:rc, labelsize=10) fig diff --git a/src/r_tipping/RateSystem.jl b/src/r_tipping/RateSystem.jl index 60ad7229..6d28d9a8 100644 --- a/src/r_tipping/RateSystem.jl +++ b/src/r_tipping/RateSystem.jl @@ -7,38 +7,57 @@ # Then we give back the ContinuousTimeDynamicalSystem with the parameter # changing according to the rate protocol -mutable struct RateProtocol - λ::Function - p_lambda::Vector - r::Float64 - t_start::Float64 +mutable struct RateProtocol + λ::Function + p_lambda::Vector + r::Float64 + t_start::Float64 t_end::Float64 end # convenience functions -RateProtocol(λ::Function,p_lambda::Vector,r::Float64)=RateProtocol(λ,p_lambda,r,-Inf,Inf) -RateProtocol(λ::Function,r::Float64)=RateProtocol(λ,[],r,-Inf,Inf) +function RateProtocol(λ::Function, p_lambda::Vector, r::Float64) + RateProtocol(λ, p_lambda, r, -Inf, Inf) +end +RateProtocol(λ::Function, r::Float64)=RateProtocol(λ, [], r, -Inf, Inf) #RateProtocol(λ::Function,p_lambda::Vector,r::Float64,t_start::Float64)=RateProtocol(λ,p_lambda,r,t_start,Inf) #RateProtocol(λ::Function,r::Float64,t_start::Float64)=RateProtocol(λ,[],r,t_start,Inf) -function modified_drift(u,p,t,ds::ContinuousTimeDynamicalSystem,λ::Function,t_start::Float64,t_end::Float64,r::Float64; - kwargs...) - - if t_start > t_end +function modified_drift( + u, + p, + t, + ds::ContinuousTimeDynamicalSystem, + λ::Function, + t_start::Float64, + t_end::Float64, + r::Float64; + kwargs..., +) + if t_start > t_end error("Please ensure that t_start ≤ t_end.") end - p̃ = r*t ≤ t_start ? λ(p,t_start;kwargs...) : t_start < r*t < t_end ? λ(p,r*t;kwargs...) : λ(p,t_end;kwargs...); # the value(s) of λ(rt) - return ds.integ.f(u,p̃,t) + p̃ = if r*t ≤ t_start + λ(p, t_start; kwargs...) + elseif t_start < r*t < t_end + λ(p, r*t; kwargs...) # the value(s) of λ(rt) + else + λ(p, t_end; kwargs...) # the value(s) of λ(rt) + end; # the value(s) of λ(rt) + return ds.integ.f(u, p̃, t) end; -function RateSystem(auto_sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t0::Float64; - kwargs...) +function RateSystem( + auto_sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t0::Float64; kwargs... +) # we wish to return a continuous time dynamical system with modified drift field - f(u,p,t) = modified_drift(u,p,t,auto_sys,rp.λ,rp.t_start,rp.t_end,rp.r;kwargs...) - prob = remake(auto_sys.integ.sol.prob;f,p=rp.p_lambda,tspan=(t0,Inf)) - nonauto_sys = CoupledODEs(prob,auto_sys.diffeq) + f(u, p, t) = modified_drift( + u, p, t, auto_sys, rp.λ, rp.t_start, rp.t_end, rp.r; kwargs... + ) + prob = remake(auto_sys.integ.sol.prob; f, p=rp.p_lambda, tspan=(t0, Inf)) + nonauto_sys = CoupledODEs(prob, auto_sys.diffeq) return nonauto_sys end diff --git a/test/ratesystem/RateSystemDraft1.jl b/test/ratesystem/RateSystemDraft1.jl index 019952ad..3cc8dab2 100644 --- a/test/ratesystem/RateSystemDraft1.jl +++ b/test/ratesystem/RateSystemDraft1.jl @@ -4,23 +4,27 @@ # we write -function RateSystem(tni,tnf,f,λ,p_λ,initvals) - func(u,p,t) = combined_system(u,t,tni,tnf,f,λ,p_λ); - return CoupledODEs(func, initvals, Float64[], t0=tni) +function RateSystem(tni, tnf, f, λ, p_λ, initvals) + func(u, p, t) = combined_system(u, t, tni, tnf, f, λ, p_λ); + return CoupledODEs(func, initvals, Float64[]; t0=tni) end -function combined_system(u,t,tni,tnf,f,λ,p_λ) - lambda = t < tni ? λ(u,p_λ,tni) : tni <= t <= tnf ? λ(u,p_λ,t) : λ(u,p_λ,tnf) - return f(u,lambda,t) +function combined_system(u, t, tni, tnf, f, λ, p_λ) + lambda = if t < tni + λ(u, p_λ, tni) + elseif tni <= t <= tnf + λ(u, p_λ, t) + else + λ(u, p_λ, tnf) + end + return f(u, lambda, t) end; - ############### # user writes # ...moved to test/RateSystem.jl (Reyk) - # further ideas # function futureSyst(RateSyst) # Canard Trajectories diff --git a/test/ratesystem/RateSystemTestDraft1.jl b/test/ratesystem/RateSystemTestDraft1.jl index 4fb306f4..ca5dc60e 100644 --- a/test/ratesystem/RateSystemTestDraft1.jl +++ b/test/ratesystem/RateSystemTestDraft1.jl @@ -1,22 +1,25 @@ using Plots -function f(u::SVector{1, Float64},p::SVector{1, Float64},t::Float64) +function f(u::SVector{1,Float64}, p::SVector{1,Float64}, t::Float64) x = u[1] λ = p[1] dx = (x+λ)^2 - 1 return SVector{1}(dx) end; -function λ(p::Vector{Float64},t::Float64) - r,λ_max = p - lambda = (λ_max/2)*(tanh(λ_max*r*t/2) +1) +function λ(p::Vector{Float64}, t::Float64) + r, λ_max = p + lambda = (λ_max/2)*(tanh(λ_max*r*t/2) + 1) return SVector{1}(lambda) end; -tni=-10.; tnf=10.; p_λ = [1.0,3.0]; initvals = [-4.]; -sys = RateSystem(tni,tnf,f,λ,p_λ,initvals) +tni=-10.0; +tnf=10.0; +p_λ = [1.0, 3.0]; +initvals = [-4.0]; +sys = RateSystem(tni, tnf, f, λ, p_λ, initvals) -traj=trajectory(sys,40.,t0=-20.) +traj=trajectory(sys, 40.0; t0=-20.0) trajic=[traj[1][i][1] for i in 1:length(traj[1])] -plot(collect(traj[2]),trajic) +plot(collect(traj[2]), trajic) diff --git a/test/ratesystem/RateSystemTestDraft2.jl b/test/ratesystem/RateSystemTestDraft2.jl index d97ce6ab..f15fe8f5 100644 --- a/test/ratesystem/RateSystemTestDraft2.jl +++ b/test/ratesystem/RateSystemTestDraft2.jl @@ -3,50 +3,62 @@ ############################################################################################################## # autonomous system -function f(u,p,t) +function f(u, p, t) x = u[1] λ = p[1] dx = (x+λ)^2 - 1 return SVector{1}(dx) end -lambda = 0. +lambda = 0.0 p = [lambda] -x0 = [-1.] +x0 = [-1.0] Δt = 1.e-3 #diffeq = (alg = AutoVern9(Rodas5(autodiff=false)),abstol=1.e-16,reltol=1.e-16) #auto_sys = CoupledODEs(f,x0,p;diffeq) -auto_sys = CoupledODEs(f,x0,p) - +auto_sys = CoupledODEs(f, x0, p) # now we make the system non-autonomous # We define a time-dependent parameter function -function λ(p,t) +function λ(p, t) λ_max = p[1] lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) return SVector{1}(lambda) end # we set the parameters for RateProtocol -λ_max = 3. +λ_max = 3.0 p_lambda = [λ_max] r = 4/3-0.02 # r just below critical rate t_start = -Inf # start time of non-autonomous part t_end = Inf # end time of non-autonomous part # And the initial time of the system -t0 = -10. +t0 = -10.0 -rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end) -nonauto_sys = RateSystem(auto_sys,rp,t0) +rp = CriticalTransitions.RateProtocol(λ, p_lambda, r, t_start, t_end) +nonauto_sys = RateSystem(auto_sys, rp, t0) -T = 20. # final simulation time -auto_traj = trajectory(auto_sys,T,x0) -nonauto_traj = trajectory(nonauto_sys,T,x0) +T = 20.0 # final simulation time +auto_traj = trajectory(auto_sys, T, x0) +nonauto_traj = trajectory(nonauto_sys, T, x0) -fig = Figure(); axs = Axis(fig[1,1]) -lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}") -lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}") -axislegend(axs,position=:rc,labelsize=10) -fig \ No newline at end of file +fig = Figure(); +axs = Axis(fig[1, 1]) +lines!( + axs, + t0 .+ auto_traj[2], + auto_traj[1][:, 1]; + linewidth=2, + label=L"\text{Autonomous system}", +) +lines!( + axs, + nonauto_traj[2], + nonauto_traj[1][:, 1]; + linewidth=2, + label=L"\text{Nonautonomous system}", +) +axislegend(axs; position=:rc, labelsize=10) +fig From f8d73effe0ed45d92efca2c9b925b6e1e65e0c8d Mon Sep 17 00:00:00 2001 From: reykboerner Date: Fri, 25 Jul 2025 15:05:31 +0200 Subject: [PATCH 41/76] remove CairoMakie dep --- Project.toml | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 3cd7362f..6960e4c4 100644 --- a/Project.toml +++ b/Project.toml @@ -5,7 +5,6 @@ repo = "https://github.com/juliadynamics/CriticalTransitions.jl.git" version = "0.6.0" [deps] -CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4" @@ -47,16 +46,16 @@ DynamicalSystemsBase = "3.13" Format = "1" Interpolations = "0.15.1, 0.16" IntervalArithmetic = "0.22.22" +LinearAlgebra = "1.10" LinearSolve = "2.37.0, 3" Optimization = "4" OptimizationOptimisers = "0.3" +Printf = "1.10" ProgressMeter = "1.10.0" +Random = "1.10" Reexport = "1.2.2" SciMLBase = "2.70" SparseArrays = "1.10" -Random = "1.10" -LinearAlgebra = "1.10" -Printf = "1.10" StateSpaceSets = "2.3.0" StaticArrays = "1.9.8" Statistics = ">=1.9" From cf5afb2c8cdd874e39d50af2cff02a8c01bfbbee Mon Sep 17 00:00:00 2001 From: reykboerner Date: Fri, 25 Jul 2025 15:35:48 +0200 Subject: [PATCH 42/76] added RateSystem test --- test/r_tipping/RateSystem.jl | 0 test/ratesystem/RateSystem.jl | 1 + test/runtests.jl | 4 ++++ 3 files changed, 5 insertions(+) create mode 100644 test/r_tipping/RateSystem.jl create mode 100644 test/ratesystem/RateSystem.jl diff --git a/test/r_tipping/RateSystem.jl b/test/r_tipping/RateSystem.jl new file mode 100644 index 00000000..e69de29b diff --git a/test/ratesystem/RateSystem.jl b/test/ratesystem/RateSystem.jl new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/test/ratesystem/RateSystem.jl @@ -0,0 +1 @@ + diff --git a/test/runtests.jl b/test/runtests.jl index 0e921e57..34a47f64 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,6 +44,10 @@ end include("trajectories/transition.jl") end +@testset "R-tipping" begin + include("r_tipping/RateSystem.jl") +end + @testset "Extensions" begin @testset "ChaosToolsExt" begin include("ext/ChaosToolsExt.jl") From e7f8580344d0b35aa438428ebca067daa81f3940 Mon Sep 17 00:00:00 2001 From: reykboerner Date: Fri, 25 Jul 2025 15:49:34 +0200 Subject: [PATCH 43/76] small fix --- docs/src/examples/RateSystem.md | 4 ++-- test/r_tipping/RateSystem.jl | 41 +++++++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+), 2 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index 586bc169..279b73b6 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -2,8 +2,8 @@ ```@example RateSystem using CriticalTransitions -using DifferentialEquations -using ModelingToolkit +#using DifferentialEquations +#using ModelingToolkit using CairoMakie using CairoMakie.Makie.MathTeXEngine: get_font font = (; diff --git a/test/r_tipping/RateSystem.jl b/test/r_tipping/RateSystem.jl index e69de29b..621d6693 100644 --- a/test/r_tipping/RateSystem.jl +++ b/test/r_tipping/RateSystem.jl @@ -0,0 +1,41 @@ +using CriticalTransitions +using Test + +@testset "RateSystem" begin + function f(u,p,t) # out-of-place + x = u[1] + λ = p[1] + dx = (x+λ)^2 - 1 + return SVector{1}(dx) + end + + lambda = 0.0 + p = [lambda] + x0 = [-1.] + auto_sys = CoupledODEs(f,x0,p) + + function λ(p,t) + λ_max = p[1] + lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) + return SVector{1}(lambda) + end + + λ_max = 3. + p_lambda = [λ_max] # parameter of the function lambda + + r = 4/3-0.02 # r just below critical rate + t_start = -Inf # start time of non-autonomous part + t_end = Inf # end time of non-autonomous part + + rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end) + + t0 = -10. # initial time of the system + nonauto_sys = RateSystem(auto_sys,rp,t0) + + T = 20. # final simulation time + auto_traj = trajectory(auto_sys,T,x0) + nonauto_traj = trajectory(nonauto_sys,T,x0) + + @test isapprox(auto_traj[1][end,1], -1; atol=1e-4) + @test isapprox(nonauto_traj[1][end,1], -4; atol=1e-4) +end \ No newline at end of file From 063b84b1399ea3bd3d2b125e1a26877a417d78d9 Mon Sep 17 00:00:00 2001 From: reykboerner Date: Fri, 25 Jul 2025 16:11:05 +0200 Subject: [PATCH 44/76] applied Formatter --- test/r_tipping/RateSystem.jl | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/test/r_tipping/RateSystem.jl b/test/r_tipping/RateSystem.jl index 621d6693..f056b938 100644 --- a/test/r_tipping/RateSystem.jl +++ b/test/r_tipping/RateSystem.jl @@ -2,40 +2,40 @@ using CriticalTransitions using Test @testset "RateSystem" begin - function f(u,p,t) # out-of-place + function f(u, p, t) # out-of-place x = u[1] λ = p[1] dx = (x+λ)^2 - 1 return SVector{1}(dx) end - lambda = 0.0 + lambda = 0.0 p = [lambda] - x0 = [-1.] - auto_sys = CoupledODEs(f,x0,p) + x0 = [-1.0] + auto_sys = CoupledODEs(f, x0, p) - function λ(p,t) + function λ(p, t) λ_max = p[1] lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) return SVector{1}(lambda) end - λ_max = 3. + λ_max = 3.0 p_lambda = [λ_max] # parameter of the function lambda r = 4/3-0.02 # r just below critical rate t_start = -Inf # start time of non-autonomous part t_end = Inf # end time of non-autonomous part - rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end) + rp = CriticalTransitions.RateProtocol(λ, p_lambda, r, t_start, t_end) - t0 = -10. # initial time of the system - nonauto_sys = RateSystem(auto_sys,rp,t0) + t0 = -10.0 # initial time of the system + nonauto_sys = RateSystem(auto_sys, rp, t0) - T = 20. # final simulation time - auto_traj = trajectory(auto_sys,T,x0) - nonauto_traj = trajectory(nonauto_sys,T,x0) + T = 20.0 # final simulation time + auto_traj = trajectory(auto_sys, T, x0) + nonauto_traj = trajectory(nonauto_sys, T, x0) - @test isapprox(auto_traj[1][end,1], -1; atol=1e-4) - @test isapprox(nonauto_traj[1][end,1], -4; atol=1e-4) -end \ No newline at end of file + @test isapprox(auto_traj[1][end, 1], -1; atol=1e-4) + @test isapprox(nonauto_traj[1][end, 1], -4; atol=1e-4) +end From 1feff1ae24b5ba666fb19734fc472e817f79b906 Mon Sep 17 00:00:00 2001 From: reykboerner Date: Fri, 25 Jul 2025 16:30:05 +0200 Subject: [PATCH 45/76] small edits in docs and added docstring drafts --- docs/src/examples/RateSystem.md | 8 +++----- src/r_tipping/RateSystem.jl | 19 ++++++++++++++++++- 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index 279b73b6..4449dd84 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -1,9 +1,8 @@ -# RateSystem +# Setting up a `RateSystem` +Let us explore an example of how to construct a `RateSystem` from an autonomous dynamical system (e.g. a `CoupledODEs`) and a time-dependent forcing protocol called `RateProtocol`. ```@example RateSystem using CriticalTransitions -#using DifferentialEquations -#using ModelingToolkit using CairoMakie using CairoMakie.Makie.MathTeXEngine: get_font font = (; @@ -12,8 +11,6 @@ font = (; ); ``` -Let us explore an example of how to use the RateSystem setting of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl). - ## Prototypical model for R-tipping, with critical rate r = 4/3 We first consider the following simple one-dimensional autonomous system with one attractor, given by the ordinary differential equation: @@ -92,4 +89,5 @@ fig ----- Author: Raphael Roemer + Date: 30 Jun 2025 diff --git a/src/r_tipping/RateSystem.jl b/src/r_tipping/RateSystem.jl index 6d28d9a8..f3029580 100644 --- a/src/r_tipping/RateSystem.jl +++ b/src/r_tipping/RateSystem.jl @@ -7,6 +7,16 @@ # Then we give back the ContinuousTimeDynamicalSystem with the parameter # changing according to the rate protocol +""" + RateProtocol + +Time-dependent forcing protocol specified by the following fields: +- `λ::Function`: forcing function +- `p_lambda::Vector`: parameters of forcing function +- `r::Float64`: rate parameter +- `t_start::Float64`: start time of protocol +- `t_end::Float64`: end time of protocol +""" mutable struct RateProtocol λ::Function p_lambda::Vector @@ -49,8 +59,15 @@ function modified_drift( return ds.integ.f(u, p̃, t) end; +""" + RateSystem(sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t0=0.0; kwargs...) + +Applies a time-dependent [`RateProtocol`](@def) to a given autonomous dynamical system +`sys`, turning it into a non-autonomous dynamical system. Returns a [`CoupledODEs`](@ref) +with the explicit parameter time-dependence incorporated. +""" function RateSystem( - auto_sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t0::Float64; kwargs... + auto_sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t0=0.0; kwargs... ) # we wish to return a continuous time dynamical system with modified drift field From 4e0ab57692cc0780a5ff4fa93774c579f2a6d1ec Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 28 Jul 2025 17:55:47 +0200 Subject: [PATCH 46/76] deleted test/ratesystem/RateSystem.jl --- test/ratesystem/RateSystem.jl | 1 - 1 file changed, 1 deletion(-) delete mode 100644 test/ratesystem/RateSystem.jl diff --git a/test/ratesystem/RateSystem.jl b/test/ratesystem/RateSystem.jl deleted file mode 100644 index 8b137891..00000000 --- a/test/ratesystem/RateSystem.jl +++ /dev/null @@ -1 +0,0 @@ - From 41128cc9294452710ac39aa1a6159405bc03fc03 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 28 Jul 2025 18:02:55 +0200 Subject: [PATCH 47/76] deleted test/ratesystem/RateSystemDraft1.jl and test/ratesystem/RateSystemTestDraft1.jl --- test/ratesystem/RateSystemDraft1.jl | 32 ------------------------- test/ratesystem/RateSystemTestDraft1.jl | 25 ------------------- 2 files changed, 57 deletions(-) delete mode 100644 test/ratesystem/RateSystemDraft1.jl delete mode 100644 test/ratesystem/RateSystemTestDraft1.jl diff --git a/test/ratesystem/RateSystemDraft1.jl b/test/ratesystem/RateSystemDraft1.jl deleted file mode 100644 index 3cc8dab2..00000000 --- a/test/ratesystem/RateSystemDraft1.jl +++ /dev/null @@ -1,32 +0,0 @@ -# | | | | -# t_i autonomous t_ni non-autonomous t_nf autonomous t_f -# | | | | - -# we write - -function RateSystem(tni, tnf, f, λ, p_λ, initvals) - func(u, p, t) = combined_system(u, t, tni, tnf, f, λ, p_λ); - return CoupledODEs(func, initvals, Float64[]; t0=tni) -end - -function combined_system(u, t, tni, tnf, f, λ, p_λ) - lambda = if t < tni - λ(u, p_λ, tni) - elseif tni <= t <= tnf - λ(u, p_λ, t) - else - λ(u, p_λ, tnf) - end - return f(u, lambda, t) -end; - -############### -# user writes - -# ...moved to test/RateSystem.jl (Reyk) - -# further ideas -# function futureSyst(RateSyst) -# Canard Trajectories - -# \dot(x) = f(x(t),λ(t)) diff --git a/test/ratesystem/RateSystemTestDraft1.jl b/test/ratesystem/RateSystemTestDraft1.jl deleted file mode 100644 index ca5dc60e..00000000 --- a/test/ratesystem/RateSystemTestDraft1.jl +++ /dev/null @@ -1,25 +0,0 @@ -using Plots - -function f(u::SVector{1,Float64}, p::SVector{1,Float64}, t::Float64) - x = u[1] - λ = p[1] - dx = (x+λ)^2 - 1 - return SVector{1}(dx) -end; - -function λ(p::Vector{Float64}, t::Float64) - r, λ_max = p - lambda = (λ_max/2)*(tanh(λ_max*r*t/2) + 1) - return SVector{1}(lambda) -end; - -tni=-10.0; -tnf=10.0; -p_λ = [1.0, 3.0]; -initvals = [-4.0]; -sys = RateSystem(tni, tnf, f, λ, p_λ, initvals) - -traj=trajectory(sys, 40.0; t0=-20.0) - -trajic=[traj[1][i][1] for i in 1:length(traj[1])] -plot(collect(traj[2]), trajic) From 2dcffc86ff3ff54f8fcac342e770e7ec79a0ba81 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 28 Jul 2025 18:08:54 +0200 Subject: [PATCH 48/76] deleted test/ratesystem --- test/ratesystem/RateSystemTestDraft2.jl | 64 ------------------------- 1 file changed, 64 deletions(-) delete mode 100644 test/ratesystem/RateSystemTestDraft2.jl diff --git a/test/ratesystem/RateSystemTestDraft2.jl b/test/ratesystem/RateSystemTestDraft2.jl deleted file mode 100644 index f15fe8f5..00000000 --- a/test/ratesystem/RateSystemTestDraft2.jl +++ /dev/null @@ -1,64 +0,0 @@ -############################################################################################################## -# test (prototypical model used for studying R-tipping, with critical rate r = 4/3) -############################################################################################################## - -# autonomous system -function f(u, p, t) - x = u[1] - λ = p[1] - dx = (x+λ)^2 - 1 - return SVector{1}(dx) -end - -lambda = 0.0 -p = [lambda] -x0 = [-1.0] - -Δt = 1.e-3 -#diffeq = (alg = AutoVern9(Rodas5(autodiff=false)),abstol=1.e-16,reltol=1.e-16) -#auto_sys = CoupledODEs(f,x0,p;diffeq) -auto_sys = CoupledODEs(f, x0, p) - -# now we make the system non-autonomous - -# We define a time-dependent parameter function -function λ(p, t) - λ_max = p[1] - lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) - return SVector{1}(lambda) -end - -# we set the parameters for RateProtocol -λ_max = 3.0 -p_lambda = [λ_max] -r = 4/3-0.02 # r just below critical rate -t_start = -Inf # start time of non-autonomous part -t_end = Inf # end time of non-autonomous part -# And the initial time of the system -t0 = -10.0 - -rp = CriticalTransitions.RateProtocol(λ, p_lambda, r, t_start, t_end) -nonauto_sys = RateSystem(auto_sys, rp, t0) - -T = 20.0 # final simulation time -auto_traj = trajectory(auto_sys, T, x0) -nonauto_traj = trajectory(nonauto_sys, T, x0) - -fig = Figure(); -axs = Axis(fig[1, 1]) -lines!( - axs, - t0 .+ auto_traj[2], - auto_traj[1][:, 1]; - linewidth=2, - label=L"\text{Autonomous system}", -) -lines!( - axs, - nonauto_traj[2], - nonauto_traj[1][:, 1]; - linewidth=2, - label=L"\text{Nonautonomous system}", -) -axislegend(axs; position=:rc, labelsize=10) -fig From 00f0ef2fd4bde0ab08e07a3e433113b0acc03506 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 28 Jul 2025 18:10:21 +0200 Subject: [PATCH 49/76] deleted examples/RateSystem.jl --- examples/RateSystem.jl | 109 ----------------------------------------- 1 file changed, 109 deletions(-) delete mode 100644 examples/RateSystem.jl diff --git a/examples/RateSystem.jl b/examples/RateSystem.jl deleted file mode 100644 index bac867a9..00000000 --- a/examples/RateSystem.jl +++ /dev/null @@ -1,109 +0,0 @@ - -# # Rate System Example - -using CriticalTransitions -using DifferentialEquations -using ModelingToolkit -using CairoMakie -using CairoMakie.Makie.MathTeXEngine: get_font -font = (; - regular=get_font(:regular), - bold=get_font(:bold), - italic=get_font(:italic), - bold_italic=get_font(:bolditalic), -); - -# Let us explore an example of the RateSystem setting of [CriticalTransitions.jl](https://github.com/JuliaDynamics/CriticalTransitions.jl). - -# ## Prototypical model for R-tipping, with critical rate r = 4/3 - -# The following simple one-dimensional model with one attractor is given by the following ordinary differential equations: -# ```math -# \begin{aligned} -# \dot{x} &= (x+\lambda)^2 - 1 -# \end{aligned} -# ``` -# The parameter ``\lambda`` shifts the location of the extrema of the drift field. - -function f(u, p, t) # out-of-place - x = u[1] - λ = p[1] - dx = (x+λ)^2 - 1 - return SVector{1}(dx) -end -lambda = 0.0 -p = [lambda] -x0 = [-1.0] -auto_sys = CoupledODEs(f, x0, p) - -## Non-autonomous case - -# Now, we define the time-dependent parameter function ``\lambda(t)`` to make the system non-autonomous and investigate the system's behaviour under parameter rampings. - -# ```@example RateSystem -# function λ(p,t) -# λ_max = p[1] -# lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) -# return SVector{1}(lambda) -# end -# ``` - -function λ(p, t) - λ_max = p[1] - lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) - return SVector{1}(lambda) -end - -# We define the following parameters -# ```@example RateSystem -# λ_max = 3. -# p_lambda = [λ_max] # parameter of the function lambda -# r = 4/3-0.02 # r just below critical rate -# t_start = -Inf # start time of non-autonomous part -# t_end = Inf # end time of non-autonomous part -# # And the initial time of the system -# t0 = -10. -# ``` - -λ_max = 3.0 -p_lambda = [λ_max] # parameter of the function lambda -r = 4/3-0.02 # r just below critical rate -t_start = -Inf # start time of non-autonomous part -t_end = Inf # end time of non-autonomous part -# And the initial time of the system -t0 = -10.0 - -# We define the RateProtocol - -# ```@example RateSystem -# rp = RateProtocol(λ,p_lambda,r,t_start,t_end) -# ``` -rp = CriticalTransitions.RateProtocol(λ, p_lambda, r, t_start, t_end) - -# We plot the two trajectories - -# ```@example RateSystem -# fig = Figure(); axs = Axis(fig[1,1]) -# lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}") -# lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}") -# axislegend(axs,position=:rc,labelsize=10) -# fig -# ``` -fig = Figure(); -axs = Axis(fig[1, 1]) -lines!( - axs, - t0 .+ auto_traj[2], - auto_traj[1][:, 1]; - linewidth=2, - label=L"\text{Autonomous system}", -) -lines!( - axs, - nonauto_traj[2], - nonauto_traj[1][:, 1]; - linewidth=2, - label=L"\text{Nonautonomous system}", -) -axislegend(axs; position=:rc, labelsize=10) -fig From 729ea5f9457870d4e8e15e3f1ef0fbf7372450b9 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 28 Jul 2025 18:15:47 +0200 Subject: [PATCH 50/76] expanded documentaation of RateSystem --- src/r_tipping/RateSystem.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/r_tipping/RateSystem.jl b/src/r_tipping/RateSystem.jl index f3029580..6c956873 100644 --- a/src/r_tipping/RateSystem.jl +++ b/src/r_tipping/RateSystem.jl @@ -65,6 +65,13 @@ end; Applies a time-dependent [`RateProtocol`](@def) to a given autonomous dynamical system `sys`, turning it into a non-autonomous dynamical system. Returns a [`CoupledODEs`](@ref) with the explicit parameter time-dependence incorporated. + +Then, the system is first autonomous, then non-autnonmous with the parameter shift given by the RateProtocol +and again autonomous at the end: + + | | | | +t_i autonomous t_start non-autonomous t_end autonomous t_f + | | | | """ function RateSystem( auto_sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t0=0.0; kwargs... From 648c863edb14ae8282733f466d36a1810e121590 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 28 Jul 2025 19:01:05 +0200 Subject: [PATCH 51/76] expanded documentation --- src/r_tipping/RateSystem.jl | 30 ++++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/src/r_tipping/RateSystem.jl b/src/r_tipping/RateSystem.jl index 6c956873..6fe6a9ee 100644 --- a/src/r_tipping/RateSystem.jl +++ b/src/r_tipping/RateSystem.jl @@ -11,11 +11,20 @@ RateProtocol Time-dependent forcing protocol specified by the following fields: -- `λ::Function`: forcing function +- `λ::Function`: forcing function of the form `λ(p, t_start; kwargs...)`` - `p_lambda::Vector`: parameters of forcing function - `r::Float64`: rate parameter - `t_start::Float64`: start time of protocol - `t_end::Float64`: end time of protocol + +If `t_start` and `t_end` are not provided, they are set to `t_start=-Inf`, and `t_end=Inf`. +If `p_lambda` is not provided, it is set to an empty array `[]`. + +The forcing protocol contains all the information that allows to take an ODE of the form +`dxₜ/dt = f(xₜ,λ)` +with `λ ∈ Rᵐ`` containing all the system parameters, and make it an ODE of the form +`dxₜ/dt = f(xₜ,λ(rt))`` +with `λ(t) ∈ Rᵐ`. """ mutable struct RateProtocol λ::Function @@ -60,21 +69,22 @@ function modified_drift( end; """ - RateSystem(sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t0=0.0; kwargs...) + apply_ramping!(sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs...) -Applies a time-dependent [`RateProtocol`](@def) to a given autonomous dynamical system +Applies a time-dependent [`RateProtocol`](@def) to a given autonomous deterministic dynamical system `sys`, turning it into a non-autonomous dynamical system. Returns a [`CoupledODEs`](@ref) with the explicit parameter time-dependence incorporated. -Then, the system is first autonomous, then non-autnonmous with the parameter shift given by the RateProtocol -and again autonomous at the end: +The returned [`CoupledODEs`](@ref) is autonomous from `t_0` to `t_start`, +then non-autnonmous from `t_start` to `t_end` with the parameter shift given by the [`RateProtocol`](@def), +and again autonomous from `t_end` to the end of the simulation: - | | | | -t_i autonomous t_start non-autonomous t_end autonomous t_f - | | | | + | | | | +`t_0` autonomous `t_start` non-autonomous `t_end` autonomous `∞` + | | | | """ -function RateSystem( - auto_sys::ContinuousTimeDynamicalSystem, rp::RateProtocol, t0=0.0; kwargs... +function apply_ramping!( + auto_sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs... ) # we wish to return a continuous time dynamical system with modified drift field From b147bebfd7f610639c7d3c2bc00bbdf8612b125a Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Mon, 28 Jul 2025 23:39:08 +0200 Subject: [PATCH 52/76] added plot of parameter shift in RateSystem documentation --- docs/src/examples/RateSystem.md | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index 4449dd84..085ec133 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -63,11 +63,21 @@ t_end = Inf # end time of non-autonomous part rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end) ``` +We plot the function ``λ(p_lambda,t)`` + +```@example RateSystem +figλ = Figure(); axsλ = Axis(figλ[1,1]) +lines!(axsλ,-10.:0.1:20.,λ(p_lambda,-10.:0.1:20.),linewidth=2,label=L"\text{λ(p_lambda,t)}") +axislegend(axsλ,position=:rc,labelsize=10) +figλ +``` + + Now, we set up the combined system with autonomous past and future and non-autonomous ramping in between: ```@example RateSystem t0 = -10. # initial time of the system -nonauto_sys = RateSystem(auto_sys,rp,t0) +nonauto_sys = apply_ramping!(auto_sys,rp,t0) ``` We can compute trajectories of this new system in the familiar way: From 79715afd4352f49a3e6fe42d2077b6c51a7a97eb Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 00:13:59 +0200 Subject: [PATCH 53/76] deleted NLPModelsIpopt from project.toml --- docs/Project.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 25abf93f..630682cd 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -11,9 +11,8 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656" DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71" -OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" From eb43377217dc31e7c6bc3430a69833cb703753a5 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 00:38:34 +0200 Subject: [PATCH 54/76] fixed error in Project.toml --- docs/Project.toml | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 630682cd..c15b0ddb 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,25 +1,3 @@ [deps] -Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7" -CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" -ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7" -Contour = "d38c429a-6771-53c6-b99e-75d170b6e991" CriticalTransitions = "251e6cd3-3112-48a5-99dd-66efcfd18334" -DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df" -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" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" -OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948" -Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" -OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" - -[compat] -Documenter = "^1.4.1" From 4cdbd7da960c9414d96a6d15faa27a4e5eea02fc Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 00:39:21 +0200 Subject: [PATCH 55/76] fixed quotation marks in quickstart.md --- docs/src/quickstart.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md index ebe410e8..c7748ba2 100644 --- a/docs/src/quickstart.md +++ b/docs/src/quickstart.md @@ -3,7 +3,7 @@ ## Installation Although the package is not yet registered, you can install it from Github via the Julia package manager: ```julia -using Pkg; Pkg.add('https://github.com/juliadynamics/CriticalTransitions.jl.git') +using Pkg; Pkg.add("https://github.com/juliadynamics/CriticalTransitions.jl.git") ``` The package is currently tested to be compatible with Julia versions `1.10` and `1.11`. From db71cb3e125df381b0a7e3686a8c5c757424cdde Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 01:00:22 +0200 Subject: [PATCH 56/76] removed the ! from apply_ramping --- docs/src/examples/RateSystem.md | 2 +- src/CriticalTransitions.jl | 2 +- src/r_tipping/RateSystem.jl | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index 085ec133..a428bb3f 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -67,7 +67,7 @@ We plot the function ``λ(p_lambda,t)`` ```@example RateSystem figλ = Figure(); axsλ = Axis(figλ[1,1]) -lines!(axsλ,-10.:0.1:20.,λ(p_lambda,-10.:0.1:20.),linewidth=2,label=L"\text{λ(p_lambda,t)}") +lines!(axsλ,-10.:0.1:20.,λ.(p_lambda,-10.:0.1:20.),linewidth=2,label=L"\lambda(p_{\lambda}, t)") axislegend(axsλ,position=:rc,labelsize=10) figλ ``` diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index f6b66e1b..12532a88 100644 --- a/src/CriticalTransitions.jl +++ b/src/CriticalTransitions.jl @@ -88,7 +88,7 @@ export find_boundary, huniform, dunion export committor, invariant_pdf, reactive_current, probability_reactive, probability_last_A, Langevin -export RateProtocol, RateSystem +export RateProtocol, apply_ramping # Error hint for extensions stubs function __init__() diff --git a/src/r_tipping/RateSystem.jl b/src/r_tipping/RateSystem.jl index 6fe6a9ee..aa6a3a09 100644 --- a/src/r_tipping/RateSystem.jl +++ b/src/r_tipping/RateSystem.jl @@ -69,7 +69,7 @@ function modified_drift( end; """ - apply_ramping!(sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs...) + apply_ramping(sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs...) Applies a time-dependent [`RateProtocol`](@def) to a given autonomous deterministic dynamical system `sys`, turning it into a non-autonomous dynamical system. Returns a [`CoupledODEs`](@ref) @@ -83,7 +83,7 @@ and again autonomous from `t_end` to the end of the simulation: `t_0` autonomous `t_start` non-autonomous `t_end` autonomous `∞` | | | | """ -function apply_ramping!( +function apply_ramping( auto_sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs... ) # we wish to return a continuous time dynamical system with modified drift field From 923849d640eac164a8e62201da1b2d83ca7f4aab Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 01:06:05 +0200 Subject: [PATCH 57/76] fixed error in documentation of RateSystem --- docs/src/examples/RateSystem.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index a428bb3f..edf6afdb 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -77,7 +77,7 @@ Now, we set up the combined system with autonomous past and future and non-auton ```@example RateSystem t0 = -10. # initial time of the system -nonauto_sys = apply_ramping!(auto_sys,rp,t0) +nonauto_sys = apply_ramping(auto_sys,rp,t0) ``` We can compute trajectories of this new system in the familiar way: From 8634b8d8c177e068191966d496d5cab72dee354c Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 01:17:35 +0200 Subject: [PATCH 58/76] resolving documentation issue with RateProtocol-plot --- docs/src/examples/RateSystem.md | 8 -------- 1 file changed, 8 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index edf6afdb..3bd89718 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -63,14 +63,6 @@ t_end = Inf # end time of non-autonomous part rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end) ``` -We plot the function ``λ(p_lambda,t)`` - -```@example RateSystem -figλ = Figure(); axsλ = Axis(figλ[1,1]) -lines!(axsλ,-10.:0.1:20.,λ.(p_lambda,-10.:0.1:20.),linewidth=2,label=L"\lambda(p_{\lambda}, t)") -axislegend(axsλ,position=:rc,labelsize=10) -figλ -``` Now, we set up the combined system with autonomous past and future and non-autonomous ramping in between: From 1c0999d06491c229a1005102b809f8d8eec223ed Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 01:22:03 +0200 Subject: [PATCH 59/76] improved documentation of apply_ramping --- src/r_tipping/RateSystem.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/r_tipping/RateSystem.jl b/src/r_tipping/RateSystem.jl index aa6a3a09..ceae2de7 100644 --- a/src/r_tipping/RateSystem.jl +++ b/src/r_tipping/RateSystem.jl @@ -79,9 +79,7 @@ The returned [`CoupledODEs`](@ref) is autonomous from `t_0` to `t_start`, then non-autnonmous from `t_start` to `t_end` with the parameter shift given by the [`RateProtocol`](@def), and again autonomous from `t_end` to the end of the simulation: - | | | | `t_0` autonomous `t_start` non-autonomous `t_end` autonomous `∞` - | | | | """ function apply_ramping( auto_sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs... From 44f11fb206c55d0e5c3dacaf20e058638dc7fc69 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 01:27:02 +0200 Subject: [PATCH 60/76] correction in tes of RaateSystem --- test/r_tipping/RateSystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/r_tipping/RateSystem.jl b/test/r_tipping/RateSystem.jl index f056b938..be24c286 100644 --- a/test/r_tipping/RateSystem.jl +++ b/test/r_tipping/RateSystem.jl @@ -30,7 +30,7 @@ using Test rp = CriticalTransitions.RateProtocol(λ, p_lambda, r, t_start, t_end) t0 = -10.0 # initial time of the system - nonauto_sys = RateSystem(auto_sys, rp, t0) + nonauto_sys = apply_ramping(auto_sys, rp, t0) T = 20.0 # final simulation time auto_traj = trajectory(auto_sys, T, x0) From 722cc748364db352a17161d67eefc0fdeba95571 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 09:42:18 +0200 Subject: [PATCH 61/76] corrected typo in RateProtocol docstring --- src/r_tipping/RateSystem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/r_tipping/RateSystem.jl b/src/r_tipping/RateSystem.jl index ceae2de7..f1f668c7 100644 --- a/src/r_tipping/RateSystem.jl +++ b/src/r_tipping/RateSystem.jl @@ -22,8 +22,8 @@ If `p_lambda` is not provided, it is set to an empty array `[]`. The forcing protocol contains all the information that allows to take an ODE of the form `dxₜ/dt = f(xₜ,λ)` -with `λ ∈ Rᵐ`` containing all the system parameters, and make it an ODE of the form -`dxₜ/dt = f(xₜ,λ(rt))`` +with `λ ∈ Rᵐ` containing all the system parameters, and make it an ODE of the form +`dxₜ/dt = f(xₜ,λ(rt))` with `λ(t) ∈ Rᵐ`. """ mutable struct RateProtocol From 345865f348a97a6b2480101611ec333aed13fa76 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 10:34:17 +0100 Subject: [PATCH 62/76] added plot of lambda to example of RateSystem --- docs/src/examples/RateSystem.md | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index 3bd89718..b1d0ad8a 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -53,7 +53,20 @@ end p_lambda = [λ_max] # parameter of the function lambda ``` -Now, we define the RateProtocol that describes the non-autonomous period: + +We plot the function ``λ(p_lambda,t)`` + +```@example RateSystem +lambda_plotvals = [λ(p_lambda, t)[1] for t in -10.0:0.1:10.0] +figλ = Figure(); axsλ = Axis(figλ[1,1],xlabel="t",ylabel=L"\lambda") +lines!(axsλ,-10.0:0.1:10.0,lambda_plotvals,linewidth=2,label=L"\lambda(p_{\lambda}, t)") +axislegend(axsλ,position=:rc,labelsize=10) +figλ +``` + + +Now, we define the RateProtocol that describes how and when the function ``λ(p_lambda,t)`` is used to +make autonomous system non-autonomous: ```@example RateSystem r = 4/3-0.02 # r just below critical rate @@ -82,7 +95,7 @@ nonauto_traj = trajectory(nonauto_sys,T,x0) We plot the two trajectories ```@example RateSystem -fig = Figure(); axs = Axis(fig[1,1]) +fig = Figure(); axs = Axis(fig[1,1],xlabel="t",ylabel="x") lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}") lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}") axislegend(axs,position=:rc,labelsize=10) From a6e80295e0fab9b0a0eb9812796fd602cbaafffa Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 10:45:53 +0100 Subject: [PATCH 63/76] removed unnecessary formatting information and improved beginning of RateSystem example --- docs/src/examples/RateSystem.md | 20 ++++++-------------- 1 file changed, 6 insertions(+), 14 deletions(-) diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md index b1d0ad8a..55f2cce7 100644 --- a/docs/src/examples/RateSystem.md +++ b/docs/src/examples/RateSystem.md @@ -1,17 +1,6 @@ -# Setting up a `RateSystem` -Let us explore an example of how to construct a `RateSystem` from an autonomous dynamical system (e.g. a `CoupledODEs`) and a time-dependent forcing protocol called `RateProtocol`. +# Setting up a `RateSystem`: a prototypical model for R-tipping -```@example RateSystem -using CriticalTransitions -using CairoMakie -using CairoMakie.Makie.MathTeXEngine: get_font -font = (; - regular=get_font(:regular), bold=get_font(:bold), - italic=get_font(:italic), bold_italic=get_font(:bolditalic) -); -``` - -## Prototypical model for R-tipping, with critical rate r = 4/3 +Let us explore an example of how to construct a `RateSystem` from an autonomous deterministic dynamical system (i.e. a `CoupledODEs`) and a time-dependent forcing protocol called `RateProtocol`. We first consider the following simple one-dimensional autonomous system with one attractor, given by the ordinary differential equation: ```math @@ -23,6 +12,9 @@ The parameter ``\lambda`` shifts the location of the extrema of the drift field. We implement this system as follows: ```@example RateSystem +using CriticalTransitions +using CairoMakie + function f(u,p,t) # out-of-place x = u[1] λ = p[1] @@ -69,7 +61,7 @@ Now, we define the RateProtocol that describes how and when the function ``λ(p_ make autonomous system non-autonomous: ```@example RateSystem -r = 4/3-0.02 # r just below critical rate +r = 4/3-0.02 # r just below critical rate 4/3 t_start = -Inf # start time of non-autonomous part t_end = Inf # end time of non-autonomous part From 3c9427915793700a6ad7cc65b3ffd82e9a473f87 Mon Sep 17 00:00:00 2001 From: raphael-roemer Date: Tue, 29 Jul 2025 14:57:09 +0100 Subject: [PATCH 64/76] correction in docstring --- src/r_tipping/RateSystem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/r_tipping/RateSystem.jl b/src/r_tipping/RateSystem.jl index f1f668c7..88dc7517 100644 --- a/src/r_tipping/RateSystem.jl +++ b/src/r_tipping/RateSystem.jl @@ -72,8 +72,8 @@ end; apply_ramping(sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs...) Applies a time-dependent [`RateProtocol`](@def) to a given autonomous deterministic dynamical system -`sys`, turning it into a non-autonomous dynamical system. Returns a [`CoupledODEs`](@ref) -with the explicit parameter time-dependence incorporated. +`sys`, turning it into a non-autonomous dynamical system. The returned [`CoupledODEs`](@ref) +has the explicit parameter time-dependence incorporated. The returned [`CoupledODEs`](@ref) is autonomous from `t_0` to `t_start`, then non-autnonmous from `t_start` to `t_end` with the parameter shift given by the [`RateProtocol`](@def), From 5f1f6bc4ab11525519e768fced831f7f0e4ea614 Mon Sep 17 00:00:00 2001 From: reykboerner Date: Tue, 29 Jul 2025 18:47:48 +0200 Subject: [PATCH 65/76] recovered docs/Project.toml --- docs/Project.toml | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index c15b0ddb..39a06969 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,26 @@ [deps] +Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7" +Contour = "d38c429a-6771-53c6-b99e-75d170b6e991" CriticalTransitions = "251e6cd3-3112-48a5-99dd-66efcfd18334" +DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df" +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" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71" +OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948" +Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" + +[compat] +Documenter = "^1.4.1" \ No newline at end of file From ef732865cc3ed2416572b1ce45f7dcb906c1ea0f Mon Sep 17 00:00:00 2001 From: reykboerner Date: Tue, 29 Jul 2025 18:54:21 +0200 Subject: [PATCH 66/76] updated .toml and applied Formatter --- Project.toml | 6 +++--- src/r_tipping/RateSystem.jl | 4 +--- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index 6960e4c4..7340b9e5 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "CriticalTransitions" uuid = "251e6cd3-3112-48a5-99dd-66efcfd18334" authors = ["Reyk Börner, Orjan Ameye, Ryan Deeley, Raphael Römer"] repo = "https://github.com/juliadynamics/CriticalTransitions.jl.git" -version = "0.6.0" +version = "0.6.1" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" @@ -45,7 +45,7 @@ DocStringExtensions = "0.9.3" DynamicalSystemsBase = "3.13" Format = "1" Interpolations = "0.15.1, 0.16" -IntervalArithmetic = "0.22.22" +IntervalArithmetic = "0.22.22, 0.23" LinearAlgebra = "1.10" LinearSolve = "2.37.0, 3" Optimization = "4" @@ -64,4 +64,4 @@ julia = "1.10" [extras] Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7" -ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7" +ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7" \ No newline at end of file diff --git a/src/r_tipping/RateSystem.jl b/src/r_tipping/RateSystem.jl index 88dc7517..06f4a029 100644 --- a/src/r_tipping/RateSystem.jl +++ b/src/r_tipping/RateSystem.jl @@ -81,9 +81,7 @@ and again autonomous from `t_end` to the end of the simulation: `t_0` autonomous `t_start` non-autonomous `t_end` autonomous `∞` """ -function apply_ramping( - auto_sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs... -) +function apply_ramping(auto_sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs...) # we wish to return a continuous time dynamical system with modified drift field f(u, p, t) = modified_drift( From 7b79902de909a4325b4724bc90245395e6b08bab Mon Sep 17 00:00:00 2001 From: reykboerner Date: Tue, 29 Jul 2025 21:40:48 +0200 Subject: [PATCH 67/76] added source files --- src/r_tipping/critical_rate.jl | 0 src/r_tipping/moving_sinks.jl | 0 2 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 src/r_tipping/critical_rate.jl create mode 100644 src/r_tipping/moving_sinks.jl diff --git a/src/r_tipping/critical_rate.jl b/src/r_tipping/critical_rate.jl new file mode 100644 index 00000000..e69de29b diff --git a/src/r_tipping/moving_sinks.jl b/src/r_tipping/moving_sinks.jl new file mode 100644 index 00000000..e69de29b From be7aca80a306031e7742039f5f39db6af620a0a9 Mon Sep 17 00:00:00 2001 From: reykboerner Date: Fri, 1 Aug 2025 00:45:00 +0200 Subject: [PATCH 68/76] added moving_sinks draft --- src/r_tipping/moving_sinks.jl | 36 ++++++++++++++++++++++++++++++++++ test/r_tipping/moving_sinks.jl | 30 ++++++++++++++++++++++++++++ 2 files changed, 66 insertions(+) create mode 100644 test/r_tipping/moving_sinks.jl diff --git a/src/r_tipping/moving_sinks.jl b/src/r_tipping/moving_sinks.jl index e69de29b..715a91ba 100644 --- a/src/r_tipping/moving_sinks.jl +++ b/src/r_tipping/moving_sinks.jl @@ -0,0 +1,36 @@ +""" + moving_sinks(ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; kwargs...) + +Calculates fixed points of a nonautonomous dynamical system at snapshots in time, giving +insight into how the equilibria change as the forcing changes over time. + +## Input +`ds` is a dynamical system of type `ContinuousTimeDynamicalSystem`, `rp` a `RateProtocol`, +and `box` is a vector of intervals (e.g. `[interval(0,1), interval(-1,1)]`) +that delimits the phase space volume in which to look for equilibria +(see [`fixedpoints`](@ref)). + +## Keyword arguments +- `times = 0:0.1:1`: time points (relative to the system time `t`) + +## output +Returns a triple of vectors: +1. Fixed points found at each time point +2. Eigenvalues associated with the fixed points +3. Stability of the fixed points (1 - stable; 0 - unstable) + +The term "moving sinks" refers to Wieczorek et al. (2023). +""" +function moving_sinks(ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; + times=0:0.1:1) + + fp, eig, stab = [], [], [] + for t in times + rate_sys = apply_ramping(ds, rp, t) + _fp, _eig, _stab = fixedpoints(rate_sys, box) + push!(fp, _fp) + push!(eig, _eig) + push!(stab, _stab) + end + return fp, eig, stab +end \ No newline at end of file diff --git a/test/r_tipping/moving_sinks.jl b/test/r_tipping/moving_sinks.jl new file mode 100644 index 00000000..036a4d83 --- /dev/null +++ b/test/r_tipping/moving_sinks.jl @@ -0,0 +1,30 @@ + +@testset "moving_sinks" begin + using ChaosTools + # Dynamical system + function fhn(u,p,t) + eps, b = p + x, y = u + dx = (-x^3 + x - y)/eps + dy = -b*y + x + return SA[dx, dy] + end + + p = [1,6] + sys = CoupledODEs(fhn, [1.,1], p) + + # Forcing + function λ(p, t) + λ_max = p[2] + lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) + return SVector{2}([p[1], lambda]) + end + + r = 0.1 + rp = CriticalTransitions.RateProtocol(λ, p, r, -10, 10) + + # Calculate moving sinks + box = [interval(-2,2), interval(-1,1)] + fp, eig, stab = moving_sinks(sys, rp, box; times=0:0.1:1) + @test length(fp[1]) == 3 +end \ No newline at end of file From 35b6c3070b026f0f4952af669c950ef7d82f1c5d Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Fri, 1 Aug 2025 08:21:09 +0200 Subject: [PATCH 69/76] apply formatter and enable spelling again --- .github/workflows/SpellCheck.yml | 1 + src/r_tipping/critical_rate.jl | 1 + src/r_tipping/moving_sinks.jl | 8 ++++---- test/r_tipping/moving_sinks.jl | 14 +++++++------- 4 files changed, 13 insertions(+), 11 deletions(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 17be9084..94597ca3 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -5,6 +5,7 @@ on: types: - opened - reopened + - synchronize - ready_for_review jobs: diff --git a/src/r_tipping/critical_rate.jl b/src/r_tipping/critical_rate.jl index e69de29b..8b137891 100644 --- a/src/r_tipping/critical_rate.jl +++ b/src/r_tipping/critical_rate.jl @@ -0,0 +1 @@ + diff --git a/src/r_tipping/moving_sinks.jl b/src/r_tipping/moving_sinks.jl index 715a91ba..98715aed 100644 --- a/src/r_tipping/moving_sinks.jl +++ b/src/r_tipping/moving_sinks.jl @@ -21,9 +21,9 @@ Returns a triple of vectors: The term "moving sinks" refers to Wieczorek et al. (2023). """ -function moving_sinks(ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; - times=0:0.1:1) - +function moving_sinks( + ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; times=0:0.1:1 +) fp, eig, stab = [], [], [] for t in times rate_sys = apply_ramping(ds, rp, t) @@ -33,4 +33,4 @@ function moving_sinks(ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; push!(stab, _stab) end return fp, eig, stab -end \ No newline at end of file +end diff --git a/test/r_tipping/moving_sinks.jl b/test/r_tipping/moving_sinks.jl index 036a4d83..efd0fba5 100644 --- a/test/r_tipping/moving_sinks.jl +++ b/test/r_tipping/moving_sinks.jl @@ -2,16 +2,16 @@ @testset "moving_sinks" begin using ChaosTools # Dynamical system - function fhn(u,p,t) - eps, b = p + function fhn(u, p, t) + eps, b = p x, y = u dx = (-x^3 + x - y)/eps dy = -b*y + x return SA[dx, dy] end - p = [1,6] - sys = CoupledODEs(fhn, [1.,1], p) + p = [1, 6] + sys = CoupledODEs(fhn, [1.0, 1], p) # Forcing function λ(p, t) @@ -19,12 +19,12 @@ lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) return SVector{2}([p[1], lambda]) end - + r = 0.1 rp = CriticalTransitions.RateProtocol(λ, p, r, -10, 10) # Calculate moving sinks - box = [interval(-2,2), interval(-1,1)] + box = [interval(-2, 2), interval(-1, 1)] fp, eig, stab = moving_sinks(sys, rp, box; times=0:0.1:1) @test length(fp[1]) == 3 -end \ No newline at end of file +end From 92740dc69a89859a3989a6be800ab3f90cbd0571 Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Fri, 1 Aug 2025 22:40:56 +0200 Subject: [PATCH 70/76] Added draft version of the critical_rate function. --- src/r_tipping/critical_rate.jl | 138 +++++++++++++++++++++++++++++++++ 1 file changed, 138 insertions(+) diff --git a/src/r_tipping/critical_rate.jl b/src/r_tipping/critical_rate.jl index 8b137891..7715b11a 100644 --- a/src/r_tipping/critical_rate.jl +++ b/src/r_tipping/critical_rate.jl @@ -1 +1,139 @@ +## the following function integrates the nonautonomous system created by the ContinuousTimeDynamicalSystem ds and RateProtocol rp, with rate-parameter r +## the initial condition is e_start at time t_start +## the system is integrated for T = t_end-t_start time units and the trajectory position at time T is then compared with the location of e_end +## if the trajectory position at time T is within a neighbourhood of radius rad of e_end, the function returns true, and false otherwise +function end_point_tracking( + ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, r::Float64, e_start::Vector, t_start::Float64, e_end::Vector, t_end::Float64, rad::Float64 +) + + # note that should ensure that Δt << dλ(rt)/dt + # more generally, the suitability of the numerical integrator diffeq has to be ensured + + T = t_end - t_start + rp.r = r # set the rate-parameter to the minimum value + rate_sys_min_rate = apply_ramping(ds, rp, t_start) # create the nonautonomous system + traj_min_rate = trajectory(rate_sys_min_rate, T, e_start) # simulate the nonautonomous system + + return abs(traj_min_rate[1][end,:]-e_end) < rad + +end + +# the following function finds via brute force an interval [a,b] satisfying b-a ≤ tol satisfying f(a)*f(b) == false for a function f returning boolean values +# i.e. a small interval [a,b] for which the value of the function f evaluated at the points a and b returns both true and false + +function binary_bisection(f::Function, a::Float64, b::Float64; + tol::Float64=1e-12, maxiter::Int64=10000) + a < b || error("Please enter a, b such that a < b.") + fa = f(a) + fa*f(b) == false || error("The function does not satisfy the initial requirement that f(a)*f(b) == false") + ii = 0 + local c + while b-a > tol + ii += 1 + ii ≤ maxiter || error("The maximum number of iterations was exceeded.") + c = (a+b)/2 + fc = f(c) + if fa*fc == true # reduce to the shrunken interval [c,b], since we must have f(c)*f(b) == false + a = c + fa = fc + else # reduce to the shrunken interval [a,c], since we have f(a)*f(c) == false + b = c + end + end + return [a,b] +end; + +## the following function aims to find a critical value for the rate-parameter between successful and non-successful "end-point tracking" + +## the end-point tracking is currently defined with respect to a stable equilibria of the frozen-in system attained at time t = t_start... +##...and a stable equilibria of the frozen-in system attained at time t = t_end + +## there is currently no check that e_start and e_end are connected via a continuous branch of moving sinks attained across [t_start,t_end],... +##... which would be a desirable check to perform although requires more work to implement (for a first draft I skipped this) +## before doing this it probably makes sense to modify the moving_sinks function so as to return continuous branches of grouped equilibria... +##...rather than just a vector of vectors containing all the (non-grouped) equilibria found for a range of frozen-in times + +## there are check for ensuring that e_start and e_end are indeed stable equilibria of the frozen-in systems attained at times t_start and t_end respectively +## there are also checks that the system succesfully end-point tracks for r = rmin and fails to end-point track for r = rmax (or vice versa, for exotic cases) + +## the function ultimately passes the end_point_tracking function to the binary_bisection function to find a critical rate for end-point tracking + + +function critical_end_point_tracking_rate( + ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, e_start::Vector, t_start::Float64, e_end::Vector, t_end::Float64; + rmin::Float64 = 1.e-2, + rmax::Float64 = 1.e2, + rad::Float64 = 1.e-1, + tol::Float64 = 1.e-3, + maxiter::Int64 = Int64(1.e4) +) + + rmin < rmax || error("Please enter rmin, rmax such that rmin < rmax.") + + ##-----check no.1-----## + + ## that e_start is a stable equilibrium of the frozen-in system attained at t = t_start + + rate_sys_start = apply_ramping(ds, rp, t_start) + box_start = [interval(x-0.1,x+0.1) for x ∈ e_start] + fp_start, ~, stab_start = fixedpoints(rate_sys_start, box_start) + if any(e -> isapprox(e,e_start;atol=1e-4),fp_start) + idx = findmin([abs(e-e_start) for e ∈ fp_start])[2] # the index in fp_start of the fixed point which is approximately equal to e_end + stab_start[idx] || @warn("The vector e_start = $(e_start) does not describe a stable equilibrium of the frozen-in system attained at t = t_start = $(t_start).") + else + @warn("The vector e_start = $(e_start) does not describe a stable equilibrium of the frozen-in system attained at t = t_start = $(t_start).") + end + + ##-----check no.2-----## + + ## that e_end is a stable equilibrium of the frozen-in system attained at t = t_end + + rate_sys_end = apply_ramping(ds, rp, t_end) + box_end = [interval(x-0.1,x+0.1) for x ∈ e_end] + fp_end, ~, stab_end = fixedpoints(rate_sys_end, box_end) + if any(e -> isapprox(e,e_end;atol=1e-4),fp_end) + idx = findmin([abs(e-e_end) for e ∈ fp_end])[2] # the index in fp_end of the fixed point which is approximately equal to e_end + stab_end[idx] || @warn("The vector e_end = $(e_end) does not describe a stable equilibrium of the frozen-in system attained at t = t_end = $(t_end).") + else + @warn("The vector e_end = $(e_end) does not describe a stable equilibrium of the frozen-in system attained at t = t_end = $(t_end).") + end + + ##-----check no.3-----## + + ## that + ## the system successfully end-point tracks for r = rmin and fails to end-point track for r = rmax + ## or + ## the system fails to end-point track for r = rmin and successfully end point tracks for r = rmax + + min_rate_track = end_point_tracking(ds,rp,rmin,e_start,t_start,e_end,t_end,rad) # true if the system end-point tracks for that rate + + max_rate_track = end_point_tracking(ds,rp,rmax,e_start,t_start,e_end,t_end,rad) # true if the system end-point tracks for that rate + + if min_rate_track && max_rate_track + error("The system successfully end-point tracks for both r = rmin = $(rmin) and r = rmax = $(rmax).") + elseif ~min_rate_track && ~max_rate_track + error("The system fails to end-point track for both r = rmin = $(rmin) and r = rmax = $(rmax).") + elseif ~min_rate_track && max_rate_track + @warn("The system fails to end-point track for r = rmin = $(rmin) but succesfully end-point tracks for r = rmax = $(rmax).") + end + + ## all checks complete ## + + ## performing the bisection to find a critical rate within the given tolerance ## + + func(r) = end_point_tracking(ds,rp,r,e_start,t_start,e_end,t_end,rad) + rcrit_interval = binary_bisection(func,rmin,rmax;tol,maxiter) + + return (rcrit_interval[1]+rcrit_interval[2])/2 + +end + +# beginning notes + +# you are trying to locate the critical rate r for which the system end-point tracks the moving branch of sinks + +# the first case is where after some travel time you are within a defined neighbourhood of a sink belonging to the frozen-in system at that time +# for this you should provide the tuple of tuples ((e₋,t₋),(e₊,t₊)) +# then you start at the position e⁻ at the time t⁻ and simulate t⁺-t⁻ time-units and compute the distance to e⁺ +# we assume that there exists some pairing (r₁,r₂) whereby the system lies within the neighbourhood of e⁺ under the rate r₁ and outside of the neighbourhood of e⁺ under the rate r₂ From f65f61d5c4720df3beac36362b3bce9c064a0561 Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Mon, 4 Aug 2025 00:00:47 +0200 Subject: [PATCH 71/76] Slight corrections to align with how fixed_points works. --- src/r_tipping/moving_sinks.jl | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/src/r_tipping/moving_sinks.jl b/src/r_tipping/moving_sinks.jl index 98715aed..6b3ca97e 100644 --- a/src/r_tipping/moving_sinks.jl +++ b/src/r_tipping/moving_sinks.jl @@ -21,16 +21,34 @@ Returns a triple of vectors: The term "moving sinks" refers to Wieczorek et al. (2023). """ +# function moving_sinks( +# ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; times=0:0.1:1 +# ) +# fp, eig, stab = [], [], [] +# for t in times +# rate_sys = apply_ramping(ds, rp, t) +# _fp, _eig, _stab = fixedpoints(rate_sys, box) +# push!(fp, _fp) +# push!(eig, _eig) +# push!(stab, _stab) +# end +# return fp, eig, stab +# end + +# modified version function moving_sinks( - ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; times=0:0.1:1 + ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; + times=rp.t_start/rp.r:(rp.t_end/rp.r-rp.t_start/rp.r)/100:rp.t_end/rp.r ) fp, eig, stab = [], [], [] + λ_mod = rp.λ for t in times - rate_sys = apply_ramping(ds, rp, t) + rp.λ = (p,τ) -> λ_mod(p,τ+t) + rate_sys = apply_ramping(ds, rp, rp.t_start) _fp, _eig, _stab = fixedpoints(rate_sys, box) push!(fp, _fp) push!(eig, _eig) push!(stab, _stab) end return fp, eig, stab -end +end \ No newline at end of file From 4243f90ac3e03fbe54df872b79bc0c3ecd7ffc18 Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Mon, 4 Aug 2025 00:13:27 +0200 Subject: [PATCH 72/76] Made corrections to the critical_rate function in the process of getting it to work on a test example. --- src/r_tipping/critical_rate.jl | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/src/r_tipping/critical_rate.jl b/src/r_tipping/critical_rate.jl index 7715b11a..30893283 100644 --- a/src/r_tipping/critical_rate.jl +++ b/src/r_tipping/critical_rate.jl @@ -1,3 +1,5 @@ +using LinearAlgebra: norm + ## the following function integrates the nonautonomous system created by the ContinuousTimeDynamicalSystem ds and RateProtocol rp, with rate-parameter r ## the initial condition is e_start at time t_start ## the system is integrated for T = t_end-t_start time units and the trajectory position at time T is then compared with the location of e_end @@ -12,10 +14,10 @@ function end_point_tracking( T = t_end - t_start rp.r = r # set the rate-parameter to the minimum value - rate_sys_min_rate = apply_ramping(ds, rp, t_start) # create the nonautonomous system - traj_min_rate = trajectory(rate_sys_min_rate, T, e_start) # simulate the nonautonomous system + nonauto_sys = apply_ramping(ds, rp, t_start) # create the nonautonomous system + traj = trajectory(nonauto_sys, T, e_start) # simulate the nonautonomous system - return abs(traj_min_rate[1][end,:]-e_end) < rad + return norm(traj[1][end,:]-e_end) < rad end @@ -60,7 +62,7 @@ end; ## the function ultimately passes the end_point_tracking function to the binary_bisection function to find a critical rate for end-point tracking -function critical_end_point_tracking_rate( +function critical_rate( ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, e_start::Vector, t_start::Float64, e_end::Vector, t_end::Float64; rmin::Float64 = 1.e-2, rmax::Float64 = 1.e2, @@ -75,11 +77,14 @@ function critical_end_point_tracking_rate( ## that e_start is a stable equilibrium of the frozen-in system attained at t = t_start - rate_sys_start = apply_ramping(ds, rp, t_start) + λ_mod = rp.λ + + rp.λ = (p,τ) -> λ_mod(p,τ+t_start) + rate_sys_start = apply_ramping(ds, rp, rp.t_start) box_start = [interval(x-0.1,x+0.1) for x ∈ e_start] fp_start, ~, stab_start = fixedpoints(rate_sys_start, box_start) if any(e -> isapprox(e,e_start;atol=1e-4),fp_start) - idx = findmin([abs(e-e_start) for e ∈ fp_start])[2] # the index in fp_start of the fixed point which is approximately equal to e_end + idx = findmin([norm(e-e_start) for e ∈ fp_start])[2] # the index in fp_start of the fixed point which is approximately equal to e_end stab_start[idx] || @warn("The vector e_start = $(e_start) does not describe a stable equilibrium of the frozen-in system attained at t = t_start = $(t_start).") else @warn("The vector e_start = $(e_start) does not describe a stable equilibrium of the frozen-in system attained at t = t_start = $(t_start).") @@ -89,11 +94,12 @@ function critical_end_point_tracking_rate( ## that e_end is a stable equilibrium of the frozen-in system attained at t = t_end - rate_sys_end = apply_ramping(ds, rp, t_end) + rp.λ = (p,τ) -> λ_mod(p,τ+t_end) + rate_sys_end = apply_ramping(ds, rp, rp.t_start) box_end = [interval(x-0.1,x+0.1) for x ∈ e_end] fp_end, ~, stab_end = fixedpoints(rate_sys_end, box_end) if any(e -> isapprox(e,e_end;atol=1e-4),fp_end) - idx = findmin([abs(e-e_end) for e ∈ fp_end])[2] # the index in fp_end of the fixed point which is approximately equal to e_end + idx = findmin([norm(e-e_end) for e ∈ fp_end])[2] # the index in fp_end of the fixed point which is approximately equal to e_end stab_end[idx] || @warn("The vector e_end = $(e_end) does not describe a stable equilibrium of the frozen-in system attained at t = t_end = $(t_end).") else @warn("The vector e_end = $(e_end) does not describe a stable equilibrium of the frozen-in system attained at t = t_end = $(t_end).") @@ -106,15 +112,17 @@ function critical_end_point_tracking_rate( ## or ## the system fails to end-point track for r = rmin and successfully end point tracks for r = rmax + rp.λ = (p,τ) -> λ_mod(p,τ) + min_rate_track = end_point_tracking(ds,rp,rmin,e_start,t_start,e_end,t_end,rad) # true if the system end-point tracks for that rate max_rate_track = end_point_tracking(ds,rp,rmax,e_start,t_start,e_end,t_end,rad) # true if the system end-point tracks for that rate if min_rate_track && max_rate_track error("The system successfully end-point tracks for both r = rmin = $(rmin) and r = rmax = $(rmax).") - elseif ~min_rate_track && ~max_rate_track + elseif !min_rate_track && !max_rate_track error("The system fails to end-point track for both r = rmin = $(rmin) and r = rmax = $(rmax).") - elseif ~min_rate_track && max_rate_track + elseif !min_rate_track && max_rate_track @warn("The system fails to end-point track for r = rmin = $(rmin) but succesfully end-point tracks for r = rmax = $(rmax).") end @@ -136,4 +144,4 @@ end # the first case is where after some travel time you are within a defined neighbourhood of a sink belonging to the frozen-in system at that time # for this you should provide the tuple of tuples ((e₋,t₋),(e₊,t₊)) # then you start at the position e⁻ at the time t⁻ and simulate t⁺-t⁻ time-units and compute the distance to e⁺ -# we assume that there exists some pairing (r₁,r₂) whereby the system lies within the neighbourhood of e⁺ under the rate r₁ and outside of the neighbourhood of e⁺ under the rate r₂ +# we assume that there exists some pairing (r₁,r₂) whereby the system lies within the neighbourhood of e⁺ under the rate r₁ and outside of the neighbourhood of e⁺ under the rate r₂ \ No newline at end of file From eea5bb50b116816c19521c494a216a96c22d3b7d Mon Sep 17 00:00:00 2001 From: Ryan Deeley Date: Mon, 4 Aug 2025 01:15:32 +0200 Subject: [PATCH 73/76] Slight corrections to align with how fixed_points works. --- src/r_tipping/moving_sinks.jl | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/src/r_tipping/moving_sinks.jl b/src/r_tipping/moving_sinks.jl index 6b3ca97e..b46bc4d6 100644 --- a/src/r_tipping/moving_sinks.jl +++ b/src/r_tipping/moving_sinks.jl @@ -42,8 +42,24 @@ function moving_sinks( ) fp, eig, stab = [], [], [] λ_mod = rp.λ + t_start = rp.t_start + t_end = rp.t_end for t in times - rp.λ = (p,τ) -> λ_mod(p,τ+t) + + t̃ = if r*t ≤ t_start + t_start + elseif t_start < r*t < t_end + r*t + else + t_end + end; # the value of the argument of the λ function in the drift field at time t + + # the fixed point function computes fixed points for the autonomous system attained at time t = 0 + # ensuring that rp.t_start < 0 < rp.t_end and shifting the rp.λ function so that modified drift calls λ_mod(p,t̃) = λ_mod(p,r*0+t̃), rather than λ_mod(p,t_start+t̃) or λ_mod(p,t_end+t̃) for instance + rp.t_start = -1.0 + rp.t_end = 1.0 + rp.λ = (p,τ) -> λ_mod(p,τ+t̃) + rate_sys = apply_ramping(ds, rp, rp.t_start) _fp, _eig, _stab = fixedpoints(rate_sys, box) push!(fp, _fp) From 7b77fc0f977ea622989475c9e3f72557feb7a8a4 Mon Sep 17 00:00:00 2001 From: reykboerner Date: Wed, 6 Aug 2025 21:23:45 +0200 Subject: [PATCH 74/76] fixed docs error --- docs/pages.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/pages.jl b/docs/pages.jl index 771cd9fb..6e203e4d 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -5,7 +5,7 @@ pages = [ "Tutorial" => "examples/tutorial.md", "Examples" => Any[ "Defining stochastic systems" => "examples/stochastic-dynamics.md", - "Setting up a RateSystem" => "examples/RateSystem.md" + "Setting up a RateSystem" => "examples/RateSystem.md", "Large deviations: Maier-Stein system" => "examples/gMAM_Maierstein.md", "Simple gMAM: Kerr Parametric Oscillator" => "examples/sgMAM_KPO.md", "Transition Path Theory: Finite element method" => "examples/transition_path_theory_double_well.md", From a55936f27c539ceb18fea5773e552cccc2354334 Mon Sep 17 00:00:00 2001 From: reykboerner Date: Tue, 4 Nov 2025 17:08:32 +0100 Subject: [PATCH 75/76] delete deprecated files --- docs/src/examples/RateSystem.md | 100 -------------------------------- src/r_tipping/RateSystem.jl | 93 ----------------------------- test/r_tipping/RateSystem.jl | 41 ------------- 3 files changed, 234 deletions(-) delete mode 100644 docs/src/examples/RateSystem.md delete mode 100644 src/r_tipping/RateSystem.jl delete mode 100644 test/r_tipping/RateSystem.jl diff --git a/docs/src/examples/RateSystem.md b/docs/src/examples/RateSystem.md deleted file mode 100644 index 55f2cce7..00000000 --- a/docs/src/examples/RateSystem.md +++ /dev/null @@ -1,100 +0,0 @@ -# Setting up a `RateSystem`: a prototypical model for R-tipping - -Let us explore an example of how to construct a `RateSystem` from an autonomous deterministic dynamical system (i.e. a `CoupledODEs`) and a time-dependent forcing protocol called `RateProtocol`. - -We first consider the following simple one-dimensional autonomous system with one attractor, given by the ordinary differential equation: -```math -\begin{aligned} - \dot{x} &= (x+\lambda)^2 - 1 -\end{aligned} -``` -The parameter ``\lambda`` shifts the location of the extrema of the drift field. -We implement this system as follows: - -```@example RateSystem -using CriticalTransitions -using CairoMakie - -function f(u,p,t) # out-of-place - x = u[1] - λ = p[1] - dx = (x+λ)^2 - 1 - return SVector{1}(dx) -end - -lambda = 0.0 -p = [lambda] -x0 = [-1.] -auto_sys = CoupledODEs(f,x0,p) -``` - -## Non-autonomous case - -Now, we want to explore a non-autonomous version of the system. -We consider a setting where in the past and in the future the system is autnonomous and in between there is a non-autonomous period ``[t_start, t_end]`` with a time-dependent parameter ramping given by the function ``\lambda(rt)``. Choosing different values of the parameter ``r`` allows us to vary the speed of the parameter ramping. - -We start by defining the function ``\lambda(t)``: -```@example RateSystem -function λ(p,t) - λ_max = p[1] - lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) - return SVector{1}(lambda) -end - -λ_max = 3. -p_lambda = [λ_max] # parameter of the function lambda -``` - - -We plot the function ``λ(p_lambda,t)`` - -```@example RateSystem -lambda_plotvals = [λ(p_lambda, t)[1] for t in -10.0:0.1:10.0] -figλ = Figure(); axsλ = Axis(figλ[1,1],xlabel="t",ylabel=L"\lambda") -lines!(axsλ,-10.0:0.1:10.0,lambda_plotvals,linewidth=2,label=L"\lambda(p_{\lambda}, t)") -axislegend(axsλ,position=:rc,labelsize=10) -figλ -``` - - -Now, we define the RateProtocol that describes how and when the function ``λ(p_lambda,t)`` is used to -make autonomous system non-autonomous: - -```@example RateSystem -r = 4/3-0.02 # r just below critical rate 4/3 -t_start = -Inf # start time of non-autonomous part -t_end = Inf # end time of non-autonomous part - -rp = CriticalTransitions.RateProtocol(λ,p_lambda,r,t_start,t_end) -``` - - - -Now, we set up the combined system with autonomous past and future and non-autonomous ramping in between: - -```@example RateSystem -t0 = -10. # initial time of the system -nonauto_sys = apply_ramping(auto_sys,rp,t0) -``` - -We can compute trajectories of this new system in the familiar way: -```@example RateSystem -T = 20. # final simulation time -auto_traj = trajectory(auto_sys,T,x0) -nonauto_traj = trajectory(nonauto_sys,T,x0) -``` - -We plot the two trajectories - -```@example RateSystem -fig = Figure(); axs = Axis(fig[1,1],xlabel="t",ylabel="x") -lines!(axs,t0.+auto_traj[2],auto_traj[1][:,1],linewidth=2,label=L"\text{Autonomous system}") -lines!(axs,nonauto_traj[2],nonauto_traj[1][:,1],linewidth=2,label=L"\text{Nonautonomous system}") -axislegend(axs,position=:rc,labelsize=10) -fig -``` - ------ -Author: Raphael Roemer - -Date: 30 Jun 2025 diff --git a/src/r_tipping/RateSystem.jl b/src/r_tipping/RateSystem.jl deleted file mode 100644 index 06f4a029..00000000 --- a/src/r_tipping/RateSystem.jl +++ /dev/null @@ -1,93 +0,0 @@ -# we consider the ODE dxₜ/dt = f(xₜ,λ(rt)) -# here λ = λ(t) ∈ Rᵐ is a function containing all the system parameters - -# We ask the user to define: -# 1) a ContinuousTimeDynamicalSystem that should be investigated and -# 2) a protocol for the time-dependent forcing with the struct RateProtocol - -# Then we give back the ContinuousTimeDynamicalSystem with the parameter -# changing according to the rate protocol -""" - RateProtocol - -Time-dependent forcing protocol specified by the following fields: -- `λ::Function`: forcing function of the form `λ(p, t_start; kwargs...)`` -- `p_lambda::Vector`: parameters of forcing function -- `r::Float64`: rate parameter -- `t_start::Float64`: start time of protocol -- `t_end::Float64`: end time of protocol - -If `t_start` and `t_end` are not provided, they are set to `t_start=-Inf`, and `t_end=Inf`. -If `p_lambda` is not provided, it is set to an empty array `[]`. - -The forcing protocol contains all the information that allows to take an ODE of the form -`dxₜ/dt = f(xₜ,λ)` -with `λ ∈ Rᵐ` containing all the system parameters, and make it an ODE of the form -`dxₜ/dt = f(xₜ,λ(rt))` -with `λ(t) ∈ Rᵐ`. -""" -mutable struct RateProtocol - λ::Function - p_lambda::Vector - r::Float64 - t_start::Float64 - t_end::Float64 -end - -# convenience functions - -function RateProtocol(λ::Function, p_lambda::Vector, r::Float64) - RateProtocol(λ, p_lambda, r, -Inf, Inf) -end -RateProtocol(λ::Function, r::Float64)=RateProtocol(λ, [], r, -Inf, Inf) -#RateProtocol(λ::Function,p_lambda::Vector,r::Float64,t_start::Float64)=RateProtocol(λ,p_lambda,r,t_start,Inf) -#RateProtocol(λ::Function,r::Float64,t_start::Float64)=RateProtocol(λ,[],r,t_start,Inf) - -function modified_drift( - u, - p, - t, - ds::ContinuousTimeDynamicalSystem, - λ::Function, - t_start::Float64, - t_end::Float64, - r::Float64; - kwargs..., -) - if t_start > t_end - error("Please ensure that t_start ≤ t_end.") - end - - p̃ = if r*t ≤ t_start - λ(p, t_start; kwargs...) - elseif t_start < r*t < t_end - λ(p, r*t; kwargs...) # the value(s) of λ(rt) - else - λ(p, t_end; kwargs...) # the value(s) of λ(rt) - end; # the value(s) of λ(rt) - return ds.integ.f(u, p̃, t) -end; - -""" - apply_ramping(sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs...) - -Applies a time-dependent [`RateProtocol`](@def) to a given autonomous deterministic dynamical system -`sys`, turning it into a non-autonomous dynamical system. The returned [`CoupledODEs`](@ref) -has the explicit parameter time-dependence incorporated. - -The returned [`CoupledODEs`](@ref) is autonomous from `t_0` to `t_start`, -then non-autnonmous from `t_start` to `t_end` with the parameter shift given by the [`RateProtocol`](@def), -and again autonomous from `t_end` to the end of the simulation: - -`t_0` autonomous `t_start` non-autonomous `t_end` autonomous `∞` -""" -function apply_ramping(auto_sys::CoupledODEs, rp::RateProtocol, t0=0.0; kwargs...) - # we wish to return a continuous time dynamical system with modified drift field - - f(u, p, t) = modified_drift( - u, p, t, auto_sys, rp.λ, rp.t_start, rp.t_end, rp.r; kwargs... - ) - prob = remake(auto_sys.integ.sol.prob; f, p=rp.p_lambda, tspan=(t0, Inf)) - nonauto_sys = CoupledODEs(prob, auto_sys.diffeq) - return nonauto_sys -end diff --git a/test/r_tipping/RateSystem.jl b/test/r_tipping/RateSystem.jl deleted file mode 100644 index be24c286..00000000 --- a/test/r_tipping/RateSystem.jl +++ /dev/null @@ -1,41 +0,0 @@ -using CriticalTransitions -using Test - -@testset "RateSystem" begin - function f(u, p, t) # out-of-place - x = u[1] - λ = p[1] - dx = (x+λ)^2 - 1 - return SVector{1}(dx) - end - - lambda = 0.0 - p = [lambda] - x0 = [-1.0] - auto_sys = CoupledODEs(f, x0, p) - - function λ(p, t) - λ_max = p[1] - lambda = (λ_max/2)*(tanh(λ_max*t/2)+1) - return SVector{1}(lambda) - end - - λ_max = 3.0 - p_lambda = [λ_max] # parameter of the function lambda - - r = 4/3-0.02 # r just below critical rate - t_start = -Inf # start time of non-autonomous part - t_end = Inf # end time of non-autonomous part - - rp = CriticalTransitions.RateProtocol(λ, p_lambda, r, t_start, t_end) - - t0 = -10.0 # initial time of the system - nonauto_sys = apply_ramping(auto_sys, rp, t0) - - T = 20.0 # final simulation time - auto_traj = trajectory(auto_sys, T, x0) - nonauto_traj = trajectory(nonauto_sys, T, x0) - - @test isapprox(auto_traj[1][end, 1], -1; atol=1e-4) - @test isapprox(nonauto_traj[1][end, 1], -4; atol=1e-4) -end From 304e2cace32419a00493b707f924d2a51f791fcf Mon Sep 17 00:00:00 2001 From: reykboerner Date: Tue, 4 Nov 2025 17:36:36 +0100 Subject: [PATCH 76/76] update moving_sinks, fix spelling and bump version --- Project.toml | 2 +- src/r_tipping/critical_rate.jl | 123 ++++++++++++++++++++------------- src/r_tipping/moving_sinks.jl | 89 ++++++++++++------------ 3 files changed, 120 insertions(+), 94 deletions(-) diff --git a/Project.toml b/Project.toml index df06ad55..17998899 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "CriticalTransitions" uuid = "251e6cd3-3112-48a5-99dd-66efcfd18334" authors = ["Reyk Börner, Orjan Ameye, Ryan Deeley, Raphael Römer", "George Datseris"] repo = "https://github.com/juliadynamics/CriticalTransitions.jl.git" -version = "0.7.0" +version = "0.7.1" [deps] ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" diff --git a/src/r_tipping/critical_rate.jl b/src/r_tipping/critical_rate.jl index 30893283..ad09a7d5 100644 --- a/src/r_tipping/critical_rate.jl +++ b/src/r_tipping/critical_rate.jl @@ -3,10 +3,17 @@ using LinearAlgebra: norm ## the following function integrates the nonautonomous system created by the ContinuousTimeDynamicalSystem ds and RateProtocol rp, with rate-parameter r ## the initial condition is e_start at time t_start ## the system is integrated for T = t_end-t_start time units and the trajectory position at time T is then compared with the location of e_end -## if the trajectory position at time T is within a neighbourhood of radius rad of e_end, the function returns true, and false otherwise +## if the trajectory position at time T is within a neighborhood of radius rad of e_end, the function returns true, and false otherwise function end_point_tracking( - ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, r::Float64, e_start::Vector, t_start::Float64, e_end::Vector, t_end::Float64, rad::Float64 + ds::ContinuousTimeDynamicalSystem, + rp::RateProtocol, + r::Float64, + e_start::Vector, + t_start::Float64, + e_end::Vector, + t_end::Float64, + rad::Float64, ) # note that should ensure that Δt << dλ(rt)/dt @@ -17,18 +24,20 @@ function end_point_tracking( nonauto_sys = apply_ramping(ds, rp, t_start) # create the nonautonomous system traj = trajectory(nonauto_sys, T, e_start) # simulate the nonautonomous system - return norm(traj[1][end,:]-e_end) < rad - + return norm(traj[1][end, :]-e_end) < rad end # the following function finds via brute force an interval [a,b] satisfying b-a ≤ tol satisfying f(a)*f(b) == false for a function f returning boolean values # i.e. a small interval [a,b] for which the value of the function f evaluated at the points a and b returns both true and false -function binary_bisection(f::Function, a::Float64, b::Float64; - tol::Float64=1e-12, maxiter::Int64=10000) - a < b || error("Please enter a, b such that a < b.") +function binary_bisection( + f::Function, a::Float64, b::Float64; tol::Float64=1e-12, maxiter::Int64=10000 +) + a < b || error("Please enter a, b such that a < b.") fa = f(a) - fa*f(b) == false || error("The function does not satisfy the initial requirement that f(a)*f(b) == false") + fa*f(b) == false || error( + "The function does not satisfy the initial requirement that f(a)*f(b) == false" + ) ii = 0 local c while b-a > tol @@ -37,13 +46,13 @@ function binary_bisection(f::Function, a::Float64, b::Float64; c = (a+b)/2 fc = f(c) if fa*fc == true # reduce to the shrunken interval [c,b], since we must have f(c)*f(b) == false - a = c - fa = fc + a = c + fa = fc else # reduce to the shrunken interval [a,c], since we have f(a)*f(c) == false - b = c + b = c end end - return [a,b] + return [a, b] end; ## the following function aims to find a critical value for the rate-parameter between successful and non-successful "end-point tracking" @@ -57,91 +66,107 @@ end; ##...rather than just a vector of vectors containing all the (non-grouped) equilibria found for a range of frozen-in times ## there are check for ensuring that e_start and e_end are indeed stable equilibria of the frozen-in systems attained at times t_start and t_end respectively -## there are also checks that the system succesfully end-point tracks for r = rmin and fails to end-point track for r = rmax (or vice versa, for exotic cases) +## there are also checks that the system successfully end-point tracks for r = rmin and fails to end-point track for r = rmax (or vice versa, for exotic cases) ## the function ultimately passes the end_point_tracking function to the binary_bisection function to find a critical rate for end-point tracking - function critical_rate( - ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, e_start::Vector, t_start::Float64, e_end::Vector, t_end::Float64; - rmin::Float64 = 1.e-2, - rmax::Float64 = 1.e2, - rad::Float64 = 1.e-1, - tol::Float64 = 1.e-3, - maxiter::Int64 = Int64(1.e4) + ds::ContinuousTimeDynamicalSystem, + rp::RateProtocol, + e_start::Vector, + t_start::Float64, + e_end::Vector, + t_end::Float64; + rmin::Float64=1.e-2, + rmax::Float64=1.e2, + rad::Float64=1.e-1, + tol::Float64=1.e-3, + maxiter::Int64=Int64(1.e4), ) - - rmin < rmax || error("Please enter rmin, rmax such that rmin < rmax.") + rmin < rmax || error("Please enter rmin, rmax such that rmin < rmax.") ##-----check no.1-----## ## that e_start is a stable equilibrium of the frozen-in system attained at t = t_start - λ_mod = rp.λ + λ_mod = rp.λ - rp.λ = (p,τ) -> λ_mod(p,τ+t_start) + rp.λ = (p, τ) -> λ_mod(p, τ+t_start) rate_sys_start = apply_ramping(ds, rp, rp.t_start) - box_start = [interval(x-0.1,x+0.1) for x ∈ e_start] + box_start = [interval(x-0.1, x+0.1) for x in e_start] fp_start, ~, stab_start = fixedpoints(rate_sys_start, box_start) - if any(e -> isapprox(e,e_start;atol=1e-4),fp_start) - idx = findmin([norm(e-e_start) for e ∈ fp_start])[2] # the index in fp_start of the fixed point which is approximately equal to e_end - stab_start[idx] || @warn("The vector e_start = $(e_start) does not describe a stable equilibrium of the frozen-in system attained at t = t_start = $(t_start).") + if any(e -> isapprox(e, e_start; atol=1e-4), fp_start) + idx = findmin([norm(e-e_start) for e in fp_start])[2] # the index in fp_start of the fixed point which is approximately equal to e_end + stab_start[idx] || @warn( + "The vector e_start = $(e_start) does not describe a stable equilibrium of the frozen-in system attained at t = t_start = $(t_start)." + ) else - @warn("The vector e_start = $(e_start) does not describe a stable equilibrium of the frozen-in system attained at t = t_start = $(t_start).") + @warn( + "The vector e_start = $(e_start) does not describe a stable equilibrium of the frozen-in system attained at t = t_start = $(t_start)." + ) end ##-----check no.2-----## - + ## that e_end is a stable equilibrium of the frozen-in system attained at t = t_end - rp.λ = (p,τ) -> λ_mod(p,τ+t_end) + rp.λ = (p, τ) -> λ_mod(p, τ+t_end) rate_sys_end = apply_ramping(ds, rp, rp.t_start) - box_end = [interval(x-0.1,x+0.1) for x ∈ e_end] + box_end = [interval(x-0.1, x+0.1) for x in e_end] fp_end, ~, stab_end = fixedpoints(rate_sys_end, box_end) - if any(e -> isapprox(e,e_end;atol=1e-4),fp_end) - idx = findmin([norm(e-e_end) for e ∈ fp_end])[2] # the index in fp_end of the fixed point which is approximately equal to e_end - stab_end[idx] || @warn("The vector e_end = $(e_end) does not describe a stable equilibrium of the frozen-in system attained at t = t_end = $(t_end).") + if any(e -> isapprox(e, e_end; atol=1e-4), fp_end) + idx = findmin([norm(e-e_end) for e in fp_end])[2] # the index in fp_end of the fixed point which is approximately equal to e_end + stab_end[idx] || @warn( + "The vector e_end = $(e_end) does not describe a stable equilibrium of the frozen-in system attained at t = t_end = $(t_end)." + ) else - @warn("The vector e_end = $(e_end) does not describe a stable equilibrium of the frozen-in system attained at t = t_end = $(t_end).") + @warn( + "The vector e_end = $(e_end) does not describe a stable equilibrium of the frozen-in system attained at t = t_end = $(t_end)." + ) end ##-----check no.3-----## - + ## that ## the system successfully end-point tracks for r = rmin and fails to end-point track for r = rmax ## or ## the system fails to end-point track for r = rmin and successfully end point tracks for r = rmax - rp.λ = (p,τ) -> λ_mod(p,τ) + rp.λ = (p, τ) -> λ_mod(p, τ) - min_rate_track = end_point_tracking(ds,rp,rmin,e_start,t_start,e_end,t_end,rad) # true if the system end-point tracks for that rate + min_rate_track = end_point_tracking(ds, rp, rmin, e_start, t_start, e_end, t_end, rad) # true if the system end-point tracks for that rate - max_rate_track = end_point_tracking(ds,rp,rmax,e_start,t_start,e_end,t_end,rad) # true if the system end-point tracks for that rate + max_rate_track = end_point_tracking(ds, rp, rmax, e_start, t_start, e_end, t_end, rad) # true if the system end-point tracks for that rate - if min_rate_track && max_rate_track - error("The system successfully end-point tracks for both r = rmin = $(rmin) and r = rmax = $(rmax).") + if min_rate_track && max_rate_track + error( + "The system successfully end-point tracks for both r = rmin = $(rmin) and r = rmax = $(rmax).", + ) elseif !min_rate_track && !max_rate_track - error("The system fails to end-point track for both r = rmin = $(rmin) and r = rmax = $(rmax).") + error( + "The system fails to end-point track for both r = rmin = $(rmin) and r = rmax = $(rmax).", + ) elseif !min_rate_track && max_rate_track - @warn("The system fails to end-point track for r = rmin = $(rmin) but succesfully end-point tracks for r = rmax = $(rmax).") + @warn( + "The system fails to end-point track for r = rmin = $(rmin) but succesfully end-point tracks for r = rmax = $(rmax)." + ) end ## all checks complete ## ## performing the bisection to find a critical rate within the given tolerance ## - func(r) = end_point_tracking(ds,rp,r,e_start,t_start,e_end,t_end,rad) - rcrit_interval = binary_bisection(func,rmin,rmax;tol,maxiter) + func(r) = end_point_tracking(ds, rp, r, e_start, t_start, e_end, t_end, rad) + rcrit_interval = binary_bisection(func, rmin, rmax; tol, maxiter) return (rcrit_interval[1]+rcrit_interval[2])/2 - end # beginning notes # you are trying to locate the critical rate r for which the system end-point tracks the moving branch of sinks -# the first case is where after some travel time you are within a defined neighbourhood of a sink belonging to the frozen-in system at that time +# the first case is where after some travel time you are within a defined neighborhood of a sink belonging to the frozen-in system at that time # for this you should provide the tuple of tuples ((e₋,t₋),(e₊,t₊)) # then you start at the position e⁻ at the time t⁻ and simulate t⁺-t⁻ time-units and compute the distance to e⁺ -# we assume that there exists some pairing (r₁,r₂) whereby the system lies within the neighbourhood of e⁺ under the rate r₁ and outside of the neighbourhood of e⁺ under the rate r₂ \ No newline at end of file +# we assume that there exists some pairing (r₁,r₂) whereby the system lies within the neighborhood of e⁺ under the rate r₁ and outside of the neighborhood of e⁺ under the rate r₂ diff --git a/src/r_tipping/moving_sinks.jl b/src/r_tipping/moving_sinks.jl index b46bc4d6..06bd66e5 100644 --- a/src/r_tipping/moving_sinks.jl +++ b/src/r_tipping/moving_sinks.jl @@ -1,17 +1,17 @@ """ - moving_sinks(ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; kwargs...) + moving_fixedpoints(rs::RateSystem, box; kwargs...) Calculates fixed points of a nonautonomous dynamical system at snapshots in time, giving insight into how the equilibria change as the forcing changes over time. ## Input -`ds` is a dynamical system of type `ContinuousTimeDynamicalSystem`, `rp` a `RateProtocol`, +`rs` is a dynamical system of type [`RateSystem`](@ref) and `box` is a vector of intervals (e.g. `[interval(0,1), interval(-1,1)]`) that delimits the phase space volume in which to look for equilibria -(see [`fixedpoints`](@ref)). +(see [`Attractors.fixedpoints`](@ref)). ## Keyword arguments -- `times = 0:0.1:1`: time points (relative to the system time `t`) +- `Δt = 0.1`: Time step between fixed point calculations (in system time units) ## output Returns a triple of vectors: @@ -21,50 +21,51 @@ Returns a triple of vectors: The term "moving sinks" refers to Wieczorek et al. (2023). """ -# function moving_sinks( -# ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; times=0:0.1:1 -# ) -# fp, eig, stab = [], [], [] -# for t in times -# rate_sys = apply_ramping(ds, rp, t) -# _fp, _eig, _stab = fixedpoints(rate_sys, box) -# push!(fp, _fp) -# push!(eig, _eig) -# push!(stab, _stab) -# end -# return fp, eig, stab -# end - -# modified version -function moving_sinks( - ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; - times=rp.t_start/rp.r:(rp.t_end/rp.r-rp.t_start/rp.r)/100:rp.t_end/rp.r -) +function moving_fixedpoints(rs::RateSystem, box; Δt=0.1) fp, eig, stab = [], [], [] - λ_mod = rp.λ - t_start = rp.t_start - t_end = rp.t_end + t_start = rs.forcing.forcing_start_time + t_end = rs.forcing_forcing_start_time + rs.forcing.forcing_duration + times = t_start:Δt:t_end + for t in times - - t̃ = if r*t ≤ t_start - t_start - elseif t_start < r*t < t_end - r*t - else - t_end - end; # the value of the argument of the λ function in the drift field at time t - - # the fixed point function computes fixed points for the autonomous system attained at time t = 0 - # ensuring that rp.t_start < 0 < rp.t_end and shifting the rp.λ function so that modified drift calls λ_mod(p,t̃) = λ_mod(p,r*0+t̃), rather than λ_mod(p,t_start+t̃) or λ_mod(p,t_end+t̃) for instance - rp.t_start = -1.0 - rp.t_end = 1.0 - rp.λ = (p,τ) -> λ_mod(p,τ+t̃) - - rate_sys = apply_ramping(ds, rp, rp.t_start) - _fp, _eig, _stab = fixedpoints(rate_sys, box) + _fp, _eig, _stab = fixedpoints(frozen_system(rs, t), box) push!(fp, _fp) push!(eig, _eig) push!(stab, _stab) end return fp, eig, stab -end \ No newline at end of file +end + +# modified version +#function moving_sinks( +# ds::ContinuousTimeDynamicalSystem, rp::RateProtocol, box; +# times=rp.t_start/rp.r:(rp.t_end/rp.r-rp.t_start/rp.r)/100:rp.t_end/rp.r +#) +# fp, eig, stab = [], [], [] +# λ_mod = rp.λ +# t_start = rp.t_start +# t_end = rp.t_end +# for t in times +# +# t̃ = if r*t ≤ t_start +# t_start +# elseif t_start < r*t < t_end +# r*t +# else +# t_end +# end; # the value of the argument of the λ function in the drift field at time t +# +# # the fixed point function computes fixed points for the autonomous system attained at time t = 0 +# # ensuring that rp.t_start < 0 < rp.t_end and shifting the rp.λ function so that modified drift calls λ_mod(p,t̃) = λ_mod(p,r*0+t̃), rather than λ_mod(p,t_start+t̃) or λ_mod(p,t_end+t̃) for instance +# rp.t_start = -1.0 +# rp.t_end = 1.0 +# rp.λ = (p,τ) -> λ_mod(p,τ+t̃) +# +# rate_sys = apply_ramping(ds, rp, rp.t_start) +# _fp, _eig, _stab = fixedpoints(rate_sys, box) +# push!(fp, _fp) +# push!(eig, _eig) +# push!(stab, _stab) +# end +# return fp, eig, stab +#end