Merge branch 'master' into krylov

.github/workflows/Invalidations.yml

@@ -0,0 +1,40 @@
+name: Invalidations
+
+on:
+  pull_request:
+
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: always.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  evaluate:
+    # Only run on PRs to the default branch.
+    # In the PR trigger above branches can be specified only explicitly whereas this check should work for master, main, or any other default branch
+    if: github.base_ref == github.event.repository.default_branch
+    runs-on: ubuntu-latest
+    steps:
+    - uses: julia-actions/setup-julia@latest
+      with:
+        version: '1'
+    - uses: actions/checkout@v4
+    - uses: julia-actions/julia-buildpkg@latest
+    - uses: julia-actions/julia-invalidations@v1
+      id: invs_pr
+
+    - uses: actions/checkout@v4
+      with:
+        ref: ${{ github.event.repository.default_branch }}
+    - uses: julia-actions/julia-buildpkg@latest
+    - uses: julia-actions/julia-invalidations
+      id: invs_default
+
+    - name: Report invalidation counts
+      run: |
+        echo "Invalidations on default branch: ${{ steps.invs_default.outputs.total }} (${{ steps.invs_default.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY
+        echo "This branch: ${{ steps.invs_pr.outputs.total }} (${{ steps.invs_pr.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY
+    - name: Check if the PR does increase number of invalidations
+      if: steps.invs_pr.outputs.total > steps.invs_default.outputs.total
+      run: exit 1

Project.toml

@@ -1,14 +1,15 @@
 name = "HarmonicBalance"
 uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 authors = ["Jan Kosata <[email protected]>", "Javier del Pino <[email protected]>", "Orjan Ameye <[email protected]>"]
-version = "0.10.11"
+version = "0.11.1"
 
 [deps]
 BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+EndpointRanges = "340492b5-2a47-5f55-813d-aca7ddf97656"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
@@ -42,6 +43,7 @@ DelimitedFiles = "1.9"
 Distances = "0.10.11"
 DocStringExtensions = "0.9.3"
 Documenter = "1.4"
+EndpointRanges = "0.2.2"
 ExplicitImports = "1.6"
 FFTW = "1.8"
 HomotopyContinuation = "2.9"
@@ -51,15 +53,15 @@ Latexify = "0.16"
 ModelingToolkit = "9.34"
 NonlinearSolve = "3.14"
 OrderedCollections = "1.6"
-OrdinaryDiffEqTsit5 = "1.1"
 OrdinaryDiffEqRosenbrock = "1.1"
+OrdinaryDiffEqTsit5 = "1.1"
 Peaks = "0.5"
 Plots = "1.39"
 PrecompileTools = "1.2"
 ProgressMeter = "1.7.2"
 SteadyStateDiffEq = "2.3.2"
-SymbolicUtils = "3.7"
-Symbolics = "~6.14"
+SymbolicUtils = "3.5"
+Symbolics = "6.4"
 julia = "1.10.0"
 
 [extras]
@@ -69,8 +71,8 @@ ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

benchmark/parametron.jl

@@ -1,4 +1,5 @@
 using HarmonicBalance
+using BenchmarkTools
 
 using Random
 const SEED = 0xd8e5d8df
@@ -18,8 +19,13 @@ dEOM = DifferentialEquation(natural_equation + forces, x)
 add_harmonic!(dEOM, x, ω)
 
 @time harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t);
-@time prob = HarmonicBalance.Problem(harmonic_eq)
 
 fixed = (Ω => 1.0, γ => 1e-2, λ => 5e-2, F => 0, α => 1.0, η => 0.3, θ => 0, ψ => 0)
-varied = ω => range(0.9, 1.1, 20)
-@time res = get_steady_states(prob, varied, fixed; show_progress=false)
+varied = ω => range(0.9, 1.1, 100)
+
+@btime res = get_steady_states(harmonic_eq, WarmUp(), varied, fixed; show_progress=false)
+@btime res = get_steady_states(harmonic_eq, WarmUp(;compile=true), varied, fixed; show_progress=false)
+@btime res = get_steady_states(harmonic_eq, Polyhedral(;only_non_zero=true), varied, fixed; show_progress=false)
+@btime res = get_steady_states(harmonic_eq, Polyhedral(;only_non_zero=false), varied, fixed; show_progress=false)
+@btime res = get_steady_states(harmonic_eq, TotalDegree(;compile=true), varied, fixed; show_progress=false)
+@btime res = get_steady_states(harmonic_eq, TotalDegree(), varied, fixed; show_progress=false)

