build: migrate to Symbolics 7 / SciMLBase 3 / MTKBase stack

.github/codecov.yml

@@ -0,0 +1,10 @@
+coverage:
+  status:
+    project:
+      default:
+        informational: true
+    patch:
+      default:
+        informational: true
+
+comment: false

.github/workflows/Documentation.yml

@@ -32,7 +32,7 @@ jobs:
       - uses: actions/checkout@v6
       - uses: julia-actions/setup-julia@v2
         with:
-          version: 'lts'
+          version: '1.12'
       - uses: julia-actions/cache@v3
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-docdeploy@v1

.github/workflows/Tests.yml

@@ -29,8 +29,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - 'lts'
-          - '1'
+          - '1.12'
         os:
           - ubuntu-latest
         arch:
@@ -60,4 +59,4 @@ jobs:
           files: lcov.info
           name: codecov-umbrella
           fail_ci_if_error: false
-        if: ${{ matrix.version == 'lts' }}
+        if: ${{ matrix.version == '1.12' }}

Project.toml

@@ -1,7 +1,7 @@
 name = "HarmonicSteadyState"
 uuid = "1158f75c-a779-4b85-8bfb-8fcf6bf02ced"
+version = "0.5.0"
 authors = ["Orjan Ameye <orjan.ameye@hotmail.com>", "Jan Kosata <kosataj@phys.ethz.ch>", "Javier del Pino <jdelpino@phys.ethz.ch>"]
-version = "0.4.1"
 
 [deps]
 BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"
@@ -23,7 +23,9 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [weakdeps]
 HarmonicBalance = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
+Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+ModelingToolkitBase = "7771a370-6774-4173-bd38-47e70ca0b839"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -32,47 +34,48 @@ SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 
 [extensions]
 HarmonicBalanceExt = "HarmonicBalance"
-PlotsExt = "Plots"
+PlotsExt = ["Plots", "Latexify"]
 QuantumCumulantsExt = "QuantumCumulants"
 SteadyStateDiffEqExt = "SteadyStateDiffEq"
 TimeEvolution = "OrdinaryDiffEqTsit5"
 
 [compat]
-Peaks = "0.5"
 Aqua = "0.8.11"
-BijectiveHilbert = "0.3, 0.5"
+BijectiveHilbert = "0.3, 0.5, 0.6"
 CheckConcreteStructs = "0.1.0"
 Distances = "0.10.11"
 DocStringExtensions = "0.9.4"
 Documenter = "1.4"
 ExplicitImports = "1.11"
 FunctionWrappers = "1.1.3"
-HarmonicBalance = "0.16"
+HarmonicBalance = "0.17"
 HomotopyContinuation = "2.12"
 JET = "0.9.18, 0.10"
+Latexify = "0.16"
 LinearAlgebra = "1.10"
-LinearSolve = "3.7.2"
-ModelingToolkit = "10"
+LinearSolve = "3.75"
+ModelingToolkitBase = "1"
 NaNMath = "1.1"
-NonlinearSolve = "4.5.1"
-OrderedCollections = "1.8"
-OrdinaryDiffEqRosenbrock = "1.9"
-OrdinaryDiffEqTsit5 = "1.1"
+NonlinearSolve = "4.17"
+OrderedCollections = "1.8, 2"
+OrdinaryDiffEqRosenbrock = "2"
+OrdinaryDiffEqTsit5 = "2"
+Peaks = "0.5"
 Plots = "1.40.9"
 PrecompileTools = "1.2"
 Printf = "1.10"
 ProgressMeter = "1.7.2"
-QuantumCumulants = "0.4"
-QuestBase = "0.3.4"
+QuantumCumulants = "0.7"
+QuestBase = "0.4.1"
 Random = "1.10"
 RuntimeGeneratedFunctions = "0.5.5"
-SciMLBase = "2.83.1"
-SteadyStateDiffEq = "2.5"
-SymbolicUtils = "3.25"
-Symbolics = "6.34"
+SciMLBase = "3"
+SteadyStateDiffEq = "2.10"
+SymbolicUtils = "4"
+Symbolics = "7"
 Test = "1.10"
 TestExtras = "0.3"
