Using new CoupledSDEs

Project.toml

@@ -45,7 +45,7 @@ Dierckx = "0.5.3"
 DiffEqNoiseProcess = "5.22"
 DocStringExtensions = "0.9.3"
 Documenter = "^1.4.1"
-DynamicalSystemsBase = "3.7.1"
+DynamicalSystemsBase = "3.11"
 ExplicitImports = "1.9"
 Format = "1"
 ForwardDiff = "0.10.36"

src/CriticalTransitions.jl

@@ -2,16 +2,10 @@ module CriticalTransitions
 
 # Base
 using Statistics: Statistics, mean
-using LinearAlgebra: LinearAlgebra, Diagonal, I, norm, tr, dot
+using LinearAlgebra: LinearAlgebra, I, norm, dot, tr #, Diagonal, tr
 using StaticArrays: StaticArrays, SVector
 
-# core
-using DynamicalSystemsBase:
-    DynamicalSystemsBase,
-    ContinuousTimeDynamicalSystem,
-    StateSpaceSets,
-    dimension,
-    dynamic_rule
+# Core
 using DiffEqNoiseProcess: DiffEqNoiseProcess
 using OrdinaryDiffEq: OrdinaryDiffEq, Tsit5
 using StochasticDiffEq:
@@ -25,6 +19,19 @@ using StochasticDiffEq:
     step!,
     terminate!,
     u_modified!
+using DynamicalSystemsBase:
+    DynamicalSystemsBase,
+    #ContinuousTimeDynamicalSystem,
+    CoupledODEs,
+    CoupledSDEs,
+    #StateSpaceSets,
+    #dimension,
+    dynamic_rule,
+    current_state,
+    set_state!
+
+#StochasticSystemsBase = Base.get_extension(DynamicalSystemsBase, :StochasticSystemsBase)
+#using DynamicalSystemsBase.StochasticSystemsBase: CoupledSDEs
 
 using ForwardDiff: ForwardDiff
 using IntervalArithmetic: IntervalArithmetic, interval
@@ -49,10 +56,15 @@ using Reexport: @reexport
 @reexport using StochasticDiffEq
 @reexport using DiffEqNoiseProcess
 
+StochasticSystemsBase = Base.get_extension(DynamicalSystemsBase, :StochasticSystemsBase)
+covariance_matrix = StochasticSystemsBase.covariance_matrix
+diffusion_matrix = StochasticSystemsBase.diffusion_matrix
+
 include("extention_functions.jl")
 include("utils.jl")
-include("CoupledSDEs.jl")
-include("CoupledSDEs_utils.jl")
+include("system_utils.jl")
+#include("CoupledSDEs.jl")
+#include("CoupledSDEs_utils.jl")
 include("io.jl")
 include("trajectories/simulation.jl")
 include("trajectories/transition.jl")
@@ -70,9 +82,9 @@ export CoupledSDEs,
     idfunc!,
     idfunc,
     add_noise_strength,
-    noise_process,
-    covariance_matrix,
-    noise_strength,
+    #noise_process,
+    #covariance_matrix,
+    #noise_strength,
     CoupledODEs
 
 # Methods
@@ -82,6 +94,7 @@ export transition, transitions
 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 basins, basinboundary
 # export edgetracking, bisect_to_edge, AttractorsViaProximity
 # export fixedpoints

src/largedeviations/action.jl

@@ -20,7 +20,7 @@ generalized norm ``||a||_Q^2 := \\langle a, Q^{-1} b \\rangle`` (see `anorm``).
 function fw_action(sys::CoupledSDEs, path, time)
     @assert all(diff(time) .≈ diff(time[1:2])) "Freidlin-Wentzell action is only defined for equispaced time"
     # Inverse of covariance matrix
-    A = inv(covariance_matrix(sys)) # Do we have to take the inverse here?
+    A = inv(covariance_matrix(sys)) 
 
     # Compute action integral
     integrand = fw_integrand(sys, path, time, A)
@@ -35,7 +35,7 @@ end;
 """
 $(TYPEDSIGNATURES)
 
-Calculates the Onsager-Machlup action of a given `path` with time points `time` for the drift field `sys.f` at given `sys.noise_strength`.
+Calculates the Onsager-Machlup action of a given `path` with time points `time` for the drift field `sys.f` at given `noise_strength`.
 
 The path must be a `(D x N)` matrix, where `D` is the dimensionality of the system `sys` and
 `N` is the number of path points. The `time` array must have length `N`.
@@ -50,10 +50,10 @@ time of the path, and ``\\sigma`` the noise strength. The subscript ``Q`` refers
 generalized norm ``||a||_Q^2 := \\langle a, Q^{-1} b \\rangle`` (see `anorm``). Here
 ``Q`` is the noise covariance matrix.
 """
-function om_action(sys::CoupledSDEs, path, time)
+function om_action(sys::CoupledSDEs, path, time, noise_strength)
     @assert all(diff(time) .≈ diff(time[1:2])) "Fw_action is only defined for equispaced time"
 
-    σ = noise_strength(sys)
+    σ = noise_strength
     # Compute action integral
     S = 0
     for i in 1:(size(path, 2) - 1)
@@ -75,11 +75,11 @@ Computes the action functional specified by `functional` for a given CoupledSDEs
 * `functional = "FW"`: Returns the Freidlin-Wentzell action ([`fw_action`](@ref))
 * `functional = "OM"`: Returns the Onsager-Machlup action ([`om_action`](@ref))
 """
