QuantumEngineeredSystems · oameye · Jun 28, 2025 · Jun 28, 2025 · Jun 28, 2025 · Jun 28, 2025
diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
@@ -16,6 +16,7 @@ permissions:
   pages: write
   id-token: write
   statuses: write
+  pull-requests: read  # Required when using `push_preview=true`
 
 # Allow only one concurrent deployment, skipping runs queued between the run in-progress and latest queued.
 # However, do NOT cancel in-progress runs as we want to allow these production deployments to complete.

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -13,7 +13,6 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

diff --git a/docs/make.jl b/docs/make.jl
@@ -16,6 +16,7 @@ using SteadyStateDiffEq
 TimeEvolution = Base.get_extension(HarmonicSteadyState, :TimeEvolution)
 ModelingToolkitExt = Base.get_extension(HarmonicBalance, :ModelingToolkitExt)
 SteadyStateDiffEqExt = Base.get_extension(HarmonicSteadyState, :SteadyStateDiffEqExt)
+HarmonicBalanceExt = Base.get_extension(HarmonicSteadyState, :HarmonicBalanceExt)
 
 bib = CitationBibliography(
     joinpath(@__DIR__, "src", "refs.bib");
@@ -34,6 +35,17 @@ include("make_md_examples.jl")
 
 include("pages.jl")
 
+# Create remotes for Documenter
+using Pkg
+remotes = Dict()
+for ext in [QuestBase, HarmonicSteadyState]
+    str = string(ext)
+    status = sprint(io -> Pkg.status(str; io=io))
+    version = match(r"(v[0-9].[0-9]+.[0-9]+)", status)[1]
+    gh_moi = Documenter.Remotes.GitHub("QuantumEngineeredSystems", str * ".jl")
+    remotes[pkgdir(ext)] = (gh_moi, version)
+end
+
 makedocs(;
     sitename="HarmonicBalance.jl",
     authors="Quest group",
@@ -46,6 +58,7 @@ makedocs(;
         SteadyStateDiffEqExt,
         HarmonicSteadyState.LinearResponse,
         PlotsExt,
+        HarmonicBalanceExt,
     ],
     format=DocumenterVitepress.MarkdownVitepress(;
         repo="github.com/QuantumEngineeredSystems/HarmonicBalance.jl",
@@ -54,6 +67,7 @@ makedocs(;
     ),
     checkdocs=:exports,
     pages=pages,
+    remotes=remotes,
     source="src",
     build="build",
     draft=(!CI),

diff --git a/docs/src/examples/cumulants_KPO.md b/docs/src/examples/cumulants_KPO.md
@@ -36,7 +36,9 @@ To compute the steady states of the KPO, we use `HarmonicBalance`. We define the
 ````@example cumulants_KPO
 fixed = (U => 0.001, κ => 0.002)
 varied = (Δ => range(-0.03, 0.03, 100), G => range(1e-5, 0.02, 100))
-problem_c1 = HarmonicSteadyState.HomotopyContinuationProblem(eqs_completed_RWA, param, varied, fixed)
+problem_c1 = HarmonicSteadyState.HomotopyContinuationProblem(
+    eqs_completed_RWA, param, varied, fixed
+)
 ````
 
 This gives us the phase

diff --git a/docs/src/examples/harmonic_oscillator_KB_vs_HB.md b/docs/src/examples/harmonic_oscillator_KB_vs_HB.md
@@ -28,9 +28,10 @@ harmonic_eq = rearrange_standard(harmonic_eq)
 ````@example harmonic_oscillator_KB_vs_HB
 varied = (ω => range(0.1, 1.9, 200)) # range of parameter values
 fixed = (ω0 => 1.0, γ => 0.05, F => 0.1) # fixed parameters
-result_krylov1 = get_steady_states(krylov_eq1, varied, fixed)
-result_krylov2 = get_steady_states(krylov_eq2, varied, fixed)
-result_harmonic = get_steady_states(harmonic_eq, varied, fixed);
+show_progress = false # show progress bar
+result_krylov1 = get_steady_states(krylov_eq1, varied, fixed; show_progress)
+result_krylov2 = get_steady_states(krylov_eq2, varied, fixed; show_progress)
+result_harmonic = get_steady_states(harmonic_eq, varied, fixed; show_progress);
 nothing #hide
 ````
 
@@ -80,9 +81,13 @@ plot(
     ),
     plot_linear_response(
         result_harmonic, x, 1; Ω_range=range(0.1, 1.9, 200), title="Harmonic"
+    ),
+    plot_linear_response(
+        result_harmonic, x, 1; Ω_range=range(0.1, 1.9, 200), title="Exact", order=2
     );
-    layout=(3, 1),
+    layout=(4, 1),
     clims=(0, 250),
+    size=(800, 600),
 )
 ````
 

diff --git a/docs/src/examples/state_dependent_perturbation.md b/docs/src/examples/state_dependent_perturbation.md
@@ -11,6 +11,7 @@ Let us load the following packages into our environment:
 ````@example state_dependent_perturbation
 using HarmonicBalance, Plots;
 HB = HarmonicBalance;
+HSS = HarmonicSteadyState;
 crange = (0, 9);
 nothing #hide
 ````
@@ -19,15 +20,15 @@ Later, we will need to classify the solutions of our perturbation. To make this
 
 ````@example state_dependent_perturbation
 function my_classify_default!(result)
-    my_classify_solutions!(result, HB.is_physical, "physical")
-    my_classify_solutions!(result, HB.is_stable, "stable")
-    return HB.order_branches!(result, ["physical", "stable"]) # shuffle the branches to have relevant
+    my_classify_solutions!(result, HSS.is_physical, "physical")
+    my_classify_solutions!(result, HSS.is_stable, "stable")
+    return HSS.order_branches!(result, ["physical", "stable"]) # shuffle the branches to have relevant
 end
-function my_classify_solutions!(res::HB.Result, f::Function, name::String)
+function my_classify_solutions!(res::HSS.Result, f::Function, name::String)
     values = my_classify_solutions(res, f)
     return res.classes[name] = values
 end
-function my_classify_solutions(res::HB.Result, f::Function)
+function my_classify_solutions(res::HSS.Result, f::Function)
     values = similar(res.solutions, BitVector)
     for (idx, soln) in enumerate(res.solutions)
         values[idx] = [
@@ -72,8 +73,8 @@ We sweep over the system where we both increase the drive frequency $\omega$ and
 
 ````@example state_dependent_perturbation
 res = 80
-fixed = HB.OrderedDict(ω0 => 1.0, α => 1.0, J => 0.005, γ => 0.005)
-varied = HB.OrderedDict((ω => range(0.99, 1.01, res), λ => range(1e-6, 0.03, res)))
+fixed = HSS.OrderedDict(ω0 => 1.0, α => 1.0, J => 0.005, γ => 0.005)
+varied = HSS.OrderedDict((ω => range(0.99, 1.01, res), λ => range(1e-6, 0.03, res)))
 method = TotalDegree()
 result_ωλ = get_steady_states(harmonic_normal, method, varied, fixed; show_progress=false);
 plot_phase_diagram(result_ωλ; class="stable")
@@ -206,7 +207,7 @@ plot_phase_diagram(result_ωλ_antisym; class=["stable", "not_zero"])
 We filter the non-zero amplitude solution and store it in a matrix $A$:
 
 ````@example state_dependent_perturbation
-branch_mat = findfirst.(HB._get_mask(result_ωλ_antisym, ["stable", "not_zero"], []))
+branch_mat = findfirst.(HSS._get_mask(result_ωλ_antisym, ["stable", "not_zero"], []))
 A = map(CartesianIndices(result_ωλ_antisym.solutions)) do idx
     branch = branch_mat[idx]
     if isnothing(branch)
@@ -242,7 +243,7 @@ permutation = first.(
 param_ranges = collect(values(varied))
 input_array = collect(Iterators.product(param_ranges..., values(fixed)...))
 input_array = getindex.(input_array, [permutation])
-input_array = HB.tuple_to_vector.(input_array)
+input_array = HSS.tuple_to_vector.(input_array)
 input_array = map(idx -> push!(input_array[idx], A[idx]...), CartesianIndices(input_array));
 nothing #hide
 ````
@@ -251,8 +252,8 @@ Solving for the steady states of the dressed symmetric mode:
 
 ````@example state_dependent_perturbation
 function solve_perturbed_system(prob, input)
-    result_full = HB.ProgressMeter.@showprogress map(input_array) do input
-        HB.HomotopyContinuation.solve(
+    result_full = HSS.ProgressMeter.@showprogress map(input_array) do input
+        HSS.HomotopyContinuation.solve(
             prob.system;
             start_system=:total_degree,
             target_parameters=input,
@@ -261,22 +262,22 @@ function solve_perturbed_system(prob, input)
         )
     end
 
-    rounded_solutions = HB.HomotopyContinuation.solutions.(result_full)
-    solutions = HB.pad_solutions(rounded_solutions)
+    rounded_solutions = HSS.HomotopyContinuation.solutions.(result_full)
+    solutions = HSS.pad_solutions(rounded_solutions)
 
     jacobian = HarmonicBalance.get_Jacobian(harmonic_tmp)
     J_variables = cat(prob.variables, collect(keys(varied)), [ua, va]; dims=1)
-    compiled_J = HB.compile_matrix(jacobian, J_variables; rules=fixed)
-    compiled_J = HB.JacobianFunction(HB.solution_type(solutions))(compiled_J)
-    result = HB.Result(
+    compiled_J = HSS.compile_matrix(jacobian, J_variables; rules=fixed)
+    compiled_J = HSS.JacobianFunction(HSS.solution_type(solutions))(compiled_J)
+    result = HSS.Result(
         solutions,
         varied,
         fixed,
         prob,
         Dict(),
         zeros(Int64, size(solutions)...),
         compiled_J,
-        HB.seed(method),
+        HSS.seed(method),
     )
 
     my_classify_default!(result)

diff --git a/docs/src/examples/steady_state_sweep.md b/docs/src/examples/steady_state_sweep.md
@@ -6,7 +6,8 @@ EditURL = "../../../examples/steady_state_sweep.jl"
 
 ````@example steady_state_sweep
 using HarmonicBalance, SteadyStateDiffEq, ModelingToolkit
-using BenchmarkTools, Plots, StaticArrays, OrdinaryDiffEq, LinearAlgebra
+using BenchmarkTools, Plots, StaticArrays, LinearAlgebra
+using OrdinaryDiffEqRosenbrock, OrdinaryDiffEqTsit5
 using HarmonicBalance: OrderedDict
 
 @variables α ω ω0 F γ η t x(t);

diff --git a/docs/src/manual/API.md b/docs/src/manual/API.md
@@ -89,6 +89,12 @@ Private = false
 Order = [:function]
 ```
 
+```@autodocs
+Modules = [Base.get_extension(HarmonicSteadyState, :HarmonicBalanceExt)]
+Private = false
+Order = [:function]
+```
+
 ```@docs
 HarmonicBalance.get_Jacobian
 ```

diff --git a/docs/src/manual/linear_response.md b/docs/src/manual/linear_response.md
@@ -34,8 +34,13 @@ HarmonicSteadyState.LinearResponse.Lorentzian
 
 Setting `order > 1` increases the accuracy of the response spectra. However, unlike for the Jacobian, here we must perform a matrix inversion for each response frequency.  
 
+```@autodocs; canonical=false
+Modules = [Base.get_extension(HarmonicSteadyState, :HarmonicBalanceExt)]
+Private = false
+Order = [:function]
+```
+
 ```@docs; canonical=false
-HarmonicSteadyState.LinearResponse.get_linear_response
 HarmonicSteadyState.LinearResponse.ResponseMatrix
 HarmonicSteadyState.LinearResponse.get_response
 ```

diff --git a/examples/Project.toml b/examples/Project.toml
@@ -8,7 +8,6 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 NonlinearSolveHomotopyContinuation = "2ac3b008-d579-4536-8c91-a1a5998c2f8b"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

diff --git a/examples/steady_state_sweep.jl b/examples/steady_state_sweep.jl
@@ -1,7 +1,8 @@
 # # Steady state sweeps
 
 using HarmonicBalance, SteadyStateDiffEq, ModelingToolkit
-using BenchmarkTools, Plots, StaticArrays, OrdinaryDiffEq, LinearAlgebra
+using BenchmarkTools, Plots, StaticArrays, LinearAlgebra
+using OrdinaryDiffEqRosenbrock, OrdinaryDiffEqTsit5
 using HarmonicBalance: OrderedDict
 
 @variables α ω ω0 F γ η t x(t);