diff --git a/benchmark/KPO_sgMAM.jl b/benchmark/KPO_sgMAM.jl index 110108eb..4c6b01b6 100644 --- a/benchmark/KPO_sgMAM.jl +++ b/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. diff --git a/benchmark/evaluate_drift.jl b/benchmark/evaluate_drift.jl index 434eeb07..ac8adc44 100644 --- a/benchmark/evaluate_drift.jl +++ b/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) diff --git a/benchmark/string_method.jl b/benchmark/string_method.jl index 027e1467..a09f76e7 100644 --- a/benchmark/string_method.jl +++ b/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 diff --git a/docs/src/examples/sgMAM_KPO.md b/docs/src/examples/sgMAM_KPO.md index 8d65fb06..d28caac2 100644 --- a/docs/src/examples/sgMAM_KPO.md +++ b/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 ```` diff --git a/docs/src/man/largedeviations.md b/docs/src/man/largedeviations.md index 889903b7..e71db65e 100644 --- a/docs/src/man/largedeviations.md +++ b/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 diff --git a/examples/sgMAM_KPO.jl b/examples/sgMAM_KPO.jl index 4d51cffa..30ae9383 100644 --- a/examples/sgMAM_KPO.jl +++ b/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: diff --git a/src/CriticalTransitions.jl b/src/CriticalTransitions.jl index 0d77b7dd..56aa0be0 100644 --- a/src/CriticalTransitions.jl +++ b/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 diff --git a/src/largedeviations/sgMAM.jl b/src/largedeviations/sgMAM.jl index 4ea19be0..5ddabc89 100644 --- a/src/largedeviations/sgMAM.jl +++ b/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) diff --git a/src/largedeviations/string_method.jl b/src/largedeviations/string_method.jl index 772003d1..2e6f33c9 100644 --- a/src/largedeviations/string_method.jl +++ b/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) diff --git a/test/largedeviations/API.jl b/test/largedeviations/API.jl index 420f7edf..c1ba313c 100644 --- a/test/largedeviations/API.jl +++ b/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 diff --git a/test/largedeviations/sgMAM.jl b/test/largedeviations/sgMAM.jl index 2af7ba07..ce3a117d 100644 --- a/test/largedeviations/sgMAM.jl +++ b/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) diff --git a/test/largedeviations/string_method.jl b/test/largedeviations/string_method.jl index c4ce195e..2dd9c2ad 100644 --- a/test/largedeviations/string_method.jl +++ b/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]')