R-tipping functionality

.github/workflows/SpellCheck.yml

@@ -5,7 +5,6 @@ on:
     types:
       - opened
       - reopened
-      - synchronize
       - ready_for_review
 
 jobs:

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"
@@ -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"

docs/Project.toml

@@ -23,4 +23,4 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 
 [compat]
-Documenter = "^1.4.1"
+Documenter = "^1.4.1"

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",

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:
 

docs/src/examples/RateSystem.md

@@ -0,0 +1,100 @@
+# 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

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`.

src/CriticalTransitions.jl

@@ -62,6 +62,8 @@ 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
 
@@ -86,6 +88,8 @@ export find_boundary, huniform, dunion
 export committor,
     invariant_pdf, reactive_current, probability_reactive, probability_last_A, Langevin
 
+export RateProtocol, apply_ramping
+
 # Error hint for extensions stubs
 function __init__()
     Base.Experimental.register_error_hint(

src/r_tipping/RateSystem.jl

@@ -0,0 +1,93 @@
+# 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

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

systems/CTLibrary.jl

@@ -13,8 +13,8 @@ include("stommel.jl")
 include("rivals.jl")
 
 export fitzhugh_nagumo!, fitzhugh_nagumo, stommel, rivals!, rivals, cessi, rooth_smooth
-modifiedtruscottbrindleywithdimensions!, modifiedtruscottbrindleywithdimensions
-originaltruscottbrindley!, originaltruscottbrindley
-rampedoriginaltruscottbrindley!, rampedoriginaltruscottbrindley
+#modifiedtruscottbrindleywithdimensions!, modifiedtruscottbrindleywithdimensions
+#originaltruscottbrindley!, originaltruscottbrindley
+#rampedoriginaltruscottbrindley!, rampedoriginaltruscottbrindley
 
 end

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.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

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