-function action(sys::CoupledSDEs, path::Matrix, time, functional)
+function action(sys::CoupledSDEs, path::Matrix, time, functional; noise_strength=nothing)
     if functional == "FW"
         action = fw_action(sys, path, time)
     elseif functional == "OM"
-        action = om_action(sys, path, time)
+        action = om_action(sys, path, time, noise_strength)
     end
     return action
 end;

src/system_utils.jl

@@ -0,0 +1,27 @@
+"""
+$(TYPEDSIGNATURES)
+
+Returns the deterministic drift ``f(x)`` of the CoupledSDEs `sys` at state `x`. For time-
+dependent systems, the time can be specified as a keyword argument `t` (by default `t=0`).
+"""
+function drift(sys::CoupledSDEs{IIP}, x; t=0) where {IIP}
+    f = dynamic_rule(sys)
+    if IIP
+        dx = similar(x)
+        f(dx, x, sys.p0, t)
+        return dx
+    else
+        return f(x, sys.p0, t)
+    end
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Computes the divergence of the drift field ``f(x)`` at state `x`. For time-
+dependent systems, the time can be specified as a keyword argument `t` (by default `t=0`).
+"""
+function div_drift(sys::CoupledSDEs, x; t=0)
+    b(x) = drift(sys, x; t)
+    return tr(ForwardDiff.jacobian(b, x))
+end;

src/trajectories/equib.jl

@@ -1,5 +1,5 @@
 """
-    equilib(sys::CoupledSDEs, state; kwargs...)
+    equilib(sys::CoupledSDEs [, state]; kwargs...)
 Returns the equilibrium solution of the system `sys` for given initial condition `state`.
 
 > Warning: This algorithm simply evolves the deterministic system forward in time until a steady-state condition is satisfied.
@@ -12,13 +12,12 @@ Returns the equilibrium solution of the system `sys` for given initial condition
 * `dt = 0.01`: time step of the ODE solver.
 * `solver = Euler()`: ODE solver used for evolving the state.
 """
-function equilib(sys::CoupledSDEs, state; dt=0.01, tmax=1e5, abstol=1e-5, solver=Tsit5())
+function equilib(sys::CoupledSDEs, state=nothing; dt=0.01, tmax=1e5, abstol=1e-5, solver=Tsit5())
     condition(u, t, integrator) = norm(integrator.uprev - u) < abstol
     affect!(integrator) = terminate!(integrator)
     equilib_cond = DiscreteCallback(condition, affect!)
-
-    prob = ODEProblem(sys.integ.f, state, (0, tmax), sys.p0)
+    isnothing(state) ? nothing : set_state!(sys, state)
+    prob = sys.integ.sol.prob
     sol = solve(prob, solver; dt=dt, callback=equilib_cond, save_on=false, save_start=false)
-
     return sol.u[1]
 end;

src/trajectories/simulation.jl

@@ -15,7 +15,7 @@ For more info, see [`SDEProblem`](https://diffeq.sciml.ai/stable/types/sde_types
 function simulate(
     sys::CoupledSDEs, T, init=current_state(sys); alg=sys.integ.alg, kwargs...
 )
-    prob = remake(referrenced_sciml_prob(sys); u0=init, tspan=(0, T))
+    prob = remake(sys.integ.sol.prob; u0=init, tspan=(0, T))
     return solve(prob, alg; kwargs...)
 end;
 
@@ -35,7 +35,7 @@ For stochastic integration, see [`simulate`](@ref).
 
 """
 function relax(sys::CoupledSDEs, T, init=current_state(sys); alg=Tsit5(), kwargs...)
-    sde_prob = referrenced_sciml_prob(sys)
+    sde_prob = sys.integ.sol.prob
     prob = ODEProblem{isinplace(sde_prob)}(dynamic_rule(sys), init, (0, T), sys.p0)
     return solve(prob, alg; kwargs...)
 end;

