added keyword to control the display of the progress bar (resolves #60)

src/DifferentialEquation.jl

@@ -16,7 +16,7 @@ julia> add_harmonic!(diff_eq, x, ω) # expand x using ω
 
 System of 1 differential equations
 Variables:       x(t)
-Harmonic ansatz: x(t) => ω;   
+Harmonic ansatz: x(t) => ω;
 
 (ω0^2)*x(t) + Differential(t)(Differential(t)(x(t))) ~ F*cos(t*ω)
 ```
@@ -36,7 +36,7 @@ get_variables(diff_eom::DifferentialEquation) = collect(keys(diff_eom.equations)
 
 is_harmonic(diff_eom::DifferentialEquation, t::Num) = all([is_harmonic(eq, t) for eq in values(diff_eom.equations)])
 
-"Pretty printing of the newly defined types" 
+"Pretty printing of the newly defined types"
 function show_fields(object)
     for field in fieldnames(typeof(object)) # display every field
         display(string(field)); display(getfield(object, field))

src/HarmonicBalance.jl

@@ -32,7 +32,7 @@ module HarmonicBalance
     import Base: exp
     exp(x::Complex{Num}) = x.re.val == 0 ? exp(im*x.im.val) : exp(x.re.val + im*x.im.val)
 
-    include("types.jl")    
+    include("types.jl")
 
     include("utils.jl")
     include("Symbolics_customised.jl")

src/HarmonicEquation.jl

@@ -11,7 +11,7 @@ show(eom::HarmonicEquation) = show_fields(eom)
     harmonic_ansatz(eom::DifferentialEquation, time::Num; coordinates="Cartesian")
 
 Expand each variable of `diff_eom` using the harmonics assigned to it with `time` as the time variable.
-For each harmonic of each variable, instance(s) of `HarmonicVariable` are automatically created and named. 
+For each harmonic of each variable, instance(s) of `HarmonicVariable` are automatically created and named.
 
 """
 function harmonic_ansatz(diff_eom::DifferentialEquation, time::Num)
@@ -20,7 +20,7 @@ function harmonic_ansatz(diff_eom::DifferentialEquation, time::Num)
     rules, vars = Dict(), []
 
     # keep count to label new variables
-    uv_idx = 1 
+    uv_idx = 1
     a_idx = 1
 
     for nvar in get_variables(diff_eom) # sum over natural variables
@@ -39,12 +39,12 @@ function harmonic_ansatz(diff_eom::DifferentialEquation, time::Num)
             end
             to_substitute += rule
         end
-        
+
         rules[nvar] = to_substitute # total sub rule for nvar
     end
     eqs = substitute_all(eqs, rules)
     HarmonicEquation(eqs, Vector{HarmonicVariable}(vars), diff_eom)
-end 
+end
 
 
 function slow_flow!(eom::HarmonicEquation; fast_time::Num, slow_time::Num, degree=2)
@@ -126,7 +126,7 @@ end
 
 """
 $(TYPEDSIGNATURES)
-Get the internal symbols of the independent variables of `eom`. 
+Get the internal symbols of the independent variables of `eom`.
 """
 function get_variables(eom::HarmonicEquation)
     return flatten(get_variables.(eom.variables))
@@ -140,7 +140,7 @@ get_variables(p::Problem) = get_variables(p.eom)
 function _parameters(eom::HarmonicEquation)
     all_symbols = flatten([cat(get_variables(eq.lhs), get_variables(eq.rhs), dims=1) for eq in eom.equations])
     # subtract the set of independent variables (i.e., time) from all free symbols
-    setdiff(all_symbols, get_variables(eom), get_independent_variables(eom)) 
+    setdiff(all_symbols, get_variables(eom), get_independent_variables(eom))
 end
 
 
@@ -149,11 +149,11 @@ end
 ###
 
 "Apply `rules` to both `equations` and `variables` field of `eom`"
-function substitute_all(eom::HarmonicEquation, rules::Union{Dict, Pair}) 
+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 
+end
 
 
 "Simplify the equations in HarmonicEquation."
@@ -200,21 +200,21 @@ function fourier_transform!(eom::HarmonicEquation, time::Num)
         # "type" is usually "u" or "v" (harmonic) or ["a"] (zero-harmonic)
         if hvar.type == "u"
             avg_eqs[i] = fourier_cos_term(eq, hvar.ω, time)
-        elseif hvar.type == "v" 
+        elseif hvar.type == "v"
             avg_eqs[i] = fourier_sin_term(eq, hvar.ω, time)
         elseif hvar.type == "a"
             avg_eqs[i] = fourier_cos_term(eq, 0, time) # pick out the constants
         end
     end
-    
+
     eom.equations = avg_eqs
 end
 
 
 """
     get_harmonic_equations(diff_eom::DifferentialEquation; fast_time=nothing, slow_time=nothing)
 
-Apply the harmonic ansatz, followed by the slow-flow, Fourier transform and dropping 
+Apply the harmonic ansatz, followed by the slow-flow, Fourier transform and dropping
 higher-order derivatives to obtain
 a set of ODEs (the harmonic equations) governing the harmonics of `diff_eom`.
 
@@ -242,7 +242,7 @@ A set of 2 harmonic equations
 Variables: u1(T), v1(T)
 Parameters: ω0, ω, F
 
-Harmonic ansatz: 
+Harmonic ansatz:
 x(t) = u1*cos(ωt) + v1*sin(ωt)
 
 Harmonic equations:
@@ -263,7 +263,7 @@ function get_harmonic_equations(diff_eom::DifferentialEquation; fast_time=nothin
     eom = slow_flow(eom, fast_time=fast_time, slow_time=slow_time; degree=degree); # drop 2nd order time derivatives
     fourier_transform!(eom, fast_time); # perform averaging over the frequencies originally specified in dEOM
     ft_eom_simplified = drop_powers(eom, d(get_variables(eom), slow_time), 2); # drop higher powers of the first-order derivatives
-    return ft_eom_simplified 
+    return ft_eom_simplified
 end
 
 

src/HarmonicVariable.jl

@@ -8,7 +8,7 @@ display(var::Vector{HarmonicVariable}) = display.(getfield.(var, Symbol("name"))
 
 function _coordinate_transform(new_var, ω, t, type)
     coords = Dict([
-        "u" =>  new_var * cos(ω*t), 
+        "u" =>  new_var * cos(ω*t),
         "v" => new_var * sin(ω*t),
         "a" => new_var])
     return coords[type]

src/Symbolics_customised.jl

@@ -31,7 +31,7 @@ expand_all(x::Num) = Num(expand_all(x.val))
 
 "Apply a function f on every member of a sum or a product"
 _apply_termwise(f, x::Add) = sum([f(arg) for arg in arguments(x)])
-_apply_termwise(f, x::Mul) = prod([f(arg) for arg in arguments(x)])	
+_apply_termwise(f, x::Mul) = prod([f(arg) for arg in arguments(x)])
 _apply_termwise(f, x::Div) = _apply_termwise(f, x.num) / _apply_termwise(f, x.den)
 _apply_termwise(f,x) = f(x)
 
@@ -77,7 +77,7 @@ function exp_to_trig(x)
 
 			# put trigarg => -trigarg the lowest alphabetic argument of trigarg is lower than that of -trigarg
 			# this is a meaningless key but gives unique signs to all sums
-			is_first = minimum(string.(arguments(trigarg))) == first_symbol			
+			is_first = minimum(string.(arguments(trigarg))) == first_symbol
 			return is_first ? cos(-trigarg) -im*sin(-trigarg) : cos(trigarg)+im* sin(trigarg)
 		end
 

src/Symbolics_utils.jl

@@ -25,7 +25,7 @@ d(funcs::Vector{Num}, x::Num, deg=1) = [d(f, x, deg) for f in funcs]
  "Declare a variable that is a function of another variable in the the current namespace"
  function declare_variable(name::String, independent_variable::Num)
     # independent_variable = declare_variable(independent_variable) convert string into Num
-    var_sym = Symbol(name)  
+    var_sym = Symbol(name)
     new_var = @variables $var_sym(independent_variable)
     @eval($(var_sym) = first($new_var)) # store the variable under "name" in this namespace
     return eval(var_sym)
@@ -160,7 +160,7 @@ function is_harmonic(x::Num, t::Num)
         return false
     else
         powers = [max_power(first(term.val.arguments), t) for term in t_terms[trigs]]
-        return unique(powers) == [1] 
+        return unique(powers) == [1]
     end
 end
 
@@ -186,7 +186,7 @@ function trig_to_exp(x::Num)
         end
         # avoid Complex{Num} where possible as this causes bugs
         # instead, the Nums store SymbolicUtils Complex types
-        term = Num(term.re.val + im*term.im.val) 
+        term = Num(term.re.val + im*term.im.val)
         append!(rules, [trig => term])
     end
 
@@ -199,7 +199,7 @@ end
 
 "Return true if `f` is a function of `var`."
 is_function(f, var) = any(isequal.(get_variables(f), var))
- 
+
 
 "Return true if `f` is a sin or cos."
 function is_trig(f::Num)
@@ -261,7 +261,7 @@ end
 function max_power(x::Num, y::Num)
     terms = get_all_terms(x)
     powers = power_of.(terms, y)
-    maximum(powers)   
+    maximum(powers)
 end
 
 max_power(x::Vector{Num}, y::Num) = maximum(max_power.(x, y))

src/classification.jl

@@ -124,7 +124,7 @@ function classify_binaries!(res::Result)
     res.classes["binary_labels"] = bin_label
 end
 
-clean_bitstrings(res::Result) = [[el for el in bit_string[phys_string]] 
+clean_bitstrings(res::Result) = [[el for el in bit_string[phys_string]]
 for (bit_string,phys_string) in zip(res.classes["stable"],res.classes["physical"])]; #remove unphysical solutions from strings
 
 function bitarr_to_int(arr)
@@ -134,7 +134,7 @@ end
 
 """
 $(TYPEDSIGNATURES)
-Removes all solution branches from `res` where NONE of the solution falls into `class`. 
+Removes all solution branches from `res` where NONE of the solution falls into `class`.
 Typically used to filter out unphysical solutions to prevent huge file sizes.
 """
 function filter_result!(res::Result, class::String)

src/hysteresis_sweep.jl

@@ -9,18 +9,18 @@ function _closest_branch_index(res::Result,state::Vector{Float64},index::Int64)
                             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))    
+    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. 
+"""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 
+- `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)
@@ -34,14 +34,14 @@ function follow_branch(starting_branch::Int64,res::Result; y="sqrt(u1^2+v1^2)",s
 
     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] 