-julia = "1.10"
+julia = "1.12"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
@@ -81,20 +84,24 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 HarmonicBalance = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+ModelingToolkitBase = "7771a370-6774-4173-bd38-47e70ca0b839"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+QuantumCumulants = "35bcea6d-e19f-57db-af74-8011de6c7255"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Peaks = "18e31ff7-3703-566c-8e60-38913d67486b"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-QuantumCumulants = "35bcea6d-e19f-57db-af74-8011de6c7255"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 
+[sources]
+HarmonicBalance = {url = "https://github.com/QuantumEngineeredSystems/HarmonicBalance.jl", rev = "source"}
+
 [targets]
-test = ["Aqua", "JET", "CheckConcreteStructs", "SymbolicUtils", "Symbolics", "ExplicitImports", "Test", "Random", "OrdinaryDiffEqTsit5", "Documenter", "Plots", "QuantumCumulants", "SteadyStateDiffEq", "LinearSolve", "TestExtras", "ModelingToolkit", "OrdinaryDiffEqRosenbrock", "NonlinearSolve", "HarmonicBalance", "Peaks"]
+test = ["Aqua", "JET", "CheckConcreteStructs", "SymbolicUtils", "Symbolics", "ExplicitImports", "Test", "Random", "OrdinaryDiffEqTsit5", "Documenter", "Plots", "SteadyStateDiffEq", "LinearSolve", "TestExtras", "ModelingToolkitBase", "OrdinaryDiffEqRosenbrock", "NonlinearSolve", "HarmonicBalance", "Peaks", "Latexify", "QuantumCumulants"]

docs/Project.toml

@@ -2,6 +2,7 @@
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
 HarmonicSteadyState = "1158f75c-a779-4b85-8bfb-8fcf6bf02ced"
+Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
@@ -13,4 +14,5 @@ SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 
 [compat]
 Documenter = "1"
-julia = "1.10"
+QuestBase = "0.4.1"
+julia = "1.12"

docs/make.jl

@@ -2,7 +2,7 @@ CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !==
 
 using HarmonicSteadyState
 using Documenter
-using Plots, SteadyStateDiffEq, OrdinaryDiffEqTsit5
+using Plots, Latexify, SteadyStateDiffEq, OrdinaryDiffEqTsit5
 
 TimeEvolution = Base.get_extension(HarmonicSteadyState, :TimeEvolution)
 SteadyStateDiffEqExt = Base.get_extension(HarmonicSteadyState, :SteadyStateDiffEqExt)
