Hopf support (#29)

* implicit Jacobian removes gauge DoF

* substitute for HarmonicVariable fixed

* duplicites in Jacobians removed

* Symbolics overload fixes

* Hopf support completed

* Jacobian keywords unified

* save/load J fixed

Author: jkosata

Project.toml

@@ -1,7 +1,7 @@
 name = "HarmonicBalance"
 uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 authors = ["Jan Kosata <kosataj@phys.ethz.ch>", "Javier del Pino <jdelpino@phys.ethz.ch>"]
-version = "0.5.0"
+version = "0.5.1"
 
 [deps]
 BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"

src/HarmonicEquation.jl

@@ -2,7 +2,7 @@ export get_independent_variables
 export get_harmonic_equations
 export is_rearranged
 export slow_flow, slow_flow!
-export _equations_without_brackets
+export _remove_brackets
 
 show(eom::HarmonicEquation) = show_fields(eom)
 
@@ -105,7 +105,9 @@ Check if `eom` is rearranged to the standard form, such that the derivatives of
 """
 function is_rearranged(eom::HarmonicEquation)
     tvar = get_independent_variables(eom)[1]
-    isequal(getfield.(eom.equations, :rhs), d(get_variables(eom), tvar))
+
+    # Hopf-containing equations are arranged by construstion (impossible to check)
+    isequal(getfield.(eom.equations, :rhs), d(get_variables(eom), tvar)) || in("Hopf", getfield.(eom.variables, :type))
 end
 
 
@@ -142,10 +144,11 @@ end
 ###
 
 "Apply `rules` to both `equations` and `variables` field of `eom`"
-function substitute_all!(eom::HarmonicEquation, rules::Dict)
-    eom.equations = replace(eom.equations, rules)
-    eom.variables = substitute_all(eom.variables, rules)
-end
+function substitute_all(eom::HarmonicEquation, rules::Union{Dict, Pair}) 
+    new_eom = deepcopy(eom)
+    new_eom.equations = expand_derivatives.(substitute_all(eom.equations, rules))
+    new_eom
+end 
 
 
 "Simplify the equations in HarmonicEquation."
@@ -264,9 +267,14 @@ _set_zero_rhs(eq::Equation) = eq.lhs - eq.rhs ~ 0
 _set_zero_rhs(eqs::Vector{Equation}) = [_set_zero_rhs(eq) for eq in eqs]
 
 
+_remove_brackets(var::Num) = declare_variable(var_name(var))
+_remove_brackets(vars::Vector{Num}) = _remove_brackets.(vars)
+_remove_brackets(hv::HarmonicVariable) = _remove_brackets(hv.symbol)
+
+
 "Returns the equation system in `eom`, dropping all argument brackets (i.e., u(T) becomes u)."
-function _equations_without_brackets(eom::HarmonicEquation)
-    variable_rules = [var => declare_variable(var_name(var)) for var in get_variables(eom)]
+function _remove_brackets(eom::HarmonicEquation)
+    variable_rules = [var => _remove_brackets(var) for var in get_variables(eom)]
     equations_lhs = Num.(getfield.(eom.equations, :lhs)  - getfield.(eom.equations, :rhs))
     substitute_all(equations_lhs, variable_rules)
 end

src/HarmonicVariable.jl

@@ -1,4 +1,5 @@
-import Symbolics.get_variables; export get_variables
+import Symbolics: get_variables; export get_variables
+import Base.isequal; export isequal
 
 # pretty-printing
 display(var::HarmonicVariable) = display(var.name)
@@ -27,6 +28,10 @@ end
 # Functions for variable substutions and manipulation of HarmonicVariable
 ###
 
+
+# when HV is used for substitute, substitute its symbol
+substitute_all(eq::Union{Num, Equation}, rules::Dict{HarmonicVariable}) = substitute(eq, Dict(zip(getfield.(keys(rules), :symbol), values(rules))))
+
 function substitute_all(var::HarmonicVariable, rules)
     sym, freq = var.symbol, var.ω
     HarmonicVariable(substitute_all(sym, rules), var.name, var.type, substitute_all(freq, rules), var.natural_variable)
@@ -41,6 +46,8 @@ get_variables(vars::Vector{Num}) = unique(flatten([Num.(get_variables(x)) for x
 
 get_variables(var::HarmonicVariable) = Num.(get_variables(var.symbol))
 
+isequal(v1::HarmonicVariable, v2::HarmonicVariable) = isequal(v1.symbol, v2.symbol)
+
 
 
 

src/modules/HC_wrapper/homotopy_interface.jl

@@ -1,3 +1,4 @@
+using LinearAlgebra
 import HomotopyContinuation: Variable
 import HarmonicBalance: Problem; export Problem
 export Num_to_Variable
@@ -36,19 +37,21 @@ declare_variable(x::Num) = declare_variable(string(x))
 
 "Constructor for the type `Problem` (to be solved by HomotopyContinuation)
 from a `HarmonicEquation`."
-function Problem(eom::HarmonicEquation; explicit_Jacobian=true)
+function Problem(eom::HarmonicEquation; Jacobian=true)
 
     S = System(eom)
     # use the rearranged system for the proper definition of the Jacobian
     # this possibly has variables in the denominator and cannot be used for solving
-    if explicit_Jacobian == true
+    if Jacobian == true || Jacobian == "explicit"
         J = HarmonicBalance.get_Jacobian(eom)
-    elseif explicit_Jacobian isa Matrix
-        size(explicit_Jacobian)[1] == size(explicit_Jacobian)[2] == length(get_variables(eom)) || error("The Jacobian shape must match the number of variables!")
-        J = explicit_Jacobian
+    elseif Jacobian == "implicit"
+        # compute the Jacobian implicitly
+        J = HarmonicBalance.LinearResponse.get_implicit_Jacobian(eom)
+    elseif Jacobian == "false"
+        dummy_J(arg) = I(1)
+        J = dummy_J
     else
-        # ignore the jacobian, is later computed implicitly
-        J = false
+        J = Jacobian
     end
     vars_orig  = get_variables(eom)
     vars_new = declare_variable.(HarmonicBalance.var_name.(vars_orig))
@@ -69,7 +72,7 @@ end
 
 
 function System(eom::HarmonicEquation)
-    eqs = expand_derivatives.(_equations_without_brackets(eom))
+    eqs = expand_derivatives.(_remove_brackets(eom))
     conv_vars = Num_to_Variable.(get_variables(eom))
     conv_para = Num_to_Variable.(eom.parameters)
     S = HomotopyContinuation.System(parse_equations(eqs),variables=conv_vars,parameters=conv_para)

src/modules/Hopf/Hopf.jl

@@ -2,6 +2,9 @@ export add_Hopf!
 export Hopf_variables
 export replace_Hopf_variable!
 
+using HarmonicBalance: is_rearranged, rearrange_standard
+using HarmonicBalance: LinearResponse.get_implicit_Jacobian
+
 
 "Add a Hopf bifurcation to the system."
 function add_Hopf!(eom::DifferentialEquation, Δω::Num)
@@ -16,39 +19,56 @@ end
 add_Hopf!(eom::DifferentialEquation; frequency::Num, multiplicity::Int) = [add_Hopf!(eom, frequency) for k in 1:multiplicity]
 
 
-
 """
     Hopf_variables(eom, freq)
 """
 Hopf_variables(eom::HarmonicEquation, freq::Num) = filter(x -> any(isequal.(freq, get_all_terms(x.ω))), eom.variables)
 
 
-function replace_Hopf_variable!(eom::HarmonicEquation, freq::Num)
-    pair_to_change = Hopf_variables(eom, freq)[1] # eliminate of the Cartesian variables, it does not matter from which VDP pair
-    to_eliminate = pair_to_change.symbols[2]
-    rules = Dict(freq => HarmonicBalance.declare_variable(string(freq), first(get_independent_variables(eom))), to_eliminate => Num(0))
-    eom.equations = expand_derivatives.(substitute_all(eom.equations, rules))
-    eom.parameters = setdiff(eom.parameters, [freq]) # freq is now NOT a parameter
-    eom.variables = setdiff(eom.variables, [pair_to_change])
-    new_variable = HarmonicBalance.declare_variable(string(freq))
-    eom.variables = cat(eom.variables, _VDP_Hopf(pair_to_change, to_eliminate, new_variable), dims=1)
-end
+"""
+    Obtain the Jacobian of `eom` with a free (Hopf) variable `free_var`.
+    `free_var` marks the variable which is free due to U(1) symmetry. Its entry in the Jacobian matrix must be discarded
+    (the discarded eigenvalue is 0, corresponding to a free phase.)
+"""
+function _Hopf_Jacobian(eom::HarmonicEquation, free_var::HarmonicVariable)
+    eom_Jac = rearrange_standard(eom)
+    free_idx = findall(x -> isequal(x, free_var), eom_Jac.variables)  
+    deleteat!(eom_Jac.equations, free_idx)
+    deleteat!(eom_Jac.variables, free_idx)
 
-function replace_Hopf_variable(eom::HarmonicEquation, freq::Num)
-    new_eom = deepcopy(eom)
-    replace_Hopf_variable!(new_eom, freq)
-    new_eom
+    # the free variable can be fixed to zero, this is also done in the corresponding HarmonicEquation later
+    eom_Jac = substitute_all(eom_Jac, free_var => 0)
+    J = get_implicit_Jacobian(eom_Jac)
 end
 
 
-function _VDP_Hopf(old::HarmonicVariable, eliminated::Num, new_sym::Num)
-    new_VDP = deepcopy(old)
-    idx = findall(x-> isequal(x, eliminated), old.symbols)[1]
-    indep = Num(first(arguments(eliminated.val)))
-    new_VDP.symbols[idx] = HarmonicBalance.declare_variable(string(new_sym), indep)
-    new_VDP.types[idx] = "Hopf"
-    new_VDP.names[new_sym] = string(new_sym)
-    delete!(new_VDP.names, Num(var_name(eliminated)))
-    return new_VDP
+"""
+    Construct a `Problem` in the case where U(1) symmetry is present
+    due to having added a Hopf variable with frequency `Hopf_ω`.
+"""
+function _Hopf_Problem(eom::HarmonicEquation, Hopf_ω::Num)
+
+    isempty(Hopf_variables(eom, Hopf_ω)) ? error("No Hopf variables found!") : nothing
+    !any(isequal.(eom.parameters, Hopf_ω)) ? error(Hopf_ω, " is not a parameter of the harmonic equation!") : nothing 
+
+    free_var = Hopf_variables(eom, Hopf_ω)[end] # eliminate one of the Cartesian variables, it does not matter which
+    
+    # get the Hopf Jacobian before altering anything - this is the usual Jacobian but the entries corresponding
+    # to the free variable are removed
+    J = _Hopf_Jacobian(eom, free_var)
+    
+    new_symbol = HarmonicBalance.declare_variable(string(Hopf_ω), first(get_independent_variables(eom)))
+    rules = Dict(Hopf_ω => new_symbol, free_var.symbol => Num(0))
+    eom.equations = expand_derivatives.(substitute_all(eom.equations, rules))
+    eom.parameters = setdiff(eom.parameters, [Hopf_ω]) # Hopf_ω is now NOT a parameter anymore
+    
+    free_var.type = "Hopf"
+    free_var.ω = Num(0)
+    free_var.symbol = new_symbol
+    free_var.name = var_name(new_symbol)
+
+    # define Problem as usual but with the Hopf Jacobian (always computed implicitly)
+    p = Problem(eom; Jacobian=J)
+    return p
 end
 

src/modules/LinearResponse/jacobian_spectrum.jl

@@ -76,7 +76,8 @@ function JacobianSpectrum(res::Result; index::Int, branch::Int)
     hvars = res.problem.eom.variables # fetch the vector of HarmonicVariable
     λs, vs = eigen(res.jacobian(solution_dict))
 
-    solution_dict = get_single_solution(res, index=index, branch=branch)
+    # convert OrderedDict to Dict - see Symbolics issue #601
+    solution_dict = Dict(get_single_solution(res, index=index, branch=branch))
 
     # blank JacobianSpectrum for each variable
     all_spectra = Dict{Num, JacobianSpectrum}([[nvar, JacobianSpectrum([])] for nvar in getfield.(hvars, :natural_variable)])

src/modules/LinearResponse/perturbations.jl

@@ -4,47 +4,30 @@ export get_Jacobian
 $(SIGNATURES)
 
 Obtain the symbolic Jacobian matrix of `eom` (either a `HarmonicEquation` or a `DifferentialEquation`).
+This is the linearised left-hand side of F(u) = du/dT.
 
 """
 function get_Jacobian(eom::HarmonicEquation)
 
-    lhs = Num.(getfield.(HarmonicBalance.rearrange_standard(eom).equations, :lhs))
-
-    old_vars = get_variables(eom)
-    bracket_rules = [var => HarmonicBalance.declare_variable(var_name(var)) for var in old_vars]
-
-    lhs = substitute_all(lhs,bracket_rules)
-    vars = getindex.(bracket_rules, 2)
-    J = Matrix{Num}(undef, length(vars), length(vars))
-    for idx in CartesianIndices(J)
-        J[idx] = expand_derivatives(d(lhs[idx[1]], vars[idx[2]]))
-    end
-    J
+    rearr = !HarmonicBalance.is_rearranged(eom) ? HarmonicBalance.rearrange_standard(eom) : eom
+    lhs = _remove_brackets(rearr)
+    vars = _remove_brackets.(eom.variables)
 
+    get_Jacobian(lhs, vars)
+
 end
 
-#   Obtain a Jacobian from a `DifferentialEquation` by first converting it into a `HarmonicEquation`.
+" Obtain a Jacobian from a `DifferentialEquation` by first converting it into a `HarmonicEquation`. "
 function get_Jacobian(diff_eom::DifferentialEquation)
     @variables T
     harmonic_eq = get_harmonic_equations(diff_eom, slow_time=T, fast_time=first(get_independent_variables(diff_eom)))
     get_Jacobian(harmonic_eq)
 end
 
 
-"""Convert a set of equations in variables `vars` into matrix form.
-If the equations are linear, this matrix does not depends on the variables."""
-function equations_to_matrix(eqs::Vector{Num}, vars::Vector{Num})
-    length(eqs) != length(vars) ? error("System under- or over-constrained!") : nothing
-    M = Matrix{Num}(undef, length(eqs), length(eqs))
-    for idx in CartesianIndices(M)
-        M[idx] = expand_derivatives(d(eqs[idx[1]], vars[idx[2]]))
-    end
-    M
-end
-
-# get the Jacobian of a set of equations `eqs` with respect to the variables `vars`
+" Get the Jacobian of a set of equations `eqs` with respect to the variables `vars`. "
 function get_Jacobian(eqs::Vector{Num}, vars::Vector{Num})
-    length(eqs) == length(vars) || error("Jacobians are only defined for square matrices!")
+    length(eqs) == length(vars) || error("Jacobians are only defined for square systems!")
     M = Matrix{Num}(undef, length(vars), length(vars))
 
     for idx in CartesianIndices(M)
@@ -55,17 +38,25 @@ end
 
 get_Jacobian(eqs::Vector{Equation}, vars::Vector{Num}) = get_Jacobian(Num.(getfield.(eqs, :lhs) .- getfield.(eqs, :rhs)), vars)
 
+
+# the Jacobian of some equations again
 function get_Jacobian_steady(eom::HarmonicEquation; differential_order=0)
+    Hopf_vars = first.(getfield.(filter(x -> x.type == "Hopf", eom.variables), :symbol))
+    Hopf_idx = findall(x -> any(isequal.(x, Hopf_vars)) , get_variables(eom))
+    nonsingular = filter( x -> x ∉ Hopf_idx, 1:length(get_variables(eom)))
+
     vars_simp = Dict([var => HarmonicBalance.declare_variable(var_name(var)) for var in get_variables(eom)])
     T = get_independent_variables(eom)[1]
-    J = get_Jacobian(eom.equations, d(get_variables(eom), T, differential_order))
+    J = get_Jacobian(eom.equations, d(get_variables(eom), T, differential_order))[nonsingular, nonsingular]
+    
     expand_derivatives.(HarmonicBalance.substitute_all(J, vars_simp))
 end
 
 # COMPILE THIS!
-function get_implicit_Jacobian(eom::HarmonicEquation)
+function get_implicit_Jacobian(eom::HarmonicEquation)::Function
     M = get_Jacobian_steady(eom, differential_order=0)
     Mp = get_Jacobian_steady(eom, differential_order=1)
+
     function J(soln::OrderedDict)
         -inv(ComplexF64.(substitute_all(Mp, soln))) * ComplexF64.(substitute_all(M, soln))
     end

src/saving.jl

@@ -16,8 +16,10 @@ end
 
 function save(filename, x::Result)
     x_nofunc = deepcopy(x)
+
     # compiled functions cause problems in saving: ignore J now, compile when loading
-    x_nofunc.jacobian = false
+    null_J() = 0
+    x_nofunc.jacobian = null_J
     JLD2.save(_jld2_name(filename), Dict("object" => x_nofunc))
 end
 

src/solve_homotopy.jl

@@ -98,6 +98,7 @@ function get_steady_states(prob::Problem, swept_parameters::ParameterRange, fixe
     # extract all the information we need from results
     #rounded_solutions = unique_points.(HomotopyContinuation.solutions.(getindex.(raw, 1)); metric = EuclideanNorm(), atol=1E-14, rtol=1E-8)
     rounded_solutions = HomotopyContinuation.solutions.(getindex.(raw, 1));
+    all(isempty.(rounded_solutions)) ? error("No solutions found!") : nothing
     solutions = pad_solutions(rounded_solutions)
 
     compiled_J = _compile_Jacobian(prob, swept_parameters, unique_fixed)
@@ -125,10 +126,10 @@ get_steady_states(eom::HarmonicEquation, swept, fixed; kwargs...) = get_steady_s
 function _compile_Jacobian(prob::Problem, swept_parameters::ParameterRange, fixed_parameters::ParameterList)
     J_variables = cat(prob.variables, collect(keys(swept_parameters)), dims=1)
 
-    if prob.jacobian != false
+    if prob.jacobian isa Matrix
         compiled_J = compile_matrix(prob.jacobian, J_variables, fixed_parameters)
     else
-        compiled_J = LinearResponse.get_implicit_Jacobian(prob.eom)
+        compiled_J = prob.jacobian # leave implicit Jacobian as is
     end
     compiled_J
 end

src/types.jl

@@ -136,7 +136,7 @@ end
 
 """Gives the relation between `var` and the underlying natural variable."""
 function _show_ansatz(var::HarmonicVariable)
-    terms = Dict("u" => "cos", "v" => "sin", "a" => "")
+    terms = Dict("u" => "cos", "v" => "sin", "a" => "", Hopf => "")
     t = var.natural_variable.val.arguments
     t = length(t) == 1 ? string(t[1]) : error("more than 1 independent variable")
     ω = string(var.ω)
@@ -168,8 +168,9 @@ $(TYPEDFIELDS)
 #  Constructors
 ```julia
 Problem(eom::HarmonicEquation; Jacobian=true) # find and store the symbolic Jacobian
-Problem(eom::HarmonicEquation; Jacobian=false) # ignore the Jacobian for now
-Problem(eom::HarmonicEquation; Jacobian::Matrix) # use the given matrix as the Jacobian
+Problem(eom::HarmonicEquation; Jacobian="implicit") # ignore the Jacobian for now, compute implicitly later
+Problem(eom::HarmonicEquation; Jacobian=J) # use J as the Jacobian (a function that takes a Dict)
+Problem(eom::HarmonicEquation; Jacobian=false) # ignore the Jacobian
 ```
 """
 mutable struct Problem
@@ -222,7 +223,7 @@ mutable struct Result
     If problem.jacobian is a symbolic matrix, this holds a compiled function.
     If problem.jacobian was `false`, this holds a function that rearranges the equations to find J
     only after numerical values are inserted (preferable in cases where the symbolic J would be very large)."
-    jacobian
+    jacobian::Function
 
     Result(sol,swept, fixed, problem, classes, J) = new(sol, swept, fixed, problem, classes, J)
     Result(sol,swept, fixed, problem, classes) = new(sol, swept, fixed, problem, classes)