Workshop fixes

Project.toml

@@ -1,7 +1,7 @@
 name = "HarmonicBalance"
 uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 authors = ["Jan Kosata <[email protected]>", "Javier del Pino <[email protected]>"]
-version = "0.6.0"
+version = "0.6.1"
 
 [deps]
 BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"
@@ -36,8 +36,8 @@ Latexify = "0.15.16"
 Peaks = "0.4.0"
 Plots = "1.35.0"
 ProgressMeter = "1.7.2"
-Symbolics = "4.10.4"
-julia = "1.8.0"
+Symbolics = "4.11.0"
+julia = "1.8.2"
 
 [extras]
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

src/HarmonicBalance.jl

@@ -56,7 +56,7 @@ module HarmonicBalance
 
     include("modules/TimeEvolution.jl")
     using .TimeEvolution
-    export ParameterSweep
+    export ParameterSweep, ODEProblem, solve
 
     include("modules/LimitCycles.jl")
     using .LimitCycles

src/hysteresis_sweep.jl

@@ -89,7 +89,7 @@ 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; 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]

src/modules/LinearResponse/plotting.jl

@@ -63,7 +63,7 @@ C = order == 1 ? get_jacobian_response(res, nat_var, Ω_range, branch) : get_lin
 C = logscale ? log.(C) : C
 
 heatmap(X, Ω_range,  C; color=:viridis, 
-    xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), kwargs...)
+    xlabel=latexify(string(first(keys(res.swept_parameters)))), ylabel=latexify("Ω"), HarmonicBalance._set_Plots_default..., kwargs...)
 end
 
 

src/modules/TimeEvolution/ODEProblem.jl

@@ -1,6 +1,8 @@
 using LinearAlgebra, Latexify
-import HarmonicBalance: transform_solutions, plot
-export transform_solutions, plot
+import HarmonicBalance: transform_solutions, plot, plot!
+import DifferentialEquations: ODEProblem, solve
+export transform_solutions, plot, plot!
+export ODEProblem, solve
 
 """ 
 
@@ -86,14 +88,21 @@ Plot a function `f` of a time-dependent solution `soln` of `harm_eq`.
 
 Parametric plot of f[1] against f[2]
 
+Also callable as plot!
 """
-function plot(soln::OrdinaryDiffEq.ODECompositeSolution, funcs, harm_eq::HarmonicEquation; kwargs...)
-    HarmonicBalance._set_Plots_default()
+function plot(soln::OrdinaryDiffEq.ODECompositeSolution, funcs, harm_eq::HarmonicEquation; add=false, kwargs...)
+
+    # start a new plot if needed
+    p = add ? plot!() : plot()
+
     if funcs isa String || length(funcs) == 1
-        plot(soln.t, transform_solutions(soln, funcs, harm_eq); legend=false, xlabel="time", ylabel=latexify(funcs), kwargs...)
+        plot!(soln.t, transform_solutions(soln, funcs, harm_eq); HarmonicBalance._set_Plots_default..., xlabel="time", ylabel=latexify(funcs), legend=false, kwargs...)
     elseif length(funcs) == 2 # plot of func vs func
-        plot(transform_solutions(soln, funcs, harm_eq)...; legend=false, xlabel=latexify(funcs[1]), ylabel=latexify(funcs[2]), kwargs...)
+        plot!(transform_solutions(soln, funcs, harm_eq)...; HarmonicBalance._set_Plots_default..., xlabel=latexify(funcs[1]), ylabel=latexify(funcs[2]), legend=false, kwargs...)
     else
         error("Invalid plotting argument: ", funcs)
     end
-end
+end
+
+
+plot!(soln::OrdinaryDiffEq.ODECompositeSolution, varargs...; kwargs...) = plot(soln, varargs...; add=true, kwargs...)

src/plotting_Plots.jl

@@ -1,11 +1,14 @@
 using Plots, Latexify
 import Plots.plot, Plots.plot!; export plot, plot!, plot_phase_diagram, savefig
 
-function _set_Plots_default()
-    font = "Computer Modern"
-    default(linewidth=2, legend=:outerright)
-    default(fontfamily=font, titlefont=font , tickfont=font)
-end
+const _set_Plots_default = Dict{Symbol, Any}([
+    :fontfamily => "computer modern",
+    :titlefont => "computer modern",
+    :tickfont => "computer modern",
+    :linewidth => 2,
+    :legend => :outerright])
+
+
 
 dim(res::Result) = length(size(res.solutions)) # give solution dimensionality
 