docs/make.jl

@@ -52,6 +52,6 @@ if CI
         devbranch="master",
         target="build",
         branch="gh-pages",
-        push_preview=false,
+        push_preview=true,
     )
 end

docs/src/.vitepress/config.mts

@@ -64,6 +64,7 @@ export default defineConfig({
         text: 'Manual', items: [
           { text: 'Entering equations of motion', link: '/manual/entering_eom.md' },
           { text: 'Computing effective system', link: '/manual/extracting_harmonics' },
+          { text: 'Computing steady states', link: '/manual/methods' },
           { text: 'Krylov-Bogoliubov', link: '/manual/Krylov-Bogoliubov_method' },
           { text: 'Time evolution', link: '/manual/time_dependent' },
           { text: 'Linear response', link: '/manual/linear_response' },
@@ -100,6 +101,7 @@ export default defineConfig({
       "/manual/": [
         { text: 'Entering equations of motion', link: '/manual/entering_eom.md' },
         { text: 'Computing effective system', link: '/manual/extracting_harmonics' },
+        { text: 'Computing steady states', link: '/manual/methods' },
         { text: 'Krylov-Bogoliubov', link: '/manual/Krylov-Bogoliubov_method' },
         { text: 'Time evolution', link: '/manual/time_dependent' },
         { text: 'Linear response', link: '/manual/linear_response' },

docs/src/.vitepress/theme/style.css

@@ -28,7 +28,7 @@ https://github.com/vuejs/vitepress/blob/main/src/client/theme-default/styles/var
     Cantarell, "Fira Sans", "Droid Sans", "Helvetica Neue", sans-serif;
 
   /* Code Snippet font */
-  --vp-font-family-mono: JuliaMono-Regular, monospace;
+  --vp-font-family-mono: "Julia Mono", Menlo, Monaco, Consolas, "Courier New", monospace;
 
 }
 

docs/src/examples/parametric_via_three_wave_mixing.md

@@ -1,5 +1,5 @@
 ```@meta
-EditURL = "parametric_via_three_wave_mixing.jl"
+EditURL = "../../../examples/parametric_via_three_wave_mixing.jl"
 ```
 
 # Parametric Pumping via Three-Wave Mixing
@@ -33,7 +33,7 @@ If we both have quadratic and cubic nonlineariy, we observe the normal duffing o
 varied = (ω => range(0.99, 1.1, 200)) # range of parameter values
 fixed = (α => 1.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025) # fixed parameters
 
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 plot(result; y="u1^2+v1^2")
 ````
 
@@ -43,7 +43,7 @@ If we set the cubic nonlinearity to zero, we recover the driven damped harmonic
 varied = (ω => range(0.99, 1.1, 100))
 fixed = (α => 0.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025)
 
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 plot(result; y="u1^2+v1^2")
 ````
 
@@ -65,17 +65,16 @@ harmonic_eq2 = get_krylov_equations(diff_eq; order=2)
 varied = (ω => range(0.4, 1.1, 500))
 fixed = (α => 1.0, β => 2.0, ω0 => 1.0, γ => 0.001, F => 0.005)
 
-result = get_steady_states(harmonic_eq2, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq2, varied, fixed)
 plot(result; y="v1")
 ````
 
 ````@example parametric_via_three_wave_mixing
 varied = (ω => range(0.4, 0.6, 100), F => range(1e-6, 0.01, 50))
 fixed = (α => 1.0, β => 2.0, ω0 => 1.0, γ => 0.01)
 
-result = get_steady_states(
-    harmonic_eq2, varied, fixed; threading=true, method=:total_degree
-)
+method = TotalDegree()
+result = get_steady_states(harmonic_eq2, method, varied, fixed)
 plot_phase_diagram(result; class="stable")
 ````
 

docs/src/examples/parametron.md

@@ -1,5 +1,5 @@
 ```@meta
-EditURL = "parametron.jl"
+EditURL = "../../../examples/parametron.jl"
 ```
 
 # [Parametrically driven resonator](@id parametron)
@@ -64,7 +64,7 @@ varied = ω => range(0.9, 1.1, 100)
 result = get_steady_states(harmonic_eq, varied, fixed)
 ````
 
-In `get_steady_states`, the default value for the keyword `method=:random_warmup` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `:total_degree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower. The `threading` keyword enables parallel tracking of homotopy paths, and it's set to `false` simply because we are using a single core computer for now.
+In `get_steady_states`, the default method `WarmUp()` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `TotalDegree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower.
 
 After solving the system, we can save the full output of the simulation and the model (e.g. symbolic expressions for the harmonic equations) into a file
 
@@ -100,6 +100,7 @@ The parametrically driven oscillator boasts a stability diagram called "Arnold's
 To perform a 2D sweep over driving frequency $\omega$ and parametric drive strength $\lambda$, we keep `fixed` from before but include 2 variables in `varied`
 
 ````@example parametron
+fixed = (ω₀ => 1.0, γ => 1e-2, F => 1e-3, α => 1.0, η => 0.3)
 varied = (ω => range(0.8, 1.2, 50), λ => range(0.001, 0.6, 50))
 result_2D = get_steady_states(harmonic_eq, varied, fixed);
 nothing #hide

docs/src/examples/wave_mixing.md

@@ -1,5 +1,5 @@
 ```@meta
-EditURL = "wave_mixing.jl"
+EditURL = "../../../examples/wave_mixing.jl"
 ```
 
 # Three Wave Mixing vs four wave mixing
@@ -39,7 +39,7 @@ response with no response at $2\omega$.
 ````@example wave_mixing
 varied = (ω => range(0.9, 1.2, 200)) # range of parameter values
 fixed = (α => 1.0, β => 0.0, ω0 => 1.0, γ => 0.005, F => 0.0025) # fixed parameters
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)# compute steady states
+result = get_steady_states(harmonic_eq, varied, fixed)# compute steady states
 
 p1 = plot(result; y="√(u1^2+v1^2)", legend=:best)
 p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1))
@@ -57,7 +57,7 @@ We would like to investigate the three-wave mixing of the driven Duffing oscilla
 ````@example wave_mixing
 varied = (ω => range(0.9, 1.2, 200))
 fixed = (α => 0.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025)
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 
 p1 = plot(result; y="√(u1^2+v1^2)", legend=:best)
 p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1))