@@ -19,12 +19,14 @@ end
 makedocs(;
     sitename="HarmonicSteadyState.jl",
     authors="Quest group",
-    modules=[
-        HarmonicSteadyState,
-        TimeEvolution,
-        SteadyStateDiffEqExt,
-        HarmonicSteadyState.LinearResponse,
-        PlotsExt,
+    modules=Module[
+        m for m in [
+            HarmonicSteadyState,
+            TimeEvolution,
+            SteadyStateDiffEqExt,
+            HarmonicSteadyState.LinearResponse,
+            PlotsExt,
+        ] if !isnothing(m)
     ],
     format=Documenter.HTML(;
         canonical="https://quantumengineeredsystems.github.io/HarmonicSteadyState.jl/stable/",

examples/forward_transmission.jl

@@ -10,7 +10,7 @@ h = hm ⊗ hc  # Hilbertspace
 
 @qnumbers m::Destroy(h, 1) c::Destroy(h, 2) # Operators
 
-@rnumbers Δ Vk Ωd γm γk # Parameters
+@variables Δ::Real Vk::Real Ωd::Real γm::Real γk::Real # Parameters
 param = [Δ, Vk, Ωd, γm, γk]
 
 H_RWA_sym = (

ext/HarmonicBalanceExt.jl

@@ -63,7 +63,7 @@ function HarmonicSteadyState.LinearResponse.ResponseMatrix(
         args = cat(symbols, [Δ]; dims=1)
         f_re = Symbolics.build_function(el.re, args; expression=Val{false})
         f_im = Symbolics.build_function(el.im, args; expression=Val{false})
-        (args...) -> f_re(args...) + im * f_im(args...)
+        return (args...) -> f_re(args...) + im * f_im(args...)
     end
     return ResponseMatrix(compiled_M, symbols, res.problem.eom.variables)
 end

ext/PlotsExt/PlotsExt.jl

@@ -24,8 +24,8 @@ using HarmonicSteadyState.LinearResponse:
     get_jacobian_response, get_rotframe_jacobian_response
 
 using Plots: heatmap, theme_palette, scatter, RGB, cgrad
-using Symbolics.Latexify: Latexify, latexify
-using Symbolics.Latexify.LaTeXStrings: LaTeXStrings, @L_str
+using Latexify: Latexify, latexify
+using Latexify.LaTeXStrings: LaTeXStrings, @L_str
 using SciMLBase: SciMLBase
 
 const _set_Plots_default = Dict{Symbol,Any}([

ext/PlotsExt/linear_response.jl

@@ -203,7 +203,7 @@ function HarmonicSteadyState.plot_eigenvalues(
                 ),
             )
         end
-        eigvals(jac)
+        return eigvals(jac)
     end
     eigenvalues_filtered = map(.*, eigenvalues, filter_branch)
 

ext/QuantumCumulantsExt.jl

@@ -39,15 +39,17 @@ function replace_to_reals(eqs::Vector{String}, ops_av::Vector{String}, conjugate
     foreach(conjugate_pairs) do pair
         is_real = isequal(pair...)
         op_pair = ops_av[pair[1]]
-        op_name = replace(op_pair, "⟨" => "", "⟩" => "", "*" => "", "′" => "⁺")
+        op_name = replace(
+            op_pair, "⟨" => "", "⟩" => "", "*" => "", "′" => "⁺", "'" => "⁺", " " => ""
+        )
         op_r = op_name * "ᵣ"
         op_i = op_name * "ᵢ"
 
         push!(variables, op_r)
         is_real || push!(variables, op_i)
         conj_replace = map(eachindex(pair)) do i
             op = ops_av[pair[i]] # ∨ f(i) contains "+" or "-"
-            op => if is_real
+            return op => if is_real
                 "(" * op_r * ")"
             else
                 "(" * op_r * " $(f(i)) im*" * op_i * ")"
@@ -89,14 +91,14 @@ function compute_real_equations(eqs::MeanfieldEquations)
         eq_re = real(eq) # the factor of √(2) is needed to have Q to C limit.
         eq_im = - imag(eq)
         iszero(eq_re) || push!(eqs_real, eq_re)
-        iszero(eq_im) || push!(eqs_real, eq_im)
+        return iszero(eq_im) || push!(eqs_real, eq_im)
     end
     return eqs_real, Num.(vars)
 end # TODO: test if order of vars and eqs is correct
 
 function add_iv(eqs, vars, iv)
     var_new = map(vars) do var
-        QuestBase.declare_variable(string(var), iv)
+        return QuestBase.declare_variable(string(var), iv)
     end
     subs = Dict(zip(vars, var_new))
     eqs_new = QuestBase.substitute_all(eqs, subs)
@@ -126,12 +128,14 @@ function QuestBase.HarmonicEquation(MFeqs::MeanfieldEquations, parameters)
     jac = HarmonicSteadyState.get_Jacobian(eqs_jac, vars_jac)
 
     hvars = map(vars) do var
-        HarmonicSteadyState.QuestBase.HarmonicVariable(Num(var), "", "", Num(1), Num(0))
+        return HarmonicSteadyState.QuestBase.HarmonicVariable(
+            Num(var), "", "", Num(1), Num(0)
+        )
     end
 
     equations_lhs = map(enumerate(vars)) do (idx, var)
         dvar = QuestBase.d(var, MFeqs.iv)
-        equations_lhs[idx] ~ dvar # by convention lhs in HB
+        return equations_lhs[idx] ~ dvar # by convention lhs in HB
     end
 
     return QuestBase.HarmonicEquation(equations_lhs, hvars, Num.(parameters), jac)

ext/SteadyStateDiffEqExt.jl

@@ -31,7 +31,7 @@ function HarmonicSteadyState.steady_state_sweep(
         # make type-stable: FD.Dual or Float
         p = get_new_parameters(prob, varied_idx, value)
         sol = solve(remake(prob; p, u0), alg; kwargs...)
-        result[i] = sol.u
+        return result[i] = sol.u
     end
     return result
 end

ext/TimeEvolution/ODEProblem.jl

@@ -29,16 +29,15 @@ function OrdinaryDiffEqTsit5.ODEProblem(
     )
     vars = get_variables(eom)
 
-    # substitute the harmonic variables
-    eqs(v) = [substitute(eq, Dict(zip(vars, v))) for eq in subeqs]
-    # substitute  sweep parameters
+    # substitute the harmonic variables and sweep parameters; fold so the result is numeric
     function eqs(v, T)
-        return real.(
-            unwrap.([
-                substitute(eq, Dict(zip(keys(sweep), [sweep[p](T) for p in keys(sweep)])))
-                for eq in eqs(v)
-            ]),
+        sub_dict = merge(
+            Dict(zip(vars, v)), Dict(zip(keys(sweep), [sweep[p](T) for p in keys(sweep)]))
         )
+        return [
+            real(Symbolics.value(Symbolics.substitute(eq, sub_dict; fold=Val{true}()))) for
+            eq in subeqs
+        ]
     end
 
     function f!(du, u, p, T) # in-place

ext/TimeEvolution/TimeEvolution.jl

@@ -1,7 +1,7 @@
 module TimeEvolution
 
 using DocStringExtensions: TYPEDSIGNATURES
-using Symbolics: Num, substitute, unwrap, get_variables
+using Symbolics: Symbolics, Num, get_variables
 using OrdinaryDiffEqTsit5: OrdinaryDiffEqTsit5
 
 using HarmonicSteadyState:

src/HC_wrapper.jl

@@ -12,7 +12,11 @@ HC = HomotopyContinuation
 
 "Conversion from Symbolics.jl types to HomotopyContinuation types."
 function HomotopyContinuation.ModelKit.Variable(var::Num)
-    isterm(var.val) ? HC.Variable(string(var.val.f)) : HC.Variable(string(var_name(var)))
+    return if isterm(var.val)
+        HC.Variable(string(var.val.f))
+    else
+        HC.Variable(string(var_name(var)))
+    end
 end
 
 "Converts a Num into Variable in the active namespace."

src/LimitCycles/gauge_fixing.jl

@@ -14,7 +14,7 @@ Equivalent to adding `n` pairs of harmonics ω +- ω_lc for each existing ω.
 """
 function add_pairs!(eom::DifferentialEquation; ω_lc::Num, n::Int)
     foreach(1:n) do k
-        add_pairs!(eom, ω_lc)
+        return add_pairs!(eom, ω_lc)
     end
     return nothing
 end

src/LinearResponse/LinearResponse.jl

@@ -5,6 +5,7 @@ using DocStringExtensions: TYPEDFIELDS, TYPEDSIGNATURES, TYPEDEF
 using ProgressMeter: ProgressMeter, Progress, next!
 
 using Symbolics: Symbolics, Num, unwrap
+using SymbolicUtils: SymbolicUtils
 using LinearAlgebra: LinearAlgebra, norm, eigen, eigvals, eigvecs, I
 
 using HarmonicSteadyState:

src/LinearResponse/Lorentzian_spectrum.jl

@@ -47,7 +47,7 @@ end
 
 "Gives the numerical value of a peak at ω."
 function evaluate(peak::Lorentzian{T}, ω::T) where {T<:Real}
-    peak.A / sqrt(((peak.ω0 - ω)^2 + (peak.Γ)^2))
+    return peak.A / sqrt(((peak.ω0 - ω)^2 + (peak.Γ)^2))
 end
 
 "Gives the numerical value of a JacobianSpectrum at ω"
@@ -112,7 +112,11 @@ function JacobianSpectrum(
         for pair in _get_uv_pairs(hvars)
             u, v = hvars[pair]
             eigvec_2d = eigvec[pair] # fetch the relevant part of the Jacobian eigenvector
-            ωnum = real(unwrap(Symbolics.substitute(u.ω, Dict(solution_dict))))
+            ωnum = real(
+                SymbolicUtils.unwrap_const(
+                    Symbolics.unwrap(Symbolics.substitute(u.ω, Dict(solution_dict)))
+                ),
+            )
             # ^ the harmonic (numerical now) associated to this harmonic variable
 
             # eigvec_2d is associated to a natural variable -> this variable gets Lorentzian peaks