Progress bars (#20)

* progress bars added

* Symbolics updated

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.4.3"
+version = "0.4.4"
 
 [deps]
 BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"
@@ -18,6 +18,7 @@ Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Peaks = "18e31ff7-3703-566c-8e60-38913d67486b"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
@@ -37,10 +38,11 @@ HomotopyContinuation = "2.6.4"
 JLD2 = "0.4.19"
 Latexify = "0.15.9"
 Peaks = "0.4.0"
+ProgressMeter = "1.7.2"
 PyCall = "1.93.0"
 PyPlot = "2.10.0"
 SymbolicUtils = "0.19.7"
-Symbolics = "4.3.1"
+Symbolics = "4.4.1"
 julia = "1.7.0"
 
 [extras]

src/HarmonicBalance.jl

@@ -8,6 +8,7 @@ module HarmonicBalance
     export d
     export plot
     using Symbolics
+    using ProgressMeter
     import SymbolicUtils: Term, Add, Div, Mul, Pow, Sym
     using DocStringExtensions
 

src/modules/LinearResponse.jl

@@ -7,6 +7,7 @@ module LinearResponse
     using ..HarmonicBalance
     using ..HC_wrapper
     using PyPlot
+    using ProgressMeter
     using DocStringExtensions
     using Latexify
 

src/modules/LinearResponse/jacobian_spectrum.jl

@@ -144,9 +144,12 @@ function plot_jacobian_spectrum(res::Result, nat_var::Num; Ω_range, branch::Int
     spectra = [JacobianSpectrum(res, branch=branch, index = i) for i in (1:length(res.solutions))[stability]]
     C = Array{Float64, 2}(undef,  length(Ω_range)-1, length(X)-1)
 
+    bar = Progress(length(CartesianIndices(C)), 1, "Diagonalizing the Jacobian for each solution ... ", 50)
     for ij in CartesianIndices(C)
         C[ij] = abs(evaluate(spectra[ij[2]][nat_var], Ω_range[ij[1]] - offset[ij[2]]))
+        next!(bar)
     end
+
     x_mat = x_scale * hcat([x*ones(length(Ω_range)) for x in X]...)
     y_mat = y_scale * hcat([Ω_range for j=1:length(X)]...)
     C = logscale ? log.(C) : C

src/modules/LinearResponse/linear_response.jl

@@ -79,7 +79,10 @@ function plot_response(res::Result, Ω_range; branch::Int, logscale=false)
     Y = Array{Float64, 2}(undef,  length(Ω_range), length(X))
 
     # this could be optimized by not grabbing the entire huge dictionary every time
+
+    bar = Progress(length(Y), 1, "Solving the linear response ODE for each solution and input frequency ... ", 50)
     for j in 1:(size(Y)[2])
+        next!(bar)
         s = get_single_solution(res, branch=branch, index=j)
         for i in 1:(size(Y)[1])
             Y[i,j] = get_response(response, s, reverse(Ω_range)[i])

src/solve_homotopy.jl

@@ -36,8 +36,7 @@ get_single_solution(res::Result, index) = [get_single_solution(res, index=index,
                         fixed_parameters::ParameterList;
                         random_warmup=false,
                             threading=false,
-                            sorting="hilbert",
-                            show_progress=true)
+                            sorting="hilbert")
 
 Solves `prob` over the ranges specified by `swept_parameters`, keeping `fixed_parameters` constant.
 `swept_parameters` accepts pairs mapping symbolic variables to arrays or `LinRange`.
@@ -48,7 +47,6 @@ Keyword arguments
 This takes far longer but can be more reliable.
 - `threading`: If `true`, multithreaded support is activated. The number of available threads is set by the environment variable `JULIA_NUM_THREADS`. 
 - `sorting`: the method used by `sort_solutions` to get continuous solutions branches.  The current options are `"hilbert"` (1D sorting along a Hilbert curve), `"nearest"` (nearest-neighbor sorting) and `"none"`.
-- `show_progress`: If `true`  Indicate whether a progress bar should be displayed (see [HomotopyContinuation.jl docs](https://www.juliahomotopycontinuation.org/HomotopyContinuation.jl/stable/solve/) 
 
 Example: solving a simple harmonic oscillator ``m \\ddot{x} + γ \\dot{x} + ω_0^2 x = F \\cos(ωt)``
 to obtain the response as a function of ``ω``
@@ -86,7 +84,7 @@ A steady state result for 1000 parameter points
 ```
 
 """
-function get_steady_states(prob::Problem, swept_parameters::ParameterRange, fixed_parameters::ParameterList; random_warmup=false, threading=false, sorting="nearest",show_progress=true)   
+function get_steady_states(prob::Problem, swept_parameters::ParameterRange, fixed_parameters::ParameterList; random_warmup=false, threading=false, sorting="nearest")   
     # make sure the variables are in our namespace to make them accessible later
     declare_variable.(string.(cat(prob.parameters, prob.variables, dims=1)))
 
@@ -95,7 +93,7 @@ function get_steady_states(prob::Problem, swept_parameters::ParameterRange, fixe
     unique_fixed = filter_duplicate_parameters(swept_parameters, fixed_parameters)
     input_array = _prepare_input_params(prob, swept_parameters, unique_fixed)
     # feed the array into HomotopyContinuation, get back an similar array of solutions
-    raw = _get_raw_solution(prob, input_array, sweep=swept_parameters, random_warmup=random_warmup, threading=threading,show_progress=show_progress)
+    raw = _get_raw_solution(prob, input_array, sweep=swept_parameters, random_warmup=random_warmup, threading=threading)
 
     # 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)
@@ -228,20 +226,22 @@ end
 
 
 "Uses HomotopyContinuation to solve `problem` at specified `parameter_values`."
-function _get_raw_solution(problem::Problem, parameter_values::Array{ParameterVector}; sweep=[],random_warmup=false, threading=false,show_progress= true)
+function _get_raw_solution(problem::Problem, parameter_values::Array{ParameterVector}; sweep=[],random_warmup=false, threading=false)
     # HomotopyContinuation accepts 1D arrays of parameter sets
     params_1D = reshape(parameter_values, :, 1)
 
     if random_warmup
         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, show_progress=show_progress)
-        result_full = HomotopyContinuation.solve(problem.system, HomotopyContinuation.solutions(warmup_solution), start_parameters=warmup_parameters, target_parameters=parameter_values, threading=threading,show_progress= show_progress)
+        warmup_solution = HomotopyContinuation.solve(problem.system,  target_parameters=warmup_parameters, threading=threading)
+        result_full = HomotopyContinuation.solve(problem.system, HomotopyContinuation.solutions(warmup_solution), start_parameters=warmup_parameters, target_parameters=parameter_values, threading=threading)
     else
         result_full = Array{Vector{Any}, 1}(undef, length(parameter_values))
+        bar = Progress(length(parameter_values), 1, "Solving via total degree homotopy ...", 50)
         for i in 1:length(parameter_values) # do NOT thread this
             p = parameter_values[i]
-            result_full[i] = [HomotopyContinuation.solve(problem.system, start_system=:total_degree, target_parameters=p, threading=threading, show_progress=show_progress), p]
+            next!(bar)
+            result_full[i] = [HomotopyContinuation.solve(problem.system, start_system=:total_degree, target_parameters=p, threading=threading, show_progress=false), p]
         end
     end
 

src/sorting.jl

@@ -113,7 +113,10 @@ Go through a vector of solution and sort each according to Euclidean norm.
 function sort_1D(solns::Vector{Vector{SteadyState}})
     sorted_solns = similar(solns) # preallocate
     sorted_solns[1] = sort(solns[1], by= x->abs.(imag(x))) # prefer real solution at first position
+
+    bar = Progress(length(solns), 1, "Ordering solutions into branches ...")
     for (i, soln) in enumerate(solns[1:end-1])
+        next!(bar)
         matched_indices = align_pair(sorted_solns[i], solns[i+1]) # pairs of matching indices
         next_indices = getindex.(matched_indices, 2) # indices of the next solution
         sorted_solns[i+1] = (solns[i+1])[next_indices]
@@ -124,16 +127,16 @@ end
 function hilbert_indices(solns::Matrix{Vector{Vector{ComplexF64}}})
     """Get mapping between 2D indexes (parameter space) and a 1D Hilbert curve"""
     Lx,Ly = size(solns) 
-    mapping = [] #compute mapping between Hilbert indices and 2Ds
-    for j in 1:Ly #length of parameter sweep 1
-        for i in 1:Lx #length of parameter sweep 2
+    mapping = [] # compute mapping between Hilbert indices and 2Ds
+    for j in 1:Ly # length of parameter sweep 1
+        for i in 1:Lx # length of parameter sweep 2
             X = [i, j]
             h = encode_hilbert(Simple2D(Int), X)
             X .= 0
             push!(mapping,(h=>decode_hilbert!(Simple2D(Int), X, h)))
         end 
     end
-    idx_pairs = [el[2] for el in sort(mapping)] #sort along the Hilbert curve. Now we can iterate over these indexes
+    idx_pairs = [el[2] for el in sort(mapping)] # sort along the Hilbert curve. Now we can iterate over these indexes
 end
 
 function naive_indices(solns::Matrix{Vector{Vector{ComplexF64}}})
@@ -164,16 +167,19 @@ end
 
 function sort_2D(solns::Matrix{Vector{Vector{ComplexF64}}}; sorting="nearest") 
     """match each 2D solution with all its surrounding neighbors, including the diagonal ones"""
-     #determine a trajectory in 2D space where nodes will be visited
-    if sorting=="hilbert" #propagating matching of solutions along a hilbert_curve in 2D
+     # determine a trajectory in 2D space where nodes will be visited
+    if sorting=="hilbert" # propagating matching of solutions along a hilbert_curve in 2D
         idx_pairs = hilbert_indices(solns)
-    elseif sorting=="nearest" #propagate matching of solutions along rows
+    elseif sorting=="nearest" # propagate matching of solutions along rows
         idx_pairs = naive_indices(solns)
     end
 
-    sorted_solns = Inf.*copy(solns)  #infinite solutions are ignored by the align_pair function. This trick allows a consistent ordering "propagation"
+    sorted_solns = Inf.*copy(solns)  # infinite solutions are ignored by the align_pair function. This trick allows a consistent ordering "propagation"
     sorted_solns[1,1] = sort(solns[1,1], by= x->abs.(imag(x))) # prefer real solution at first position
+
+    bar = Progress(length(idx_pairs), 1, "Ordering solutions into branches ...")
     for i in 1:length(idx_pairs)-1 
+        next!(bar)
         neighbors =  get_nn_2D(idx_pairs[i+1],size(solns,1),size(solns,2))
         reference = [sorted_solns[ind...] for ind in neighbors]
         matched_indices = align_pair(reference, solns[idx_pairs[i+1]...]) # pairs of matching indices