Fix depricated hysteris_sweep file (#110)

* update deprecated `hysteresis_sweep file`

the function was using undefined (old) function and assumed old data structures

add tests for hysteresis_sweep functions

* move `hystesis_sweep` to TimeEvolution Module

* Add extra tests for `hysteris_sweep`

* clean linear response test file

* plot jacobian response of a parameter ODE sweep

* change print statement to logging event `@info`
logging event can be turned off on demand. This behavior is favorable when working on clusters etc.

* fix bug where the new dict is not ordered

* add OrderedCollections to the TimeEvolution module

* only perturb the point in phase space;
not phase space itself

* add constructor for HarmonicEquation

* Make `is_rearranged` compatible with by constructed arranged harmonic equations

* add the option to plot empty spectra

* get rid of `'` typo in test file

---------

Co-authored-by: Javier del Pino <[email protected]>

src/HarmonicBalance.jl

@@ -17,7 +17,6 @@ module HarmonicBalance
     export @variables
     export d
 
-
     import Base: ComplexF64, Float64; export ComplexF64, Float64
     ComplexF64(x::Complex{Num}) = ComplexF64(Float64(x.re) + im*Float64(x.im))
     Float64(x::Complex{Num}) = Float64(ComplexF64(x))
@@ -52,17 +51,18 @@ module HarmonicBalance
     include("saving.jl")
     include("transform_solutions.jl")
     include("plotting_Plots.jl")
-    include("hysteresis_sweep.jl")
 
     include("modules/HC_wrapper.jl")
     using .HC_wrapper
 
     include("modules/LinearResponse.jl")
     using .LinearResponse
+    export plot_linear_response, plot_rotframe_jacobian_response
 
     include("modules/TimeEvolution.jl")
     using .TimeEvolution
     export ParameterSweep, ODEProblem, solve
+    export plot_1D_solutions_branch, follow_branch
 
     include("modules/LimitCycles.jl")
     using .LimitCycles

src/HarmonicEquation.jl

@@ -106,12 +106,16 @@ $(TYPEDSIGNATURES)
 Check if `eom` is rearranged to the standard form, such that the derivatives of the variables are on one side.
 """
 function is_rearranged(eom::HarmonicEquation)
-    tvar = get_independent_variables(eom)[1]
+    tvar = get_independent_variables(eom)[1]; dvar = d(get_variables(eom), tvar);
+    lhs = getfield.(eom.equations, :lhs); rhs = getfield.(eom.equations, :rhs);
 
-    # 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
+    HB_bool = isequal(rhs, dvar)
+    hopf_bool = in("Hopf", getfield.(eom.variables, :type))
+    MF_bool = !any([occursin(str1,str2) for str1 in string.(dvar) for str2 in string.(lhs)])
 
+    # Hopf-containing equations or MF equation are arranged by construstion
+    HB_bool || hopf_bool || MF_bool
+end
 
 """
 $(TYPEDSIGNATURES)

src/modules/LinearResponse/Lorentzian_spectrum.jl

@@ -67,21 +67,23 @@ _get_as(hvars::Vector{HarmonicVariable}) = findall(x -> isequal(x.type, "a"), hv
 
 
 #   Returns the spectra of all variables in `res` for `index` of `branch`.
-function JacobianSpectrum(res::Result; index::Int, branch::Int)
+function JacobianSpectrum(res::Result; index::Int, branch::Int, force=false)
+    hvars = res.problem.eom.variables # fetch the vector of HarmonicVariable
+    # blank JacobianSpectrum for each variable
+    all_spectra = Dict{Num, JacobianSpectrum}([[nvar, JacobianSpectrum([])] for nvar in getfield.(hvars, :natural_variable)])
 
-    res.classes["stable"][index][branch] || error("\nThe solution is unstable - it has no JacobianSpectrum!\n")
+    if force
+        res.classes["stable"][index][branch] || return all_spectra # if the solution is unstable, return empty spectra
+    else
+        res.classes["stable"][index][branch] || error("\nThe solution is unstable - it has no JacobianSpectrum!\n")
+    end
 
     solution_dict = get_single_solution(res, branch=branch, index=index)
-
-    hvars = res.problem.eom.variables # fetch the vector of HarmonicVariable
     λs, vs = eigen(res.jacobian(solution_dict))
 
     # 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)])
-
     for (j, λ) in enumerate(λs)
         eigvec = vs[:, j] # the eigenvector
 
@@ -130,4 +132,3 @@ function _simplify_spectra!(spectra::Dict{Num, JacobianSpectrum})
         spectra[var] = _simplify_spectrum(spectra[var])
     end
 end
-

src/modules/LinearResponse/plotting.jl

@@ -1,9 +1,6 @@
 using Plots, Latexify, ProgressMeter
-using Latexify.LaTeXStrings
-import ..HarmonicBalance: dim, _get_mask
-using Plots.PlotUtils.ColorSchemes
-
-export plot_linear_response
+using HarmonicBalance: _set_Plots_default
+export plot_linear_response, plot_rotframe_jacobian_response
 
 
 function get_jacobian_response(res::Result, nat_var::Num, Ω_range, branch::Int; show_progress=true)
@@ -23,6 +20,21 @@ function get_jacobian_response(res::Result, nat_var::Num, Ω_range, branch::Int;
     end
     C
 end
+function get_jacobian_response(res::Result, nat_var::Num, Ω_range, followed_branches::Vector{Int}; show_progress=true, force=false)
+
+    spectra = [JacobianSpectrum(res, branch=branch, index = i, force=force) for (i, branch) in pairs(followed_branches)]
+    C = Array{Float64, 2}(undef,  length(Ω_range), length(spectra))
+
+    if show_progress
+        bar = Progress(length(CartesianIndices(C)), 1, "Diagonalizing the Jacobian for each solution ... ", 50)
+    end
+    # evaluate the Jacobians for the different values of noise frequency Ω
+    for ij in CartesianIndices(C)
+        C[ij] = abs(evaluate(spectra[ij[2]][nat_var], Ω_range[ij[1]]))
+        show_progress ? next!(bar) : nothing
+    end
+    C
+end
 
 
 function get_linear_response(res::Result, nat_var::Num, Ω_range, branch::Int; order, show_progress=true)
@@ -95,8 +107,23 @@ X = Vector{Float64}(collect(values(res.swept_parameters))[1][stable])
 C = order == 1 ? get_jacobian_response(res, nat_var, Ω_range, branch, show_progress=show_progress) : get_linear_response(res, nat_var, Ω_range, branch; order=order, show_progress=show_progress)
 C = logscale ? log.(C) : C
 
-heatmap(X, Ω_range,  C; color=:viridis,
-    xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), HarmonicBalance._set_Plots_default..., kwargs...)
+xlabel = latexify(string(first(keys(res.swept_parameters)))); ylabel = latexify("Ω");
+heatmap(X, Ω_range,  C; color=:viridis, xlabel=xlabel, ylabel=ylabel, _set_Plots_default..., kwargs...)
+end
+function plot_linear_response(res::Result, nat_var::Num, followed_branches::Vector{Int}; Ω_range, logscale=false, show_progress=true, switch_axis = false, force=true, kwargs...)
+    length(size(res.solutions)) != 1 && error("1D plots of not-1D datasets are usually a bad idea.")
+
+    X = Vector{Float64}(collect(first(values(res.swept_parameters))))
+
+    C = get_jacobian_response(res, nat_var, Ω_range, followed_branches; force=force)
+    C = logscale ? log.(C) : C
+
+    xlabel = latexify(string(first(keys(res.swept_parameters)))); ylabel = latexify("Ω");
+    if switch_axis
+        heatmap(Ω_range, X, C'; color=:viridis, xlabel=ylabel, ylabel=xlabel, _set_Plots_default..., kwargs...)
+    else
+        heatmap(X, Ω_range, C; color=:viridis, xlabel=xlabel, ylabel=ylabel, _set_Plots_default..., kwargs...)
+    end
 end
 
 """
@@ -123,7 +150,7 @@ function plot_rotframe_jacobian_response(res::Result; Ω_range, branch::Int, log
     C = logscale ? log.(C) : C
 
     heatmap(X, Ω_range,  C; color=:viridis,
-       xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), HarmonicBalance._set_Plots_default..., kwargs...)
+       xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), _set_Plots_default..., kwargs...)
 end
 
 function plot_eigenvalues(res; branch, type=:imag, projection= v -> 1, cscheme = :default, kwargs...)

src/modules/LinearResponse/response.jl

@@ -1,7 +1,4 @@
 export get_response
-export plot_response
-
-
 
 """
     get_response_matrix(diff_eq::DifferentialEquation, freq::Num; order=2)
@@ -95,8 +92,3 @@ end
 # rotating frame into the lab frame
 _plusamp(uv) = norm(uv)^2 - 2*(imag(uv[1])*real(uv[2]) - real(uv[1])*imag(uv[2]))
 _minusamp(uv) = norm(uv)^2 + 2*(imag(uv[1])*real(uv[2]) - real(uv[1])*imag(uv[2]))
-
-
-
-
-

src/modules/TimeEvolution.jl

@@ -2,7 +2,9 @@ module TimeEvolution
 
     using ..HarmonicBalance
     using Symbolics
+    using Plots
     using OrdinaryDiffEq
+    using OrderedCollections
     using DSP
     using FFTW
     using Peaks
@@ -12,5 +14,6 @@ module TimeEvolution
     include("TimeEvolution/ODEProblem.jl")
     include("TimeEvolution/FFT_analysis.jl")
     include("TimeEvolution/sweeps.jl")
+    include("TimeEvolution/hysteresis_sweep.jl")
 
 end

src/modules/TimeEvolution/hysteresis_sweep.jl

@@ -0,0 +1,80 @@
+export plot_1D_solutions_branch, follow_branch
+
+"""Calculate distance between a given state and a stable branch"""
+function _closest_branch_index(res::Result, state::Vector{Float64}, index::Int64)
+    #search only among stable solutions
+    stable = HarmonicBalance._apply_mask(res.solutions, HarmonicBalance._get_mask(res, ["physical", "stable"], []))
+
+    steadystates = reduce(hcat, stable[index])
+    distances    = vec(sum(abs2.(steadystates .- state), dims=1))
+    return argmin(replace(distances, NaN => Inf))
+end
+
+"""Return the indexes and values following stable branches along a 1D sweep.
+When a no stable solutions are found (e.g. in a bifurcation), the next stable solution is calculated by time evolving the previous solution (quench).
+Keyword arguments
+- `y`:  Dependent variable expression (parsed into Symbolics.jl) to evaluate the followed solution branches on .
+- `sweep`: Direction for the sweeping of solutions. A `right` (`left`) sweep proceeds from the first (last) solution, ordered as the sweeping parameter.
+- `tf`: time to reach steady
+- `ϵ`: small random perturbation applied to quenched solution, in a bifurcation in order to favour convergence in cases where multiple solutions are identically accessible (e.g. symmetry breaking into two equal amplitude states)
+"""
+function follow_branch(starting_branch::Int64, res::Result; y="u1^2+v1^2", sweep="right", tf=10000, ϵ=1e-4)
+    sweep_directions = ["left", "right"]
+    sweep ∈ sweep_directions  || error("Only the following (1D) sweeping directions are allowed:  ", sweep_directions)
+
+    # get stable solutions
+    Y  = transform_solutions(res, y)
+    Ys = HarmonicBalance._apply_mask(Y, HarmonicBalance._get_mask(res, ["physical", "stable"], []))
+    Ys = sweep == "left" ? reverse(Ys) : Ys
+
+    followed_branch    = zeros(Int64,length(Y))  # followed branch indexes
+    followed_branch[1] = starting_branch
+
+    p1 = first(keys(res.swept_parameters)) # parameter values
+
+    for i in 2:length(Ys)
+        s = Ys[i][followed_branch[i-1]] # solution amplitude in the current branch and current parameter index
+        if !isnan(s) # the solution is not unstable or unphysical
+            followed_branch[i] = followed_branch[i-1]
+        else # bifurcation found
+            next_index = sweep == "right" ? i : length(Ys)-i+1
+
+            # create a synthetic starting point out of an unphysical solution: quench and time evolve
+            # the actual solution is complex there, i.e. non physical. Take real part for the quench.
+            sol_dict  = get_single_solution(res, branch=followed_branch[i-1], index=next_index)
+
+            var = res.problem.variables
+            var_values_noise = real.(getindex.(Ref(sol_dict), var)) .+ 0.0im .+ ϵ*rand(length(var))
+            for (i, v) in enumerate(var)
+                sol_dict[v] = var_values_noise[i]
+            end
+
+            problem_t = OrdinaryDiffEq.ODEProblem(res.problem.eom, sol_dict, timespan=(0, tf))
+            res_t     = OrdinaryDiffEq.solve(problem_t, OrdinaryDiffEq.Tsit5(), saveat=tf)
+
+            followed_branch[i] = _closest_branch_index(res, res_t.u[end], next_index) # closest branch to final state
+
+            @info "bifurcation @ $p1 = $(real(sol_dict[p1])): switched branch $(followed_branch[i-1]) ➡ $(followed_branch[i])"
+        end
+    end
+    if sweep == "left"
+        Ys = reverse(Ys)
+        followed_branch = reverse(followed_branch)
+    end
+
+    return followed_branch, Ys
+end
+
+
+"""1D plot with the followed branch highlighted"""
+function plot_1D_solutions_branch(starting_branch::Int64, res::Result;
+    x::String, y::String, sweep="right", tf=10000, ϵ=1e-4, class="default", not_class=[], kwargs...)
+    p = plot(res; x=x, y=y, class=class, not_class=not_class, kwargs...)
+
+    followed_branch, Ys = follow_branch(starting_branch, res, y=y, sweep=sweep, tf=tf, ϵ=ϵ)
+    Y_followed = [Ys[param_idx][branch] for (param_idx,branch) in enumerate(followed_branch)]
+    X = real.(res.swept_parameters[HarmonicBalance._parse_expression(x)])
+
+    Plots.plot!(p, X, real.(Y_followed); linestyle=:dash, c=:gray, label = sweep*" sweep", HarmonicBalance._set_Plots_default..., kwargs...)
+    return p
+end

src/plotting_Plots.jl

@@ -94,6 +94,14 @@ function _apply_mask(solns::Array{Vector{ComplexF64}},  booleans)
     factors = replace.(booleans, 0 => NaN)
     map(.*, solns, factors)
 end
+function _apply_mask(solns::Vector{Vector{Vector{ComplexF64}}}, booleans)
+    Nan_vector = NaN.*similar(solns[1][1])
+    new_solns = [
+        [booleans[i][j] ? solns[i][j] : Nan_vector for j in eachindex(solns[i])]
+        for i in eachindex(solns)
+    ]
+    return new_solns
+end
 
 
 """ Project the array `a` into the real axis, warning if its contents are complex. """

src/types.jl

@@ -112,6 +112,7 @@ mutable struct HarmonicEquation
 
     # use a self-referential constructor with _parameters
     HarmonicEquation(equations, variables, nat_eq) = (x = new(equations, variables, Vector{Num}([]), nat_eq); x.parameters=_parameters(x); x)
+    HarmonicEquation(equations, variables, parameters, natural_equation) = new(equations, variables, parameters, natural_equation)
 end