@@ -75,7 +75,7 @@ We would like to investigate the three-wave mixing of the driven Duffing oscilla
 ````@example wave_mixing
 varied = (ω => range(0.9, 1.2, 200))
 fixed = (α => 1.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025)
-result = get_steady_states(harmonic_eq, varied, fixed; threading=true)
+result = get_steady_states(harmonic_eq, varied, fixed)
 
 p1 = plot(result; y="√(u1^2+v1^2)", legend=:best)
 p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1))

docs/src/manual/methods.md

@@ -0,0 +1,18 @@
+# Methods
+
+We offer several methods for solving the nonlinear algebraic equations that arise from the harmonic balance procedure. Each method has different tradeoffs between speed, robustness, and completeness.
+
+## Total Degree Method
+```@docs
+TotalDegree
+```
+
+## Polyhedral Method
+```@docs
+Polyhedral
+```
+
+## Warm Up Method
+```@docs
+WarmUp
+```

docs/src/manual/plotting.md

@@ -12,7 +12,7 @@ The function `plot` is multiple-dispatched to plot 1D and 2D datasets.
 In 1D, the solutions are colour-coded according to the branches obtained by `sort_solutions`. 
 
 ```@docs
-HarmonicBalance.plot(::Result, varags...)
+HarmonicBalance.plot(::HarmonicBalance.Result, varags...)
 ```
 
 

docs/src/manual/solving_harmonics.md

@@ -8,9 +8,9 @@ having called `get_harmonic_equations`, we need to set all time-derivatives to z
 Once defined, a `Problem` can be solved for a set of input parameters using `get_steady_states` to obtain `Result`.
 
 ```@docs
-Problem
+HarmonicBalance.Problem
 get_steady_states
-Result
+HarmonicBalance.Result
 ```
 
 

