Revert "Fix depricated hysteris_sweep file (#110)"

This reverts commit b43fbfa.

src/HarmonicBalance.jl

@@ -17,6 +17,7 @@ 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))
@@ -51,18 +52,17 @@ 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,17 +106,13 @@ $(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]; dvar = d(get_variables(eom), tvar);
-    lhs = getfield.(eom.equations, :lhs); rhs = getfield.(eom.equations, :rhs);
-
-    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)])
+    tvar = get_independent_variables(eom)[1]
 
-    # Hopf-containing equations or MF equation are arranged by construstion
-    HB_bool || hopf_bool || MF_bool
+    # 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
 
+
 """
 $(TYPEDSIGNATURES)
 Rearrange `eom` to the standard form, such that the derivatives of the variables are on one side.

src/hysteresis_sweep.jl

@@ -0,0 +1,101 @@
+export plot_1D_solutions_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
+    physical     = filter_solutions.(
+                            filter_solutions.(
+                            res.solutions, res.classes["physical"]),res.classes["stable"])
+
+    steadystates = reduce(hcat,physical[index]) #change structure for easier indexing
+    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="sqrt(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)
+
+    #filter stable solutions
+    Y  = transform_solutions(res, y)
+    Ys = filter_solutions.(filter_solutions.(Y, res.classes["physical"]),res.classes["stable"]);
+    Ys = real.(reduce(hcat,Ys)) #facilitate indexing. Ys is an array (length(sweep))x(#solutions)
+
+    followed_branch    = zeros(Int64,length(Y)) #followed branch indexes
+    followed_branch[1] = starting_branch
+
+    if sweep == "left"
+        Ys = reverse(Ys,dims=2)
+    end
+
+    p1 = collect(keys(res.swept_parameters))[1] #parameter values
+
+    for i in 2:size(Ys)[2]
+        s = Ys[followed_branch[i-1],i] #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
+            if sweep == "right"
+                next_index = i
+            elseif sweep == "left" #reverse order of iterators
+                next_index = length(Y)-i+1
+            end
+
+            # 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)
+            print("bifurcation @ ", p1 ," = ", real(sol_dict[p1])," ")
+            sol_dict   = Dict(zip(keys(sol_dict),real.(values(sol_dict)) .+ 0.0im .+ ϵ*rand(length(values(sol_dict)))))
+            problem_t  = TimeEvolution.ODEProblem(res.problem.eom, steady_solution=sol_dict, timespan=(0,tf))
+            res_t      = TimeEvolution.solve(problem_t,saveat=tf)
+
+            followed_branch[i] = _closest_branch_index(res,res_t.u[end],next_index) #closest branch to final state
+
+            print("switched branch ", followed_branch[i-1] ," -> ", followed_branch[i],"\n")
+        end
+    end
+
+    if sweep == "left"
+        Ys = reverse(Ys,dims=2)
+        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, xscale="linear",yscale="linear",plot_only=["physical"],marker_classification="stable",filename=nothing,kwargs...)
+    plot_1D_solutions(res; x=x, y=y, xscale=xscale,yscale=yscale,
+                        plot_only=plot_only,marker_classification=marker_classification,filename=filename,kwargs...)
+
+    followed_branch,Ys = follow_branch(starting_branch,res,y=y,sweep=sweep,tf=tf,ϵ=ϵ)
+    if sweep == "left"
+        m = "<"
+    elseif sweep == "right"
+        m = ">"
+    end
+
+    Y_followed = [Ys[branch,param_idx] for (param_idx,branch) in enumerate(followed_branch)]
+
+    lines = plot(transform_solutions(res, x),Y_followed,c="k",marker=m; _set_Plots_default..., kwargs...)
+
+    #extract plotted data and return it
+    xdata,ydata = [line.get_xdata() for line in lines], [line.get_ydata() for line in lines]
+    markers = [line.get_marker() for line in lines]
+    save_dict= Dict(string(x) => xdata,string(y)=>ydata,"markers"=>markers)
+
+    !isnothing(filename) ? JLD2.save(_jld2_name(filename), save_dict) : nothing
+    return save_dict
+end

src/modules/LinearResponse/Lorentzian_spectrum.jl

@@ -67,23 +67,21 @@ _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, 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)])
+function JacobianSpectrum(res::Result; index::Int, branch::Int)
 
-    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
+    res.classes["stable"][index][branch] || error("\nThe solution is unstable - it has no JacobianSpectrum!\n")
 
     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
 
@@ -132,3 +130,4 @@ function _simplify_spectra!(spectra::Dict{Num, JacobianSpectrum})
         spectra[var] = _simplify_spectrum(spectra[var])
     end
 end
+

src/modules/LinearResponse/plotting.jl

@@ -1,6 +1,9 @@
 using Plots, Latexify, ProgressMeter
-using HarmonicBalance: _set_Plots_default
-export plot_linear_response, plot_rotframe_jacobian_response
+using Latexify.LaTeXStrings
+import ..HarmonicBalance: dim, _get_mask
+using Plots.PlotUtils.ColorSchemes
+
+export plot_linear_response
 
 
 function get_jacobian_response(res::Result, nat_var::Num, Ω_range, branch::Int; show_progress=true)
@@ -20,21 +23,6 @@ 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)
@@ -107,23 +95,8 @@ 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
 
-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
+heatmap(X, Ω_range,  C; color=:viridis,
+    xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), HarmonicBalance._set_Plots_default..., kwargs...)
 end
 
 """
@@ -150,7 +123,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("Ω"), _set_Plots_default..., kwargs...)
+       xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), HarmonicBalance._set_Plots_default..., kwargs...)
 end
 
 function plot_eigenvalues(res; branch, type=:imag, projection= v -> 1, cscheme = :default, kwargs...)

src/modules/LinearResponse/response.jl

@@ -1,4 +1,7 @@
 export get_response
+export plot_response
+
+
 
 """
     get_response_matrix(diff_eq::DifferentialEquation, freq::Num; order=2)
@@ -92,3 +95,8 @@ 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,9 +2,7 @@ module TimeEvolution
 
     using ..HarmonicBalance
     using Symbolics
-    using Plots
     using OrdinaryDiffEq
-    using OrderedCollections
     using DSP
     using FFTW
     using Peaks
@@ -14,6 +12,5 @@ module TimeEvolution
     include("TimeEvolution/ODEProblem.jl")
     include("TimeEvolution/FFT_analysis.jl")
     include("TimeEvolution/sweeps.jl")
-    include("TimeEvolution/hysteresis_sweep.jl")
 
 end

src/plotting_Plots.jl

@@ -94,14 +94,6 @@ 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,7 +112,6 @@ 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