refactor: simple geometric minimal action method

benchmark/KPO_sgMAM.jl

@@ -43,7 +43,7 @@ function H_p(x, p) # ℜ² → ℜ²
     return Matrix([H_pu H_pv]')
 end
 
-sys = SgmamSystem{false,2}(H_x, H_p)
+sys = ExtendedHamiltonianSystem{false,2}(H_x, H_p)
 
 # setup
 Nt = 500  # number of discrete time steps
@@ -58,7 +58,7 @@ xx = @. (xb[1] - xa[1]) * s + xa[1] + 4 * s * (1 - s) * xsaddle[1]
 yy = @. (xb[2] - xa[2]) * s + xa[2] + 4 * s * (1 - s) * xsaddle[2] + 0.01 * sin(2π * s)
 x_initial = Matrix([xx yy]')
 
-MLP = sgmam(sys, x_initial; iterations=100_000, ϵ=10e2, show_progress=true)
+MLP = simple_geometric_min_action_method(sys, x_initial; iterations=100_000, ϵ=10e2, show_progress=true)
 x_min = MLP.path
 S_min = MLP.action
 
@@ -69,7 +69,7 @@ plot(x_initial[1, :], x_initial[2, :]; label="init", lw=3, c=:black)
 plot!(x_min[1, :], x_min[2, :]; label="MLP", lw=3, c=:red)
 plot!(string[1, :], string[2, :]; label="string", lw=3, c=:blue)
 
-@btime $sgmam($sys, $x_initial, iterations=100, ϵ=10e2, show_progress=false) # 25.803 ms (29024 allocations: 105.69 MiB)
-@profview sgmam(sys, x_initial, iterations=100, ϵ=10e2, show_progress=false)
+@btime $simple_geometric_min_action_method($sys, $x_initial, iterations=100, ϵ=10e2, show_progress=false) # 25.803 ms (29024 allocations: 105.69 MiB)
+@profview simple_geometric_min_action_method(sys, x_initial, iterations=100, ϵ=10e2, show_progress=false)
 
 # The bottleneck is atm at the LinearSolve call to update the x in the new iteration. So the more improve, one needs to write it own LU factorization.

benchmark/evaluate_drift.jl

@@ -51,7 +51,7 @@ function H_p(x, p) # ℜ² → ℜ²
     return Matrix([H_pu H_pv]')
 end
 
-sys_m = SgmamSystem{false,2}(H_x, H_p)
+sys_m = ExtendedHamiltonianSystem{false,2}(H_x, H_p)
 
 x_init_m = Matrix([xx yy]')
 
@@ -79,7 +79,7 @@ ds = CoupledODEs(prob)
 jac = jacobian(ds)
 jac([1, 1], (), 0.0)
 
-sgSys′ = SgmamSystem(ds);
+sgSys′ = ExtendedHamiltonianSystem(ds);
 
 p_r = rand(2, Nt)
 sgSys′.H_x(x_init_m, p_r) ≈ sys_m.H_x(x_init_m, p_r)

benchmark/string_method.jl

@@ -51,7 +51,7 @@ function sss()
         return StateSpaceSet([H_pu H_pv])
     end
 
-    sys_sss = SgmamSystem(H_x, H_p)
+    sys_sss = ExtendedHamiltonianSystem(H_x, H_p)
 
     x_init_sss = StateSpaceSet([xx yy])
     return sys_sss, x_init_sss
@@ -75,7 +75,7 @@ function m()
         return Matrix([H_pu H_pv]')
     end
 
-    sys_m = SgmamSystem(H_x, H_p)
+    sys_m = ExtendedHamiltonianSystem(H_x, H_p)
 
     x_init_m = Matrix([xx yy]')
     return sys_m, x_init_m

docs/src/examples/sgMAM_KPO.md

@@ -63,10 +63,10 @@ function H_p(x, p) # ℜ² → ℜ²
     return Matrix([H_pu H_pv]')
 end
 
-sys = SgmamSystem{false,2}(H_x, H_p)
+sys = ExtendedHamiltonianSystem{false,2}(H_x, H_p)
 ````
 
-We saved this function in the `SgmamSystem` struct. We want to find the optimal path between two attractors in the phase space. We define the initial trajectory as `wiggle` between the two attractors.
+We saved this function in the `ExtendedHamiltonianSystem` struct. We want to find the optimal path between two attractors in the phase space. We define the initial trajectory as `wiggle` between the two attractors.
 
 ````@example sgMAM_KPO
 Nt = 500  # number of discrete time steps
@@ -85,10 +85,10 @@ yy = @. (xb[2] - xa[2]) * s + xa[2] + 4 * s * (1 - s) * xsaddle[2] + 0.01 * sin(
 x_initial = Matrix([xx yy]')
 ````
 
-The optimisation is the performed by the `sgmam` function:
+The optimisation is the performed by the `simple_geometric_min_action_method` function:
 
 ````@example sgMAM_KPO
-MLP = sgmam(sys, x_initial; iterations=1_000, ϵ=10e2, show_progress=false)
+MLP = simple_geometric_min_action_method(sys, x_initial; iterations=1_000, ϵ=10e2, show_progress=false)
 x_min = MLP.path;
 nothing #hide
 ````

docs/src/man/largedeviations.md

@@ -37,8 +37,8 @@ The simple gMAM reduces the complexity of the original gMAM by requiring only fi
 
 The implementation below perform a constrained gradient descent where it assumes an autonomous system with additive Gaussian noise.
 ```@docs
-sgmam
-SgmamSystem
+simple_geometric_min_action_method
+ExtendedHamiltonianSystem
 ```
 
 ## Action functionals

examples/sgMAM_KPO.jl

@@ -54,9 +54,9 @@ function H_p(x, p) # ℜ² → ℜ²
     return Matrix([H_pu H_pv]')
 end
 
-sys = SgmamSystem{false,2}(H_x, H_p)
+sys = ExtendedHamiltonianSystem{false,2}(H_x, H_p)
 
-# We saved this function in the `SgmamSystem` struct. We want to find the optimal path between two attractors in the phase space. We define the initial trajectory as `wiggle` between the two attractors.
+# We saved this function in the `ExtendedHamiltonianSystem` struct. We want to find the optimal path between two attractors in the phase space. We define the initial trajectory as `wiggle` between the two attractors.
 
 Nt = 500  # number of discrete time steps
 s = collect(range(0; stop=1, length=Nt))
@@ -70,9 +70,9 @@ xx = @. (xb[1] - xa[1]) * s + xa[1] + 4 * s * (1 - s) * xsaddle[1]
 yy = @. (xb[2] - xa[2]) * s + xa[2] + 4 * s * (1 - s) * xsaddle[2] + 0.01 * sin(2π * s)
 x_initial = Matrix([xx yy]')
 
-# The optimisation is the performed by the `sgmam` function:
+# The optimisation is the performed by the `simple_geometric_min_action_method` function:
 
-MLP = sgmam(sys, x_initial; iterations=1_000, ϵ=10e2, show_progress=false)
+MLP = simple_geometric_min_action_method(sys, x_initial; iterations=1_000, ϵ=10e2, show_progress=false)
 x_min = MLP.path;
 
 # The function returns the optimal path `x_min`, the minimal action `S_min`, the Lagrange multipliers `lambda` associated with the optimal path, the optimal generalised momentum `p`, and the time derivative of the optimal path `xdot`. We can plot the initial trajectory and the optimal path:

src/CriticalTransitions.jl

@@ -72,7 +72,7 @@ export CoupledSDEs, CoupledODEs, noise_process, covariance_matrix, diffusion_mat
 export drift, div_drift, solver
 export StateSpaceSet
 
-export sgmam, SgmamSystem
+export simple_geometric_min_action_method, ExtendedHamiltonianSystem
 export fw_action, om_action, action, geometric_action
 export min_action_method, action_minimizer, geometric_min_action_method, string_method
 export MinimumActionPath

src/largedeviations/sgMAM.jl

@@ -5,11 +5,11 @@ This system operates in an extended phase space where the Hamiltonian is assumed
 quadratic in the extended momentum. The phase space coordinates `x` are doubled to
 form a 2n-dimensional extended phase space.
 """
-struct SgmamSystem{IIP,D,Hx,Hp}
+struct ExtendedHamiltonianSystem{IIP,D,Hx,Hp}
     H_x::Hx
     H_p::Hp
 
-    function SgmamSystem(ds::ContinuousTimeDynamicalSystem)
+    function ExtendedHamiltonianSystem(ds::ContinuousTimeDynamicalSystem)
         if ds isa CoupledSDEs
             proper_sgMAM_system(ds)
         end
@@ -38,16 +38,16 @@ struct SgmamSystem{IIP,D,Hx,Hp}
         end
         return new{isinplace(ds),dimension(ds),typeof(H_x),typeof(H_p)}(H_x, H_p)
     end
-    function SgmamSystem{IIP,D}(H_x::Function, H_p::Function) where {IIP,D}
+    function ExtendedHamiltonianSystem{IIP,D}(H_x::Function, H_p::Function) where {IIP,D}
         return new{IIP,D,typeof(H_x),typeof(H_p)}(H_x, H_p)
     end
 end
 
-function prettyprint(mlp::SgmamSystem{IIP,D}) where {IIP,D}
+function prettyprint(mlp::ExtendedHamiltonianSystem{IIP,D}) where {IIP,D}
     return "Doubled $D-dimensional phase space containing $(IIP ? "in-place" : "out-of-place") functions"
 end
 
-Base.show(io::IO, mlp::SgmamSystem) = print(io, prettyprint(mlp))
+Base.show(io::IO, mlp::ExtendedHamiltonianSystem) = print(io, prettyprint(mlp))
 
 """
 $(TYPEDSIGNATURES)
@@ -67,8 +67,8 @@ the Lagrange multipliers, the momentum, and the state derivatives. The implement
 based on the work of [Grafke et al. (2019)](https://homepages.warwick.ac.uk/staff/T.Grafke/simplified-geometric-minimum-action-method-for-the-computation-of-instantons.html.
 ).
 """
-function sgmam(
-    sys::SgmamSystem,
+function simple_geometric_min_action_method(
+    sys::ExtendedHamiltonianSystem,
     x_initial::Matrix{T};
     ϵ::Real=1e-1,
     iterations::Int64=1000,
@@ -103,11 +103,11 @@ function sgmam(
         StateSpaceSet(x'), S[end]; λ=lambda, generalized_momentum=p, path_velocity=xdot
     )
 end
-function sgmam(sys, x_initial::StateSpaceSet; kwargs...)
-    return sgmam(sys, Matrix(Matrix(x_initial)'); kwargs...)
+function simple_geometric_min_action_method(sys, x_initial::StateSpaceSet; kwargs...)
+    return simple_geometric_min_action_method(sys, Matrix(Matrix(x_initial)'); kwargs...)
 end
-function sgmam(sys::ContinuousTimeDynamicalSystem, x_initial::Matrix{<:Real}; kwargs...)
-    return sgmam(SgmamSystem(sys), Matrix(Matrix(x_initial)'); kwargs...)
+function simple_geometric_min_action_method(sys::ContinuousTimeDynamicalSystem, x_initial::Matrix{<:Real}; kwargs...)
+    return simple_geometric_min_action_method(ExtendedHamiltonianSystem(sys), Matrix(Matrix(x_initial)'); kwargs...)
 end
 
 function init_allocation(x_initial, Nt)

src/largedeviations/string_method.jl

@@ -8,7 +8,7 @@ The string method is an iterative algorithm used to find minimum energy path (ME
 This implementation allows for computation between arbitrary points, not just stable fixed points.
 
 # Arguments
-- `sys::SgmamSystem`: The doubled phase space system for which the string method is computed
+- `sys::ExtendedHamiltonianSystem`: The doubled phase space system for which the string method is computed
 - `x_initial`: Initial path discretized as a matrix where each column represents a point on the path
 - `ϵ::Real`: Step size for the evolution step
 - `iterations::Int64`: Maximum number of iterations for path convergence
@@ -18,7 +18,7 @@ This implementation allows for computation between arbitrary points, not just st
 - `x`: The final converged path representing the MEP
 """
 function string_method(
-    sys::Union{SgmamSystem,Function},
+    sys::Union{ExtendedHamiltonianSystem,Function},
     x_initial::Matrix;
     ϵ::Real=1e-1,
     iterations::Int64=1000,
@@ -68,7 +68,7 @@ function string_method(sys::ContinuousTimeDynamicalSystem, init; kwargs...)
 end
 
 function string_method(
-    b::Union{SgmamSystem,Function},
+    b::Union{ExtendedHamiltonianSystem,Function},
     x_initial::StateSpaceSet{D};
     ϵ::Real=1e-1,
     iterations::Int64=1000,
@@ -90,10 +90,10 @@ function string_method(
     return x
 end
 
-function update_x!(x::Matrix, sys::SgmamSystem, ϵ::Real)
+function update_x!(x::Matrix, sys::ExtendedHamiltonianSystem, ϵ::Real)
     return x += ϵ * sys.H_p(x, 0 * x) # euler integration
 end
-function update_x!(x::StateSpaceSet, sys::SgmamSystem, ϵ::Real)
+function update_x!(x::StateSpaceSet, sys::ExtendedHamiltonianSystem, ϵ::Real)
     return x += ϵ .* vec(sys.H_p(x, 0 * Matrix(x))) # euler integration
 end
 function update_x!(x::StateSpaceSet, b::Function, ϵ::Real)

test/largedeviations/API.jl

@@ -17,6 +17,6 @@
     for inst in [corr_alt, addit_non_autom, linear_multipli]
         @test_throws ArgumentError min_action_method(inst, init, T)
         @test_throws ArgumentError geometric_min_action_method(inst, init)
-        @test_throws ArgumentError sgmam(inst, init)
+        @test_throws ArgumentError simple_geometric_min_action_method(inst, init)
     end
 end

test/largedeviations/sgMAM.jl

@@ -2,7 +2,7 @@ using CriticalTransitions
 using ModelingToolkit
 using Test
 
-@testset "SgmamSystem KPO" begin
+@testset "ExtendedHamiltonianSystem KPO" begin
     λ = 3 / 1.21 * 2 / 295
     ω0 = 1.000
     ω = 1.000
@@ -51,8 +51,8 @@ using Test
     prob = ODEProblem(sysMTK, sts .=> zeros(2), (0.0, 100.0), (); jac=true)
     ds = CoupledODEs(prob)
 
-    sys = SgmamSystem{false,2}(H_x, H_p)
-    sys′ = SgmamSystem(ds)
+    sys = ExtendedHamiltonianSystem{false,2}(H_x, H_p)
+    sys′ = ExtendedHamiltonianSystem(ds)
 
     Nt = 500  # number of discrete time steps
     p_r = rand(2, Nt)

test/largedeviations/string_method.jl

@@ -50,7 +50,7 @@ yy = @. (xb[2] - xa[2]) * s + xa[2] + 4 * s * (1 - s) * xsaddle[2] + 0.01 * sin(
         return StateSpaceSet([H_pu H_pv])
     end
 
-    sys_sss = SgmamSystem{false,2}(H_x, H_p)
+    sys_sss = ExtendedHamiltonianSystem{false,2}(H_x, H_p)
 
     x_init_sss = StateSpaceSet([xx yy])
 
@@ -75,7 +75,7 @@ yy = @. (xb[2] - xa[2]) * s + xa[2] + 4 * s * (1 - s) * xsaddle[2] + 0.01 * sin(
         return Matrix([H_pu H_pv]')
     end
 
-    sys_m = SgmamSystem{false,2}(H_x, H_p)
+    sys_m = ExtendedHamiltonianSystem{false,2}(H_x, H_p)
 
     x_init_m = Matrix([xx yy]')