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

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;