+
+    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]
@@ -56,21 +56,21 @@ function follow_branch(starting_branch::Int64,res::Result; y="sqrt(u1^2+v1^2)",s
             #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))))) 
+            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
 
@@ -89,8 +89,8 @@ function plot_1D_solutions_branch(starting_branch::Int64,res::Result; x::String,
 
     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...)   
-    
+    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]

src/modules/HC_wrapper/homotopy_interface.jl

@@ -55,15 +55,15 @@ function Problem(eom::HarmonicEquation; Jacobian=true)
     end
     vars_orig  = get_variables(eom)
     vars_new = declare_variable.(HarmonicBalance.var_name.(vars_orig))
-    return Problem(vars_new, eom.parameters, S, J,eom)   
+    return Problem(vars_new, eom.parameters, S, J,eom)
 end
 
 
 "A constructor for Problem from explicitly entered equations, variables and parameters."
 function Problem(equations::Vector{Num},variables::Vector{Num},parameters::Vector{Num})
     conv_vars = Num_to_Variable.(variables)
     conv_para = Num_to_Variable.(parameters)
-    
+
     eqs_HC=[Expression(eval(symbol)) for symbol in [Meta.parse(s) for s in [string(eq) for eq in equations]]] #note in polar coordinates there could be imaginary factors, requiring the extra replacement "I"=>"1im"
     system = HomotopyContinuation.System(eqs_HC, variables = conv_vars, parameters = conv_para)
     J = HarmonicBalance.get_Jacobian(equations,variables) #all derivatives are assumed to be in the left hand side;