@@ -40,11 +43,10 @@ default behaviour: plot stable solutions as full lines, unstable as dashed
 the x and y axes are taken automatically from `res`
 """
 function plot(res::Result, varargs...; kwargs...)::Plots.Plot
-    HarmonicBalance._set_Plots_default()
     if dim(res) == 1
-        plot1D(res, varargs...; kwargs...)
+        plot1D(res, varargs...; _set_Plots_default..., kwargs...)
     elseif dim(res) == 2
-        plot2D(res, varargs...; kwargs...)
+        plot2D(res, varargs...; _set_Plots_default..., kwargs...)
     else
         error("Data dimension ", dim(res), " not supported")
     end
@@ -56,9 +58,9 @@ $(TYPEDSIGNATURES)
 
 Similar to `plot` but adds a plot onto an existing plot.
 """
-plot!(res::Result, varargs...; kwargs...)::Plots.Plot = plot(res, varargs...; add=true, kwargs...)
-
-
+function plot!(res::Result, varargs...; kwargs...)::Plots.Plot
+    plot(res, varargs...; add=true, _set_Plots_default..., kwargs...)
+end
 """
 $(TYPEDSIGNATURES)
 
@@ -105,7 +107,7 @@ function plot1D(res::Result; x::String="default", y::String, class="default", no
 
     if class == "default"
         if not_class == [] # plot stable full, unstable dashed
-            p = plot1D(res; x=x, y=y, class=["physical", "stable"], kwargs...)
+            p = plot1D(res; x=x, y=y, class=["physical", "stable"], add=add, kwargs...)
             plot1D(res; x=x, y=y, class="physical", not_class="stable", add=true, style=:dash, kwargs...)
             return p
         else
@@ -167,11 +169,10 @@ Class selection done by passing `String` or `Vector{String}` as kwarg:
 Other kwargs are passed onto Plots.gr()
 """
 function plot_phase_diagram(res::Result; kwargs...)::Plots.Plot
-    _set_Plots_default()
     if dim(res) == 1
-        plot_phase_diagram_1D(res; kwargs...)
+        plot_phase_diagram_1D(res; _set_Plots_default..., kwargs...)
     elseif dim(res) == 2
-        plot_phase_diagram_2D(res; kwargs...)
+        plot_phase_diagram_2D(res; _set_Plots_default..., kwargs...)
     else
         error("Data dimension ", dim(res), " not supported")
     end

src/solve_homotopy.jl

@@ -7,7 +7,7 @@ $(TYPEDSIGNATURES)
 Return an ordered dictionary specifying all variables and parameters of the solution
 in `result` on `branch` at the position `index`.
 """
-function get_single_solution(res::Result; branch::Int64, index)
+function get_single_solution(res::Result; branch::Int64, index)::OrderedDict{Num, ComplexF64}
 
     # check if the dimensionality of index matches the solutions
     if length(size(res.solutions)) !== length(index)
@@ -127,6 +127,7 @@ end
 
 get_steady_states(p::Problem, swept, fixed; kwargs...) = get_steady_states(p, ParameterRange(swept), ParameterList(fixed); kwargs...)
 get_steady_states(eom::HarmonicEquation, swept, fixed; kwargs...) = get_steady_states(Problem(eom), swept, fixed; kwargs...)
+get_steady_states(p::Union{Problem, HarmonicEquation}, fixed; kwargs...) = get_steady_states(p, [], fixed; kwargs...)
 
 
 """ Compile the Jacobian from `prob`, inserting `fixed_parameters`.
@@ -239,7 +240,7 @@ function _get_raw_solution(problem::Problem, parameter_values::Array{ParameterVe
     # HomotopyContinuation accepts 1D arrays of parameter sets
     params_1D = reshape(parameter_values, :, 1)
 
-    if random_warmup
+    if random_warmup && !isempty(sweep)
         complex_pert = [1E-2 * issubset(p, keys(sweep))*randn(ComplexF64) for p in problem.parameters] # complex perturbation of the warmup parameters
         warmup_parameters = params_1D[Int(round(length(params_1D)/2))] .* (ones(length(params_1D[1])) + complex_pert)
         warmup_solution = HomotopyContinuation.solve(problem.system,  target_parameters=warmup_parameters, threading=threading)