R-tipping functionality

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"

docs/pages.jl

@@ -5,6 +5,7 @@ pages = [
     "Tutorial" => "examples/tutorial.md",
     "Examples" => Any[
         "Defining stochastic systems" => "examples/stochastic-dynamics.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",

src/r_tipping/critical_rate.jl

@@ -0,0 +1,172 @@
+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 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,
+)
+
+    # 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 
+    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
+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 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),
+)
+    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.λ
+
+    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 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 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)."
+        )
+    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)
+    rate_sys_end = apply_ramping(ds, rp, rp.t_start)
+    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 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)."
+        )
+    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, τ)
+
+    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 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 neighborhood of e⁺ under the rate r₁ and outside of the neighborhood of e⁺ under the rate r₂ 

src/r_tipping/moving_sinks.jl

@@ -0,0 +1,71 @@
+"""
+    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
+`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 [`Attractors.fixedpoints`](@ref)).
+
+## Keyword arguments
+- `Δt = 0.1`: Time step between fixed point calculations (in system time units)
+
+## 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_fixedpoints(rs::RateSystem, box; Δt=0.1)
+    fp, eig, stab = [], [], []
+    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
+        _fp, _eig, _stab = fixedpoints(frozen_system(rs, t), 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
+#)
+#    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

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