Stochsystem overhaul

README.md

@@ -2,7 +2,6 @@
 
 [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://juliadynamics.github.io/CriticalTransitions.jl/dev/)
 [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliadynamics.github.io/CriticalTransitions.jl/stable/)
-
 [![Tests](https://github.com/JuliaDynamics/CriticalTransitions.jl/actions/workflows/ci.yml/badge.svg)](github.com/JuliaDynamics/CriticalTransitions.jl/actions/workflows/ci.yml)
 
 A Julia package for the numerical investigation of **noise- and rate-induced transitions in dynamical systems**.

src/CriticalTransitions.jl

@@ -24,7 +24,7 @@ include("trajectories/transition.jl")
 include("noiseprocesses/stochprocess.jl")
 include("largedeviations/action.jl")
 include("largedeviations/min_action_method.jl")
-# include("largedeviations/geometric_min_action_method.jl")
+include("largedeviations/geometric_min_action_method.jl")
 
 include("../systems/fitzhughnagumo.jl")
 include("../systems/truscottbrindley_mod.jl")
@@ -50,6 +50,7 @@ export equilib, fixedpoints, basins, basinboundary, basboundary
 export simulate, relax
 export transition, transitions
 export fw_integrand, fw_action, om_action, action, geometric_action
+export div_drift
 export min_action_method, geometric_min_action_method
 export edgetracking, bisect_to_edge, attractor_mapper, bisect_to_edge2
 export make_jld2, make_h5, sys_string, sys_info, intervals_to_box

src/largedeviations/action.jl

@@ -37,9 +37,9 @@ function fw_action(sys::CoupledSDEs, path, time; cov_inv=nothing)
 end;
 
 """
- om_action(sys::CoupledSDEs, path, time; kwargs...)
+ om_action(sys::CoupledSDEs, path, time, noise_strength; kwargs...)
 Calculates the Onsager-Machlup action of a given `path` with time points `time` for the
-drift field `sys.f` at noise strength `sys.σ`.
+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`.
@@ -52,13 +52,13 @@ Returns a single number, which is the value of the action functional
 where ``\\phi_t`` denotes the path in state space, ``b`` the drift field, ``T`` the total
 time of the path, and ``\\sigma`` the noise strength. The subscript ``Q`` refers to the
 generalized norm ``||a||_Q^2 := \\langle a, Q^{-1} b \\rangle`` (see [`anorm`](@ref)). Here
-``Q`` is the noise covariance matrix `sys.Σ`.
+``Q`` is the noise covariance matrix.
 
 ## Keyword arguments
 * `cov_inv = nothing`: Inverse of the covariance matrix ``\\Sigma``.
  If `nothing`, ``\\Sigma^{-1}`` is computed automatically.
 """
-function om_action(sys::CoupledSDEs, path, time; cov_inv=nothing)
+function om_action(sys::CoupledSDEs, path, time, noise_strength; cov_inv=nothing)
  if ~all(diff(time) .≈ diff(time[1:2]))
  error("Fw_action is only defined for equispaced time")
  end
@@ -67,8 +67,8 @@ function om_action(sys::CoupledSDEs, path, time; cov_inv=nothing)
  # Compute action integral
  S = 0
  for i in 1:(size(path, 2) - 1)
- S += noise_strength(sys)^2/2 * ((div_b(sys, path[:,i+1]) + div_b(sys, path[:,i]))/2 *
- (time[i+1]-time[i]))
+ S += noise_strength^2/2 * ((div_drift(sys, path[:,i+1])
+ + div_drift(sys, path[:,i]))/2 * (time[i+1]-time[i]))
  end
  fw_action(sys, path, time, cov_inv=A) + S/2
  end;
@@ -85,7 +85,8 @@ function action(sys::CoupledSDEs, path::Matrix, time, functional; kwargs...)
  if functional == "FW"
  return fw_action(sys, path, time; kwargs...)
  elseif functional == "OM"
- return om_action(sys, path, time; kwargs...)
+ error("Not yet implemented since it depends on noise strength.")
+ #return om_action(sys, path, time; kwargs...)
  end
 end;
 
@@ -151,10 +152,10 @@ function fw_integrand(sys::CoupledSDEs, path, time, A)
 end;
 
 """
- div_b(sys::CoupledSDEs, x)
+ div_drift(sys::CoupledSDEs, x)
 Computes the divergence of the drift field `sys.f` at the given point `x`.
 """
-function div_b(sys::CoupledSDEs, x)
+function div_drift(sys::CoupledSDEs, x)
  b(x) = drift(sys, x)
  tr(ForwardDiff.jacobian(b, x))
 end;

src/largedeviations/geometric_min_action_method.jl

@@ -2,14 +2,14 @@
 #include("action.jl")
 
 """
- geometric_min_action_method(sys::StochSystem, x_i::State, x_f::State, arclength=1; kwargs...)
+ geometric_min_action_method(sys::CoupledSDEs, x_i::SVector, x_f::SVector, arclength=1; kwargs...)
 
 Computes the minimizer of the Freidlin-Wentzell action using the geometric minimum
 action method (gMAM). Beta version, to be further documented.
 
 To set an initial path different from a straight line, see the multiple dispatch method
 
- - `geometric_min_action_method(sys::StochSystem, init::Matrix, arclength::Float64; kwargs...)`.
+ - `geometric_min_action_method(sys::CoupledSDEs, init::Matrix, arclength::Float64; kwargs...)`.
 
 ## Keyword arguments
 
@@ -29,7 +29,7 @@ algorithm[^1].
 [^1]: [Heymann and Vanden-Eijnden, PRL (2008)](https://link.aps.org/doi/10.1103/PhysRevLett.100.140601)
 """
 function geometric_min_action_method(
- sys::StochSystem, x_i::State, x_f::State, arclength = 1.0;
+ sys::CoupledSDEs, x_i::SVector, x_f::SVector, arclength = 1.0;
  N = 100,
  maxiter = 100,
  converge = 1e-5,
@@ -38,25 +38,25 @@ function geometric_min_action_method(
  verbose=true,
  showprogress=false)
  path = reduce(hcat, range(x_i, x_f, length = N))
- geometric_min_action_method(sys::StochSystem, path, arclength;
+ geometric_min_action_method(sys::CoupledSDEs, path, arclength;
  maxiter = maxiter, converge = converge,
  method = method, tau = tau, verbose = verbose, showprogress = showprogress)
 end
 
 """
- geometric_min_action_method(sys::StochSystem, init::Matrix, arclength::Float64; kwargs...)
+ geometric_min_action_method(sys::CoupledSDEs, init::Matrix, arclength::Float64; kwargs...)
 
 Runs the geometric Minimum Action Method (gMAM) to find the minimum action path (instanton) from an
 initial condition `init`, given a system `sys` and total arc length `arclength`.
 
-The initial path `init` must be a matrix of size `(D, N)`, where `D` is the dimension
-`length(sys.u)` of the system and `N` is the number of path points.
+The initial path `init` must be a matrix of size `(D, N)`, where `D` is the dimension of the
+system and `N` is the number of path points.
 
 For more information see the main method,
-[`geometric_min_action_method(sys::StochSystem, x_i::State, x_f::State, arclength::Float64; kwargs...)`](@ref).
+[`geometric_min_action_method(sys::CoupledSDEs, x_i::SVector, x_f::SVector, arclength::Float64; kwargs...)`](@ref).
 """
 
-function geometric_min_action_method(sys::StochSystem, init::Matrix, arclength=1.0;
+function geometric_min_action_method(sys::CoupledSDEs, init::Matrix, arclength=1.0;
  maxiter = 100,
  converge = 1e-5,
  method = LBFGS(),
@@ -66,7 +66,7 @@ function geometric_min_action_method(sys::StochSystem, init::Matrix, arclength=1
 
  verbose && println("=== Initializing gMAM action minimizer ===")
 
- A = inv(sys.Σ)
+ A = inv(covariance_matrix(sys))
  path = init
  x_i = init[:, 1]
  x_f = init[:, end]
@@ -106,15 +106,16 @@ end
 function interpolate_path(path, sys, N, arclength)
  s = zeros(N)
  for j in 2:N
- s[j] = s[j - 1] + anorm(path[:, j] - path[:, j - 1], sys.Σ) #! anorm or norm?
+ s[j] = s[j - 1] + anorm(path[:, j] - path[:, j - 1], covariance_matrix(sys))
+ #! anorm or norm?
  end
  s_length = s / s[end] * arclength
  interp = ParametricSpline(s_length, path, k = 3)
  return reduce(hcat, [interp(x) for x in range(0, arclength, length = N)])
 end
 
 """
- heymann_vandeneijnden_step(sys::StochSystem, path, N, L; kwargs...)
+ heymann_vandeneijnden_step(sys::CoupledSDEs, path, N, L; kwargs...)
 
 Solves eq. (6) of Ref.[^1] for an initial `path` with `N` points and arclength `L`.
 
@@ -126,12 +127,12 @@ Solves eq. (6) of Ref.[^1] for an initial `path` with `N` points and arclength `
 
 [^1]: [Heymann and Vanden-Eijnden, PRL (2008)](https://link.aps.org/doi/10.1103/PhysRevLett.100.140601)
 """
-function heymann_vandeneijnden_step(sys::StochSystem, path, N, L;
+function heymann_vandeneijnden_step(sys::CoupledSDEs, path, N, L;
  tau = 0.1,
  diff_order = 4,
  cov_inv = nothing)
 
- (cov_inv == nothing) ? A = inv(sys.Σ) : A = cov_inv
+ (cov_inv == nothing) ? A = inv(covariance_matrix(sys)) : A = cov_inv
 
  dx = L / (N - 1)
  update = zeros(size(path))

test/largedeviations/MAM.jl

@@ -1,12 +1,23 @@
-@testset "OuMAM" begin
- ou = StochSystem((u, p, t) -> -u, [], [1.0], 1.0)
+@testset "MAM FitzHugh-Nagumo" begin
+ p = [0.1,3,1,1,1,0]
+ fhn = CoupledSDEs(fitzhugh_nagumo, diag_noise_funtion(0.1), zeros(2), p)
+ x_i = SA[sqrt(2/3), sqrt(2/27)]
+ x_f = SA[0.001, 0.0]
+ N, T = 200, 2.0
+ inst = min_action_method(fhn, x_i, x_f, N, T; maxiter=100, save_info=false)
+ S = fw_action(fhn, inst, range(0., T, length=N))
+ @test isapprox(S, 0.18, atol=0.01)
+end
+
+@testset "MAM Ornstein-Uhlenbeck" begin
+ ou = CoupledSDEs((u, p, t) -> -u, diag_noise_funtion(0.18), SA[1.0])
  x0 = -1
  xT = 2.0
  T = 10.0
  N = 51
  t = range(0, T, N)
- inst_mam = min_action_method(ou, [x0], [xT], N, T, verbose = false, save_info = false)
+ inst_mam = min_action_method(ou, SA[x0], SA[xT], N, T, verbose = false, save_info = false)
  inst_sol = ((xT - x0 * exp(-T)) * exp.(t) .+ (x0 * exp(T) - xT) * exp.(-t)) /
  (exp(T) - exp(-T))
  @test maximum(abs.(inst_mam' .- inst_sol)) < 0.1
-end
+end

test/largedeviations/action_fhn.jl

@@ -0,0 +1,48 @@
+using Random
+Random.seed!(SEED)
+
+# System setup - FitzHugh-Nagumo model
+p = [1.0, 3.0, 1.0, 1.0, 1.0, 0.0] # Parameters (ϵ, β, α, γ, κ, I)
+σ = 0.215 # noise strength
+sys = CoupledSDEs(fitzhugh_nagumo, diag_noise_funtion(σ), zeros(2), p, seed = SEED)
+A = inv(CT.covariance_matrix(sys))
+T, N = 2.0, 100
+x_i = SA[sqrt(2/3), sqrt(2/27)]
+x_f = SA[0.0, 0.0]
+path = reduce(hcat, range(x_i, x_f, length=N))
+time = range(0.0, T, length=N)
+
+@testset "fw_action" begin
+ S = fw_action(sys, path, time)
+ @test isapprox(S, 0.32, atol=0.01)
+end
+
+# Test om_action function
+@testset "om_action" begin
+ sigma = 0.2
+ S = om_action(sys, path, time, sigma)
+ @test isapprox(S, 0.26, atol=0.01)
+end
+
+# Test action function
+@testset "action" begin
+ @test action(sys, path, time, "FW") == fw_action(sys, path, time)
+end
+
+# Test geometric_action function
+@testset "geometric_action" begin
+ S = geometric_action(sys, path)
+ @test isapprox(S, 0.23, atol=0.01)
+end
+
+# Test fw_integrand function
+@testset "fw_integrand" begin
+ integrand = fw_integrand(sys, path, time, A)
+ @test minimum(integrand) > 0.18
+end
+
+# Test div_drift function
+@testset "div_drift" begin
+ @test div_drift(sys, zeros(2)) == -2.0
+ @test div_drift(sys, x_i) == -4.0
+end

test/largedeviations/gMAM.jl

@@ -1,3 +1,16 @@
+@testset "gMAM FitzHugh-Nagumo" begin
+ p = [0.1,3,1,1,1,0]
+ fhn = CoupledSDEs(fitzhugh_nagumo, diag_noise_funtion(0.1), zeros(2), p)
+ x_i = SA[sqrt(2/3), sqrt(2/27)]
+ x_f = SA[0.001, 0.0]
+ N = 100
+ res = geometric_min_action_method(fhn, x_i, x_f; N=75, maxiter=200)
+ S = geometric_action(fhn, res[1][end])
+ @test isapprox(S, 0.18, atol=0.01)
+end
+
+
+"""
 @testset "gMAM Meier Stein" begin
  function meier_stein(u, p, t) # out-of-place
  x, y = u
@@ -36,3 +49,4 @@
  # @test all(isapprox.(action_val, 0.3375, atol = 1e-3))
  end # HeymannVandenEijnden
 end # gMAM Meier Stein
+"""

test/runtests.jl

@@ -18,10 +18,11 @@ end
 # include("stochsystem.jl")
 # end
 
-# @testset "Large Deviations" begin
-# include("largedeviations/MAM.jl")
-# include("largedeviations/gMAM.jl")
-# end
+@testset "Large Deviations" begin
+ include("largedeviations/action_fhn.jl")
+ include("largedeviations/MAM.jl")
+ include("largedeviations/gMAM.jl")
+end
 
 # @testset "utilities" begin
 # include("utils.jl")