test/CoupledSDEs.jl

@@ -21,7 +21,7 @@
         # diagonal additive noise: σ*N(0,dt)
         # a vector of random numbers dW whose size matches the output of g where the noise is applied element-wise
         prob = SDEProblem(meier_stein, diagonal_noise(σ), zeros(2), (0.0, Inf), ())
-        sde = CoupledSDEs(meier_stein, idfunc, zeros(2), (), σ)
+        sde = CoupledSDEs(meier_stein, zeros(2); noise_strength=σ)
 
         @test sde.integ.sol.prob.f isa SDEFunction
         @test sde.integ.sol.prob.f.f.f == prob.f.f
@@ -46,7 +46,7 @@
         σ = 0.25 # noise strength
         W = WienerProcess(0.0, 0.0, 0.0)
         prob = SDEProblem(f!, diagonal_noise!(σ), zeros(2), (0.0, Inf), (); noise=W)
-        sde = CoupledSDEs(f!, idfunc!, zeros(2), (), σ; noise=W)
+        sde = CoupledSDEs(f!, zeros(2); noise_process=W, noise_strength=σ)
 
         @test sde.integ.sol.prob.f isa SDEFunction
         @test sde.integ.sol.prob.f.f.f == prob.f.f
@@ -74,7 +74,7 @@
         g_sde(u, p, t) = σ .* u
         g_Csde(u, p, t) = σ .* u
         prob = SDEProblem(meier_stein, g_sde, zeros(2), (0.0, Inf), ())
-        sde = CoupledSDEs(meier_stein, g_Csde, zeros(2), (), σ)
+        sde = CoupledSDEs(meier_stein, zeros(2); g=g_Csde)
 
         @test sde.integ.sol.prob.f isa SDEFunction
         @test sde.integ.sol.prob.f.f.f == prob.f.f
@@ -108,7 +108,7 @@
         # prob = SDEProblem(f, g, ones(2), (0.0, 1.0), noise_rate_prototype = zeros(2, 4))
         # diffeq =(alg=SRIW1(), noise_rate_prototype = zeros(2, 4))
         sde = CoupledSDEs(
-            f, g, zeros(2); noise_rate_prototype=zeros(2, 4), diffeq=(alg=RKMilCommute(),)
+            f, zeros(2); g, noise_prototype=zeros(2, 4), diffeq=(alg=RKMilCommute(),)
         )
         prob = SDEProblem(f, g, zeros(2), (0.0, Inf); noise_rate_prototype=zeros(2, 4))
 
@@ -133,7 +133,7 @@
         Z0 = zeros(2)
         W = CorrelatedWienerProcess(Γ, t0, W0, Z0)
         prob = SDEProblem(f!, diagonal_noise!(σ), zeros(2), (0.0, Inf), (); noise=W)
-        sde = CoupledSDEs(f!, idfunc!, zeros(2), (), σ; noise=W)
+        sde = CoupledSDEs(f!, zeros(2); noise_process=W, noise_strength=σ)
 
         @test sde.integ.sol.prob.f isa SDEFunction
         @test sde.integ.sol.prob.f.f.f == prob.f.f
@@ -158,4 +158,4 @@
             @test W.covariance ≈ cov(sol; dims=2) / dt rtol = 1e-2
         end
     end
-end
+end

test/basin/basin_boundary-N54XDK9GMK.jl

@@ -0,0 +1,18 @@
+@testset "fitzhugh_nagumo" begin
+    # Define systems
+    g = diag_noise_funtion(0.215)
+    system = CoupledSDEs(fitzhugh_nagumo, g, [2.0, 0.1], [0.1, 3.0, 1.0, 1.0, 1.0, 0.0])
+
+    # Set up domain
+    A, B, C = [0.0, 0.0], [1.0, 0.0], [0.0, 1.0]
+    box = intervals_to_box([-2.0, -2.0], [2.00, 2.0])
+    step = 0.1
+
+    boas = basins(system, A, B, C, box; bstep = [step, step])
+    X, Y, attractors, h = boas
+    boundary = basinboundary(boas)
+    @test length(attractors) == 2
+    @test length(unique(h)) == 3
+    @test h[end÷2+1, end÷2+1] == -1
+    @test boundary[:,end÷2+1] == zeros(2)
+end

test/issue14.jl

