From a4bd550def064fba3281f1173167794ad6b86ed7 Mon Sep 17 00:00:00 2001
From: Jan Kosata <58295890+jkosata@users.noreply.github.com>
Date: Thu, 21 Apr 2022 12:14:16 +0200
Subject: [PATCH] Progress bars (#20)
* progress bars added
* Symbolics updated
---
Project.toml | 6 +++--
src/HarmonicBalance.jl | 1 +
src/modules/LinearResponse.jl | 1 +
.../LinearResponse/jacobian_spectrum.jl | 3 +++
src/modules/LinearResponse/linear_response.jl | 3 +++
src/solve_homotopy.jl | 18 +++++++--------
src/sorting.jl | 22 ++++++++++++-------
7 files changed, 35 insertions(+), 19 deletions(-)
diff --git a/Project.toml b/Project.toml
index 6ab75cfa2..5986b7d81 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
name = "HarmonicBalance"
uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
authors = ["Jan Kosata <kosataj@phys.ethz.ch>", "Javier del Pino <jdelpino@phys.ethz.ch>"]
-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]
diff --git a/src/HarmonicBalance.jl b/src/HarmonicBalance.jl
index 9cb5aaab0..c82a18bc3 100644
--- a/src/HarmonicBalance.jl
+++ b/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
diff --git a/src/modules/LinearResponse.jl b/src/modules/LinearResponse.jl
index ee982082a..5c49779b6 100644
--- a/src/modules/LinearResponse.jl
+++ b/src/modules/LinearResponse.jl
@@ -7,6 +7,7 @@ module LinearResponse
using ..HarmonicBalance
using ..HC_wrapper
using PyPlot
+ using ProgressMeter
using DocStringExtensions
using Latexify
diff --git a/src/modules/LinearResponse/jacobian_spectrum.jl b/src/modules/LinearResponse/jacobian_spectrum.jl
index 4cc7948b8..0e5942068 100644
--- a/src/modules/LinearResponse/jacobian_spectrum.jl
+++ b/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
diff --git a/src/modules/LinearResponse/linear_response.jl b/src/modules/LinearResponse/linear_response.jl
index add361858..6e839c091 100644
--- a/src/modules/LinearResponse/linear_response.jl
+++ b/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])
diff --git a/src/solve_homotopy.jl b/src/solve_homotopy.jl
index b57f81177..1ff944c11 100644
--- a/src/solve_homotopy.jl
+++ b/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
diff --git a/src/sorting.jl b/src/sorting.jl
index 571b3f780..9a72b3e35 100644
--- a/src/sorting.jl
+++ b/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