src/modules/LimitCycles/Hopf.jl

@@ -34,7 +34,7 @@ get_Hopf_variables(eom::HarmonicEquation, Δω::Num) = filter(x -> any(isequal.(
 
 """
     Obtain the Jacobian of `eom` with a gauge-fixed (Hopf) variable `fixed_var`.
-    `fixed_var` marks the variable which is fixed to zero due to U(1) symmetry. 
+    `fixed_var` marks the variable which is fixed to zero due to U(1) symmetry.
     This leaves behind a finite degeneracy of solutions (see Chapter 6 of thesis).
 
     The Jacobian is always 'implicit' - a function which only returns the numerical Jacobian when a numerical solution
@@ -48,37 +48,37 @@ end
 """
 $(TYPEDSIGNATURES)
 
-Construct a `Problem` from `eom` in the case where U(1) symmetry is present 
+Construct a `Problem` from `eom` in the case where U(1) symmetry is present
 due to having added a Hopf variable with frequency `Δω`.
 """
 function _Hopf_Problem(eom::HarmonicEquation, Δω::Num; explicit_Jacobian=false)
 
     eom = deepcopy(eom) # do not mutate eom
     isempty(get_Hopf_variables(eom, Δω)) ? error("No Hopf variables found!") : nothing
-    !any(isequal.(eom.parameters, Δω)) ? error(Δω, " is not a parameter of the harmonic equation!") : nothing 
+    !any(isequal.(eom.parameters, Δω)) ? error(Δω, " is not a parameter of the harmonic equation!") : nothing
 
     # eliminate one of the Cartesian variables
-    fixed_var = get_Hopf_variables(eom, Δω)[1] 
-    
+    fixed_var = get_Hopf_variables(eom, Δω)[1]
+
     # get the Hopf Jacobian before altering anything - this is the usual Jacobian but the entries corresponding
     # to the free variable are removed
     J = explicit_Jacobian ? substitute_all(get_Jacobian(eom), _remove_brackets(fixed_var) => 0) : _Hopf_implicit_Jacobian(eom, fixed_var)
     _fix_gauge!(eom, Δω, fixed_var)
-    
+
     # define Problem as usual but with the Hopf Jacobian (always computed implicitly)
     p = Problem(eom; Jacobian=J)
     return p
 end
 
 
 """
-    get_steady_states(eom::HarmonicEquation, swept, fixed, Δω; kwargs...)   
+    get_steady_states(eom::HarmonicEquation, swept, fixed, Δω; kwargs...)
 
 Variant of `get_steady_states` for a limit cycle problem characterised by a Hopf frequency (usually called Δω)
 
 Solutions with Δω = 0 are labelled unphysical since this contradicts the assumption of distinct harmonic variables corresponding to distinct harmonics.
 """
-function get_steady_states(eom::HarmonicEquation, swept, fixed, Δω; explicit_Jacobian=false, kwargs...)   
+function get_steady_states(eom::HarmonicEquation, swept, fixed, Δω; explicit_Jacobian=false, kwargs...)
     prob = _Hopf_Problem(eom, Δω, explicit_Jacobian=explicit_Jacobian);
     result = HarmonicBalance.get_steady_states(prob, swept, fixed; random_warmup=true, threading=true, classify_default=true, kwargs...)
 
@@ -110,15 +110,15 @@ end
 """
 $(TYPEDSIGNATURES)
 
-Fix the gauge in `eom` where `Δω` is the limit cycle frequency by constraining `fixed_var` to zero and promoting `Δω` to a variable. 
+Fix the gauge in `eom` where `Δω` is the limit cycle frequency by constraining `fixed_var` to zero and promoting `Δω` to a variable.
 """
 function _fix_gauge!(eom::HarmonicEquation, Δω::Num, fixed_var::HarmonicVariable)
 
     new_symbol = HarmonicBalance.declare_variable(string(Δω), first(get_independent_variables(eom)))
     rules = Dict(Δω => new_symbol, fixed_var.symbol => Num(0))
     eom.equations = expand_derivatives.(substitute_all(eom.equations, rules))
     eom.parameters = setdiff(eom.parameters, [Δω]) # Δω is now NOT a parameter anymore
-    
+
     fixed_var.type = "Hopf"
     fixed_var.ω = Num(0)
     fixed_var.symbol = new_symbol

src/modules/LimitCycles/analysis.jl

@@ -1,5 +1,3 @@
-
-
 import HarmonicBalance.classify_solutions
 
 function classify_unique!(res::Result, Δω; class_name="unique")