@@ -0,0 +1,22 @@
+using CriticalTransitions, StaticArrays, Test
+
+function meier_stein(u, p, t) # in-place
+    x, y = u
+    dx = x-x^3 -10*x*y^2
+    dy = -(1+x^2)*y
+    SA[dx, dy]
+end
+σ = 0.25
+sys = StochSystem(meier_stein, [], zeros(2), σ, idfunc, nothing, I(2), "WhiteGauss")
+
+# initial path: parabola
+xx = range(-1.0,1.0, length=30)
+yy = 0.3 .* (- xx .^ 2 .+ 1)
+init = Matrix([xx yy]')
+
+gm = geometric_min_action_method(sys, init, maxiter=1000)
+@test all(isapprox.(path[2,:][end-5:end], 0, atol=1e-3))
+
+gm[2][end]
+
+isapprox(gm[2][end], 0.338, atol=1e-3)

test/largedeviations/MAM.jl

@@ -1,7 +1,7 @@
 @testset "MAM FitzHugh-Nagumo" begin
     p = [0.1, 3, 1, 1, 1, 0]
     σ = 0.1
-    fhn = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ)
+    fhn = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=σ)
     x_i = SA[sqrt(2 / 3), sqrt(2 / 27)]
     x_f = SA[0.001, 0.0]
     N, T = 200, 2.0
@@ -14,7 +14,7 @@ end
 
 @testset "MAM Ornstein-Uhlenbeck" begin
     σ = 0.18
-    ou = CoupledSDEs((u, p, t) -> -u, idfunc, SA[1.0], (), σ)
+    ou = CoupledSDEs((u, p, t) -> -u, SA[1.0]; noise_strength=σ)
     x0 = -1
     xT = 2.0
     T = 10.0

test/largedeviations/action_fhn.jl

@@ -4,9 +4,9 @@ Random.seed!(SEED)
 # System setup - FitzHugh-Nagumo model
 p = [1.0, 3.0, 1.0, 1.0, 1.0, 0.0] # Parameters (ϵ, β, α, γ, κ, I)
 σ = 0.2 # noise strength
-sys = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ; seed=SEED)
+sys = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=σ)#, seed=SEED)
 
-A = inv(CT.covariance_matrix(sys))
+A = inv(covariance_matrix(sys))
 T, N = 2.0, 100
 
 x_i = SA[sqrt(2 / 3), sqrt(2 / 27)]
@@ -22,7 +22,7 @@ end
 
 # Test om_action function
 @testset "om_action" begin
-    S = om_action(sys, path, time)
+    S = om_action(sys, path, time, σ)
     @test isapprox(S, 0.26, atol=0.01)
 end
 

test/largedeviations/gMAM.jl

@@ -1,7 +1,7 @@
 @testset "gMAM FitzHugh-Nagumo" begin
     p = [0.1, 3, 1, 1, 1, 0]
     σ = 0.1
-    fhn = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ)
+    fhn = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=σ)
     x_i = SA[sqrt(2 / 3), sqrt(2 / 27)]
     x_f = SA[0.001, 0.0]
     N = 100
@@ -19,7 +19,7 @@ end
     end
     σ = 0.25
     # sys = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ)
-    sys = CoupledSDEs(meier_stein, idfunc, zeros(2), (), σ)
+    sys = CoupledSDEs(meier_stein, zeros(2); noise_strength=σ)
 
     # initial path: parabola
     xx = range(-1.0, 1.0; length=30)

test/runtests.jl

@@ -44,9 +44,9 @@ end
     JET.test_package(CriticalTransitions; target_defined_modules=true)
 end
 
-@testset "CoupledSDEs" begin
-    include("CoupledSDEs.jl")
-end
+#@testset "CoupledSDEs" begin
+#    include("CoupledSDEs.jl")
+#end
 
 @testset "ModelingToolkit" begin
     include("ModelingToolkit.jl")
@@ -58,9 +58,9 @@ end
     include("largedeviations/gMAM.jl")
 end
 
-@testset "utilities" begin
-    include("utils.jl")
-end
+#@testset "utilities" begin
+#    include("utils.jl")
+#end
 
 @testset "Trajactories" begin
     include("trajactories/simulate.jl")

test/trajactories/simulate.jl

@@ -3,9 +3,9 @@
     p = [0.1, 3.0, 1.0, 1.0, 1.0, 0.0] # Parameters (ϵ, β, α, γ, κ, I)
     σ = 0.1
     # Simulate
-    sys = CoupledSDEs(fitzhugh_nagumo, idfunc, SA[1.0, 0.0], p, σ)
+    sys = CoupledSDEs(fitzhugh_nagumo, SA[1.0, 0.0], p; noise_strength=σ)
     traj = relax(sys, 10)
-    sys = CoupledSDEs(fitzhugh_nagumo, idfunc, SA[1.0, 0.0], p, σ)
+    sys = CoupledSDEs(fitzhugh_nagumo, SA[1.0, 0.0], p; noise_strength=σ)
     sim = simulate(sys, 10)
     @test traj[1][1, 1] == 1.0
     @test sim.u[1][1] == 1.0