docs/src/manual/time_dependent.md

@@ -2,10 +2,10 @@
 
 Generally, solving the ODE of oscillatory systems in time requires numerically tracking the oscillations. This is a computationally expensive process; however, using the harmonic ansatz removes the oscillatory time-dependence. Simulating instead the harmonic variables of a `HarmonicEquation` is vastly more efficient - a steady state of the system appears as a fixed point in multidimensional space rather than an oscillatory function.
 
-The Extention `TimeEvolution` is used to interface `HarmonicEquation` with the solvers contained in `OrdinaryDiffEq.jl`. Time-dependent parameter sweeps are defined using the object `ParameterSweep`. To use the `TimeEvolution` extension, one must first load the `OrdinaryDiffEq.jl` package.
+The Extention `TimeEvolution` is used to interface `HarmonicEquation` with the solvers contained in `OrdinaryDiffEq.jl`. Time-dependent parameter sweeps are defined using the object `AdiabaticSweep`. To use the `TimeEvolution` extension, one must first load the `OrdinaryDiffEq.jl` package.
 ```@docs
 ODEProblem(::HarmonicEquation, ::Any; timespan::Tuple)
-ParameterSweep
+AdiabaticSweep
 ```
 
 ## Plotting

docs/src/tutorials/classification.md

@@ -20,7 +20,7 @@ We performe a 2d sweep in the driving frequency $\omega$ and driving strength $\
 fixed = (ω₀ => 1.0, γ => 0.002, α => 1.0)
 varied = (ω => range(0.99, 1.01, 100), λ => range(1e-6, 0.03, 100))
 
-result_2D = get_steady_states(harmonic_eq, varied, fixed, threading=true)
+result_2D = get_steady_states(harmonic_eq, varied, fixed)
 ```
 By default the steady states of the system are classified by four different catogaries:
 * `physical`: Solutions that are physical, i.e., all variables are purely real.

docs/src/tutorials/limit_cycles.md

@@ -62,8 +62,8 @@ Here we reconstruct the results of [Zambon et al., Phys Rev. A 102, 023526 (2020
 ```@example lc
 using HarmonicBalance
 @variables γ F α ω0 F0 η ω J t x(t) y(t);
-eqs = [d(x,t,2) + γ*d(x,t) + ω0^2*x + α*x^3+ 2*J*ω0*(x-y) - F0*cos(ω*t),
-       d(y,t,2) + γ * d(y,t) + ω0^2 * y + α*y^3 + 2*J*ω0*(y-x) - η*F0*cos(ω*t)]
+eqs = [d(x,t,2) + γ*d(x,t) + ω0^2*x + α*x^3 + 2*J*ω0*(x-y) - F0*cos(ω*t),
+       d(y,t,2) + γ*d(y,t) + ω0^2*y + α*y^3 + 2*J*ω0*(y-x) - η*F0*cos(ω*t)]
 diff_eq = DifferentialEquation(eqs, [x,y])
 ```
 
@@ -82,8 +82,7 @@ fixed = (
     J => 154.1e-6,   # coupling term
     α => 3.867e-7,   # Kerr nonlinearity
     ω => 1.4507941,  # pump frequency, resonant with antisymmetric mode (in paper, ħω0 + J)
-    η => -0.08,      # pumping leaking to site 2  (F2 = ηF1)
-    F0 => 0.002       # pump amplitude (overriden in sweeps)
+    η => -0.08     # pumping leaking to site 2  (F2 = ηF1)
 )
 varied = F0 => range(0.002, 0.03, 50)
 
@@ -105,7 +104,7 @@ using OrdinaryDiffEqTsit5
 initial_state = result[1][1]
 
 T = 2e6
-sweep = ParameterSweep(F0 => (0.002, 0.011), (0,T))
+sweep = AdiabaticSweep(F0 => (0.002, 0.011), (0,T))
 
 # start from initial_state, use sweep, total time is 2*T
 time_problem = ODEProblem(harmonic_eq, initial_state, sweep=sweep, timespan=(0,2*T))