split up package in an ecosystem (#419)

* QuestBase

* make QuestBase package

* new HarmonicBalance module

* update QuestBase

* rm old folders

* Make HarmonicSteadyState depending on HarmonicBalance for now

* try to get test suite working of HarmonicSteadyState

* Make HC wrapper a Module again and declare variables again

* woop woop  test working

* exported QuestBase and HarmonicSteadyState

* format

* refactor: move get_Jacobian from QuestBase

* new QuestBase version

* up compat DocStringExtensions

* up compat SymbolicUtils

* format

* CI       fail-fast: false

* up compat Aqua and JET

* refactor: update comment for clarity on future Symbolics.jacobian replacement

* bump QuestBase

* fix: turn of JET on julia prerelease

* start updating docs

* empty commit

* fix examples

* Revert "fix examples"

This reverts commit 15e72e8.

* Reexport HarmonicSteadyState

* limit to julia 1.11 for now

* fix state_dependent_perturbation

* bump Random dependency to version 1.11.0

Author: oameye

.github/workflows/Documentation.yml

@@ -23,6 +23,9 @@ concurrency:
   group: ${{ github.workflow }}-${{ github.ref }}
   cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
 
+env:
+  JULIA_NUM_THREADS: 10
+
 jobs:
   # Build job
   build:
@@ -31,7 +34,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
-          version: 'lts'
+          version: '1'
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-docdeploy@v1

.github/workflows/Tests.yml

@@ -29,8 +29,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - 'pre'
-          - 'lts'
+          # - 'pre'
+          # - 'lts'
           - '1'
         os:
           - ubuntu-latest

Project.toml

@@ -1,57 +1,52 @@
 name = "HarmonicBalance"
 uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 authors = ["Orjan Ameye <orjan.ameye@hotmail.com>", "Jan Kosata <kosataj@phys.ethz.ch>", "Javier del Pino <jdelpino@phys.ethz.ch>"]
-version = "0.14.4"
+version = "0.15.0"
 
 [deps]
-BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"
-Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
-HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+HarmonicSteadyState = "1158f75c-a779-4b85-8bfb-8fcf6bf02ced"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+QuestBase = "7e80f742-43d6-403d-a9ea-981410111d43"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [weakdeps]
-LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-QuantumCumulants = "35bcea6d-e19f-57db-af74-8011de6c7255"
-SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 
 [extensions]
 ModelingToolkitExt = "ModelingToolkit"
-PlotsExt = "Plots"
-QuantumCumulantsExt = "QuantumCumulants"
-SteadyStateDiffEqExt = "SteadyStateDiffEq"
-TimeEvolution = "OrdinaryDiffEqTsit5"
 
 [compat]
-BijectiveHilbert = "0.3.0, 0.5"
-Distances = "0.10.11"
-DocStringExtensions = "0.9.3"
-FunctionWrappers = "1.1.3"
-HomotopyContinuation = "2.12"
-LinearSolve = "2.39.0, 3"
-ModelingToolkit = "9.70.0"
-NaNMath = "1.1"
-OrderedCollections = "1.6"
-OrdinaryDiffEqTsit5 = "1.1"
-Plots = "1.40.9"
+Aqua = "0.8.11"
+DocStringExtensions = "0.9.4"
+Documenter = "1.4"
+ExplicitImports = "1.6"
+HarmonicSteadyState = "0.2.0"
+JET = "0.9.18"
+ModelingToolkit = "9.60"
 PrecompileTools = "1.2"
-ProgressMeter = "1.7.2"
-QuantumCumulants = "0.3.7"
-SciMLBase = "2.83.1"
-SteadyStateDiffEq = "2.5.0"
-SymbolicUtils = "3.25.1"
-Symbolics = "~6.34"
-julia = "1.10.0"
+QuestBase = "0.3.0"
+Random = "1.11.0"
+Reexport = "1.2.2"
+SymbolicUtils = "3.25"
+Symbolics = "~6.34.0"
+Test = "1.11.0"
+julia = "1.11.0"
+
+[extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+QuestBase = "7e80f742-43d6-403d-a9ea-981410111d43"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Aqua", "JET", "SymbolicUtils", "Symbolics", "ExplicitImports", "Test", "Random", "ModelingToolkit", "Documenter", "QuestBase"]

benchmark/3parametron.jl

@@ -23,7 +23,7 @@ add_harmonic!(dEOM, x, ω)
 fixed = (Ω => 1.0, γ => 1e-2, λ => 5e-2, F => 0, α => 1.0, η => 0.3, θ => 0, ψ => 0)
 varied = ω => range(0.9, 1.1, 100)
 
-prob = HarmonicBalance.Problem(harmonic_eq)
+prob = HarmonicSteadyState.HomotopyContinuationProblem(harmonic_eq)
 @btime res = get_steady_states(prob, WarmUp(), varied, fixed; show_progress=false) # 380.126 ms (925069 allocations: 50.84 MiB)
 @btime res = get_steady_states(
     prob, WarmUp(; compile=true), varied, fixed; show_progress=false

docs/Project.toml

@@ -6,6 +6,7 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 HarmonicBalance = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
+HarmonicSteadyState = "1158f75c-a779-4b85-8bfb-8fcf6bf02ced"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
@@ -17,6 +18,7 @@ OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 QuantumCumulants = "35bcea6d-e19f-57db-af74-8011de6c7255"
+QuestBase = "7e80f742-43d6-403d-a9ea-981410111d43"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"

docs/make.jl

@@ -1,6 +1,8 @@
 CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
 
 using HarmonicBalance
+using HarmonicSteadyState
+using QuestBase
 
 using Documenter
 using DocumenterVitepress
@@ -11,17 +13,17 @@ using ModelingToolkit
 using OrdinaryDiffEqTsit5
 using SteadyStateDiffEq
 
-TimeEvolution = Base.get_extension(HarmonicBalance, :TimeEvolution)
+TimeEvolution = Base.get_extension(HarmonicSteadyState, :TimeEvolution)
 ModelingToolkitExt = Base.get_extension(HarmonicBalance, :ModelingToolkitExt)
-SteadyStateDiffEqExt = Base.get_extension(HarmonicBalance, :SteadyStateDiffEqExt)
+SteadyStateDiffEqExt = Base.get_extension(HarmonicSteadyState, :SteadyStateDiffEqExt)
 
 bib = CitationBibliography(
     joinpath(@__DIR__, "src", "refs.bib");
     style=:numeric,  # default
 )
 
 using Plots
-PlotsExt = Base.get_extension(HarmonicBalance, :PlotsExt)
+PlotsExt = Base.get_extension(HarmonicSteadyState, :PlotsExt)
 default(; fmt=:png)
 # Gotta set this environment variable when using the GR run-time on CI machines.
 # This happens as examples will use Plots.jl to make plots and movies.
@@ -37,10 +39,12 @@ makedocs(;
     authors="Quest group",
     modules=[
         HarmonicBalance,
+        QuestBase,
+        HarmonicSteadyState,
         TimeEvolution,
         ModelingToolkitExt,
         SteadyStateDiffEqExt,
-        HarmonicBalance.LinearResponse,
+        HarmonicSteadyState.LinearResponse,
         PlotsExt,
     ],
     format=DocumenterVitepress.MarkdownVitepress(;
@@ -58,7 +62,7 @@ makedocs(;
     else
         [:linkcheck, :cross_references, :missing_docs, :docs_block]
     end,
-    doctest=false,  # We test it in the CI, no need to run it here
+    doctest=false, # We test it in the CI, no need to run it here
     plugins=[bib],
 )
 

docs/src/background/stability_response.md

@@ -128,6 +128,6 @@ The above solution applies to every eigenvalue $\lambda$ of the Jacobian. It is
 
 _Sidenote:_ As $J$ a real matrix, there is an eigenvalue $\lambda_r^*$ for each $\lambda_r$. The maximum number of peaks in the linear response is thus equal to the dimensionality of $\mathbf{u}(T)$.
 
-The linear response of the system in the state $\mathbf{u}_0$ is thus fully specified by the complex eigenvalues and eigenvectors of $J(\mathbf{u}_0)$. In HarmonicBalance.jl, the module [LinearResponse](@ref linresp_man) creates a set of plottable [`Lorentzian`](@ref HarmonicBalance.LinearResponse.Lorentzian) objects to represent this.
+The linear response of the system in the state $\mathbf{u}_0$ is thus fully specified by the complex eigenvalues and eigenvectors of $J(\mathbf{u}_0)$. In HarmonicBalance.jl, the module [LinearResponse](@ref linresp_man) creates a set of plottable [`Lorentzian`](@ref HarmonicSteadyState.LinearResponse.Lorentzian) objects to represent this.
 
 [Check out this example](@ref linresp_ex) of the linear response module of HarmonicBalance.jl

docs/src/examples/ab_initio_noise.md

@@ -2,7 +2,7 @@
 EditURL = "../../../examples/ab_initio_noise.jl"
 ```
 
-Ab Initio Noise sidebands and spectra
+# Ab Initio Noise sidebands and spectra
 
 This example demonstrates how to compute the spectra obtained from probing the system with noisy probe.
 
@@ -94,9 +94,7 @@ We will perform parameter sweep to generate noise spectra across the driving fre
 
 ````@example ab_initio_noise
 setter! = setp(sdesystem, ω)
-prob_func(prob, i, repeat) = (prob′ = remake(prob);
-setter!(prob′, ωrange[i]);
-prob′)
+prob_func(prob, i, repeat) = (prob′=remake(prob); setter!(prob′, ωrange[i]); prob′)
 output_func(sol, i) = (outputpsd(sol), false)
 prob_ensemble = EnsembleProblem(sdeproblem; prob_func=prob_func, output_func=output_func)
 sol_ensemble = solve(

docs/src/examples/cumulants_KPO.md

@@ -36,7 +36,7 @@ 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 = HarmonicBalance.Problem(eqs_completed_RWA, param, varied, fixed)
+problem_c1 = HarmonicSteadyState.HomotopyContinuationProblem(eqs_completed_RWA, param, varied, fixed)
 ````
 
 This gives us the phase
@@ -54,7 +54,7 @@ The next order cumulant can be computed by setting the `order` keyword.
 ops = [a]
 eqs_RWA = meanfield(ops, H_RWA, [a]; rates=[κ], order=2)
 eqs_c2 = complete(eqs_RWA)
-problem_c2 = HarmonicBalance.Problem(eqs_c2, param, varied, fixed)
+problem_c2 = HarmonicSteadyState.HomotopyContinuationProblem(eqs_c2, param, varied, fixed)
 ````
 
 Which gives us the phase diagram
@@ -69,7 +69,7 @@ However, the phase diagram seems to wrong. Indeed, plotting, the photon number `
 ````@example cumulants_KPO
 fixed = (U => 0.001, κ => 0.002, G => 0.01)
 varied = (Δ => range(-0.03, 0.03, 200))
-problem_c2 = HarmonicBalance.Problem(eqs_c2, param, varied, fixed)
+problem_c2 = HarmonicSteadyState.HomotopyContinuationProblem(eqs_c2, param, varied, fixed)
 result = get_steady_states(problem_c2, TotalDegree())
 plot(result; y="a⁺aᵣ", class="stable")
 ````
@@ -92,7 +92,7 @@ eqs_c3 = complete(eqs_RWA)
 
 fixed = (U => 0.001, κ => 0.002, G => 0.01)
 varied = (Δ => range(-0.03, 0.03, 50))
-problem_c3 = HarmonicBalance.Problem(eqs_c3, param, varied, fixed)
+problem_c3 = HarmonicSteadyState.HomotopyContinuationProblem(eqs_c3, param, varied, fixed)
 result = get_steady_states(problem_c3, TotalDegree())
 classify_solutions!(result, "a⁺aᵣ < 0", "neg photon number");
 plot(result; y="a⁺aᵣ", class="stable")

docs/src/examples/state_dependent_perturbation.md

@@ -228,19 +228,16 @@ harmonic_tmp.equations = HB.Symbolics.substitute(
     HB.rearrange_standard(harmonic_normal).equations[1:2], Dict(u2 => ua, v2 => va)
 )
 harmonic_tmp.parameters = push!(harmonic_tmp.parameters, ua, va)
-prob = HarmonicBalance.Problem(harmonic_tmp, varied, fixed)
+prob = HarmonicSteadyState.HomotopyContinuationProblem(harmonic_tmp, varied, fixed)
 ````
 
 We will sweep over the $\omega-\lambda$ plane and substitute the non-zero amplitude solution of the antisymmetric mode into the coupled equations of thesymmetric mode.
 
 ````@example state_dependent_perturbation
 all_keys = cat(collect(keys(varied)), collect(keys(fixed)); dims=1)
-permutation =
-    first.(
-        filter(
-            !isempty, [findall(x -> isequal(x, par), all_keys) for par in prob.parameters]
-        )
-    )
+permutation = first.(
+    filter(!isempty, [findall(x -> isequal(x, par), all_keys) for par in prob.parameters])
+)
 
 param_ranges = collect(values(varied))
 input_array = collect(Iterators.product(param_ranges..., values(fixed)...))