Fix Limit cycle bug (#186)

* Fix symbolic version

* increase timout-minutes of CI

* format limit_cycle module

* get rid of unsued kwarg `explicit_Jacobian` in `_cycle_Problem`

* change `cycle_harmonic` kwarg to `limit_cycle_harmonic`

* instead of sorting always take the first in the list

* clean test file

src/modules/LimitCycles/analysis.jl

@@ -1,6 +1,6 @@
 import HarmonicBalance: classify_solutions, _free_symbols, _symidx, _is_physical
 
-function classify_unique!(res::Result, Δω; class_name="unique_cycle")
+function classify_unique!(res::Result, Δω; class_name = "unique_cycle")
 
     # 1st degeneracy: arbitrary sign of Δω
     i1 = _symidx(Δω, res)
@@ -12,4 +12,4 @@ function classify_unique!(res::Result, Δω; class_name="unique_cycle")
     c2 = classify_solutions(res, soln -> _is_physical(soln) && real(soln[i2]) >= 0)
 
     res.classes[class_name] = map(.*, c1, c2)
-end
+end

src/modules/LimitCycles/gauge_fixing.jl

@@ -3,8 +3,8 @@ export add_pairs!
 
 using HarmonicBalance: is_rearranged, rearrange_standard, _remove_brackets
 using HarmonicBalance.LinearResponse: get_implicit_Jacobian, get_Jacobian
-import HarmonicBalance: is_stable, is_physical, is_Hopf_unstable, order_branches!, classify_binaries!, find_branch_order
-
+import HarmonicBalance: is_stable, is_physical, is_Hopf_unstable, order_branches!,
+                        classify_binaries!, find_branch_order
 
 function add_pairs!(eom::DifferentialEquation, ω_lc::Num)
     for var in get_variables(eom), ω in eom.harmonics[var]
@@ -13,27 +13,24 @@ function add_pairs!(eom::DifferentialEquation, ω_lc::Num)
     end
 end
 
-
 """
     add_pairs!(eom::DifferentialEquation; ω_lc::Num, n=1)
 
 Add a limit cycle harmonic `ω_lc` to the system
 Equivalent to adding `n` pairs of harmonics ω +- ω_lc for each existing ω.
 """
-add_pairs!(eom::DifferentialEquation; ω_lc::Num, n::Int) = [add_pairs!(eom, ω_lc) for k in 1:n]
-
+add_pairs!(eom::DifferentialEquation; ω_lc::Num, n::Int) = [add_pairs!(eom, ω_lc)
+                                                            for k in 1:n]
 
 """
 $(TYPEDSIGNATURES)
 
 Return the harmonic variables which participate in the limit cycle labelled by `ω_lc`.
 """
-function get_cycle_variables(eom::HarmonicEquation, ω_lc::Num; sorted=false)
+function get_cycle_variables(eom::HarmonicEquation, ω_lc::Num)
     vars = filter(x -> any(isequal.(ω_lc, get_all_terms(x.ω))), eom.variables)
-    sorted ? sort(vars, by = x -> x.ω / ω_lc) : vars
 end
 
-
 """
     Obtain the Jacobian of `eom` with a gauge-fixed variable `fixed_var`.
     `fixed_var` marks the variable which is fixed to zero due to U(1) symmetry.
@@ -42,29 +39,29 @@ end
     For limit cycles, we always use an 'implicit' Jacobian - a function which only returns the numerical Jacobian when a numerical solution
     is inserted. Finding the analytical Jacobian is usually extremely time-consuming.
 """
-function _gaugefixed_Jacobian(eom::HarmonicEquation, fixed_var::HarmonicVariable; explicit=false, sym_order, rules)
+function _gaugefixed_Jacobian(eom::HarmonicEquation, fixed_var::HarmonicVariable;
+        explicit = false, sym_order, rules)
     if explicit
         _fix_gauge!(get_Jacobian(eom), fixed_var) # get a symbolic explicit J, compile later
     else
         rules = Dict(rules)
         setindex!(rules, 0, _remove_brackets(fixed_var))
-        get_implicit_Jacobian(eom, rules=rules, sym_order=sym_order)
+        get_implicit_Jacobian(eom, rules = rules, sym_order = sym_order)
     end
 end
 
-
 """
 $(TYPEDSIGNATURES)
 
 Construct a `Problem` from `eom` in the case where U(1) symmetry is present
 due to having added a limit cycle frequency `ω_lc`.
 `explicit_Jacobian=true` attempts to derive a symbolic Jacobian (usually not possible).
 """
-function _cycle_Problem(eom::HarmonicEquation, ω_lc::Num; explicit_Jacobian=false)
-
+function _cycle_Problem(eom::HarmonicEquation, ω_lc::Num)
     eom = deepcopy(eom) # do not mutate eom
     isempty(get_cycle_variables(eom, ω_lc)) ? error("No Hopf variables found!") : nothing
-    !any(isequal.(eom.parameters, ω_lc)) ? error(ω_lc, " is not a parameter of the harmonic equation!") : nothing
+    !any(isequal.(eom.parameters, ω_lc)) ?
+    error(ω_lc, " is not a parameter of the harmonic equation!") : nothing
 
     # eliminate one of the u,v variables (gauge fixing)
     fixed_var = _choose_fixed(eom, ω_lc) # remove the HV of the lowest subharmonic
@@ -74,11 +71,14 @@ function _cycle_Problem(eom::HarmonicEquation, ω_lc::Num; explicit_Jacobian=fal
     _fix_gauge!(eom, ω_lc, fixed_var)
 
     # define Problem as usual but with the Hopf Jacobian (always computed implicitly)
-    p = Problem(eom; Jacobian=nothing)
+    p = Problem(eom; Jacobian = nothing)
     return p
 end
 
-_choose_fixed(eom, ω_lc) = get_cycle_variables(eom, ω_lc, sorted=true)[1]
+function _choose_fixed(eom, ω_lc)
+    vars = get_cycle_variables(eom, ω_lc)
+    first(vars) # This is arbitrary; better would be to substitute with values
+end
 
 """
     get_limit_cycles(eom::HarmonicEquation, swept, fixed, ω_lc; kwargs...)
@@ -87,44 +87,56 @@ Variant of `get_steady_states` for a limit cycle problem characterised by a Hopf
 
 Solutions with ω_lc = 0 are labelled unphysical since this contradicts the assumption of distinct harmonic variables corresponding to distinct harmonics.
 """
-function get_limit_cycles(eom::HarmonicEquation, swept, fixed, ω_lc; explicit_Jacobian=false, kwargs...)
-    prob = _cycle_Problem(eom, ω_lc, explicit_Jacobian=explicit_Jacobian);
-    prob.jacobian = _gaugefixed_Jacobian(eom, _choose_fixed(eom, ω_lc), explicit=explicit_Jacobian, sym_order=_free_symbols(prob, swept), rules=fixed)
-    result = get_steady_states(prob, swept, fixed; method=:warmup, threading=true, classify_default=true, kwargs...)
+function get_limit_cycles(
+        eom::HarmonicEquation, swept, fixed, ω_lc; explicit_Jacobian = false, kwargs...)
+
+    prob = _cycle_Problem(eom, ω_lc)
+    prob.jacobian = _gaugefixed_Jacobian(
+        eom, _choose_fixed(eom, ω_lc), explicit = explicit_Jacobian,
+        sym_order = _free_symbols(prob, swept), rules = fixed)
+
+    result = get_steady_states(prob, swept, fixed; method = :warmup,
+        threading = true, classify_default = true, kwargs...)
 
     _classify_limit_cycles!(result, ω_lc)
 
     result
 end
 
-get_limit_cycles(eom::HarmonicEquation,swept,fixed; cycle_harmonic, kwargs...) = get_limit_cycles(eom, swept, fixed, cycle_harmonic; kwargs...)
+function get_limit_cycles(eom::HarmonicEquation, swept, fixed; limit_cycle_harmonic, kwargs...)
+    get_limit_cycles(eom, swept, fixed, limit_cycle_harmonic; kwargs...)
+end
 
 # if abs(ω_lc) < tol, set all classifications to false
 # TOLERANCE HARDCODED FOR NOW
 function _classify_limit_cycles!(res::Result, ω_lc::Num)
     ω_lc_idx = findfirst(x -> isequal(x, ω_lc), res.problem.variables)
-    for idx in CartesianIndices(res.solutions), c in filter(x -> x != "binary_labels", keys(res.classes))
-        res.classes[c][idx] .*= abs.(getindex.(res.solutions[idx], ω_lc_idx)) .> 1E-10
+    for idx in CartesianIndices(res.solutions),
+        c in filter(x -> x != "binary_labels", keys(res.classes))
+
+        res.classes[c][idx] .*= abs.(getindex.(res.solutions[idx], ω_lc_idx)) .> 1e-10
     end
 
     classify_unique!(res, ω_lc)
 
-    unique_stable = find_branch_order(map(.*, res.classes["stable"], res.classes["unique_cycle"]))
+    unique_stable = find_branch_order(map(
+        .*, res.classes["stable"], res.classes["unique_cycle"]))
 
     # branches which are unique but never stable
-    unique_unstable = setdiff(find_branch_order(map(.*, res.classes["unique_cycle"], res.classes["physical"])), unique_stable)
+    unique_unstable = setdiff(
+        find_branch_order(map(.*, res.classes["unique_cycle"], res.classes["physical"])),
+        unique_stable)
     order_branches!(res, vcat(unique_stable, unique_unstable)) # shuffle to have relevant branches first
 end
 
-
 """
 $(TYPEDSIGNATURES)
 
 Fix the gauge in `eom` where `ω_lc` is the limit cycle frequency by constraining `fixed_var` to zero and promoting `ω_lc` to a variable.
 """
 function _fix_gauge!(eom::HarmonicEquation, ω_lc::Num, fixed_var::HarmonicVariable)
-
-    new_symbol = HarmonicBalance.declare_variable(string(ω_lc), first(get_independent_variables(eom)))
+    new_symbol = HarmonicBalance.declare_variable(
+        string(ω_lc), first(get_independent_variables(eom)))
     rules = Dict(ω_lc => new_symbol, fixed_var.symbol => Num(0))
     eom.equations = expand_derivatives.(substitute_all(eom.equations, rules))
     eom.parameters = setdiff(eom.parameters, [ω_lc]) # ω_lc is now NOT a parameter anymore

test/limit_cycle.jl

@@ -2,24 +2,27 @@ using HarmonicBalance
 import HarmonicBalance.LinearResponse.plot_linear_response
 import HarmonicBalance.LimitCycles.get_limit_cycles
 
-@variables ω_lc, t, ω0, x(t), μ
 
-natural_equation = d(d(x, t), t) - μ * (1 - x^2) * d(x, t) + x
-dEOM = DifferentialEquation(natural_equation, x)
+@testset "van der Pol oscillator " begin
+    @variables ω_lc, t, ω0, x(t), μ
 
-# order should NOT matter if correct gauge is fixed (that corresponding to 1*ω_lc)
-add_harmonic!(dEOM, x, [ω_lc, 3 * ω_lc])
+    natural_equation = d(d(x, t), t) - μ * (1 - x^2) * d(x, t) + x
+    dEOM = DifferentialEquation(natural_equation, x)
 
-harmonic_eq = get_harmonic_equations(dEOM)
+    # order matters for 1*ω_lc gauge to be fixed
+    add_harmonic!(dEOM, x, [ω_lc, 3 * ω_lc])
 
-fixed = ();
-varied = μ => range(1, 5, 5)
+    harmonic_eq = get_harmonic_equations(dEOM)
+    HarmonicBalance.LimitCycles._choose_fixed(harmonic_eq, ω_lc)
 
-result = get_limit_cycles(harmonic_eq, varied, fixed, cycle_harmonic=ω_lc, show_progress=false, seed=SEED)
+    fixed = ();
+    varied = μ => range(1, 5, 5)
 
-@test sum(any.(classify_branch(result, "stable"))) == 4
-@test sum(any.(classify_branch(result, "unique_cycle"))) == 1
+    result = get_limit_cycles(harmonic_eq, varied, fixed, ω_lc; show_progress=false, seed=SEED)
 
-plot(result, y="ω_lc")
+    @test sum(any.(classify_branch(result, "stable"))) == 4
+    @test sum(any.(classify_branch(result, "unique_cycle"))) == 1
 
-plot_linear_response(result, x, branch=1, Ω_range=range(0.9, 1.1, 2), order=1)
+    plot(result, y="ω_lc")
+    plot_linear_response(result, x, branch=1, Ω_range=range(0.9, 1.1, 2), order=1)
+end