build: support SymbolicUtils 4 / Symbolics 7

.github/codecov.yml

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

.github/workflows/Documentation.yml

@@ -37,7 +37,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

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

Project.toml

@@ -14,11 +14,11 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 DocStringExtensions = "0.9.4"
 SymbolicUtils = "4"
 Symbolics = "7"
-julia = "1.10"
+julia = "1.12"
 Random = "1.10"
 LinearAlgebra = "1.10"
 Test = "1.10"
-OrderedCollections = "1.8"
+OrderedCollections = "1.8, 2"
 Aqua = "0.8.11"
 ExplicitImports = "1.11"
 JET = "0.9.18, 0.10.0, 0.11"

src/DifferentialEquation.jl

@@ -65,7 +65,7 @@ mutable struct DifferentialEquation
     function DifferentialEquation(eq::Equation, vars::Vector{Num})
         typerhs = typeof(eq.rhs)
         typelhs = typeof(eq.lhs)
-        throw(
+        return throw(
             ArgumentError(
                 "The variables are of type $(typeof(vars)). Whereas you gave one equation of
                 type $(typerhs)~$(typelhs). Commenly one forgot to broadcast the equation symbol `~`.",
@@ -93,9 +93,8 @@ $(TYPEDSIGNATURES)
 
 Return the dependent variables of `diff_eom`.
 """
-Symbolics.get_variables(diff_eom::DifferentialEquation)::Vector{Num} = collect(
-    keys(diff_eom.equations)
-)
+Symbolics.get_variables(diff_eom::DifferentialEquation)::Vector{Num} =
+    collect(keys(diff_eom.equations))
 
 """
 $(TYPEDSIGNATURES)
@@ -104,9 +103,8 @@ Check if all equations in `diff_eom` are harmonic with respect to `t`. The funct
 differential equation system `diff_eom` and a variable `t`, and returns `true` if all equations
 are harmonic with respect to `t`, otherwise it returns `false`.
 """
-QuestBase.is_harmonic(diff_eom::DifferentialEquation, t::Num)::Bool = all([
-    is_harmonic(eq, t) for eq in values(diff_eom.equations)
-])
+QuestBase.is_harmonic(diff_eom::DifferentialEquation, t::Num)::Bool =
+    all([is_harmonic(eq, t) for eq in values(diff_eom.equations)])
 
 """
 $(TYPEDSIGNATURES)

src/HarmonicEquation.jl

@@ -174,6 +174,12 @@ function rearrange!(eom::HarmonicEquation, new_rhs::Vector{Num})
     soln = Symbolics.symbolic_linear_solve(
         eom.equations, new_rhs; simplify=false, check=true
     )
+    # `symbolic_linear_solve` returns nested fractions whose denominators hold the
+    # coefficient determinant; for a trigonometric ansatz it contains identities like
+    # cos(ωt)² + sin(ωt)² = 1 that SymbolicUtils 4's `simplify` no longer collapses.
+    # Flatten the nest (polynomial-level, no AC matching) and reduce the denominator,
+    # so no consumer works with superficially time-dependent denominators.
+    soln = reduce_denominator.(Symbolics.simplify_fractions.(Num.(soln)))
     eom.equations = soln .~ new_rhs
     return nothing
 end

src/HarmonicVariable.jl

@@ -110,6 +110,5 @@ end
 "Returns the symbols of a `HarmonicVariable`."
 Symbolics.get_variables(var::HarmonicVariable)::Num = Num(first(get_variables(var.symbol)))
 
-Base.isequal(v1::HarmonicVariable, v2::HarmonicVariable)::Bool = isequal(
-    v1.symbol, v2.symbol
-)
+Base.isequal(v1::HarmonicVariable, v2::HarmonicVariable)::Bool =
+    isequal(v1.symbol, v2.symbol)

src/QuestBase.jl

@@ -27,7 +27,6 @@ using Symbolics:
     Equation,
     Differential,
     arguments,
-    substitute,
     term,
     expand,
     operation,

src/Symbolics/Symbolics_utils.jl

@@ -63,8 +63,28 @@ Perform substitutions in `rules` on `x`.
 `Symbolics.substitute_in_deriv`.
 """
 Subtype = Union{Num,Equation,BasicSymbolic}
+
+# SymbolicUtils 4 / Symbolics 7: substitute() does not recurse into call arguments
+# (e.g. it leaves x(t) unchanged when substituting t => T). Walk the expression and
+# rewrite matching nodes ourselves.
+function _deep_substitute(e::BasicSymbolic, unwrap_rules::Dict)
+    # SU 4's default filter blocks recursion into the arguments of callable-symbolic
+    # terms (e.g. `x(t)`). Override with an always-true filter so substitutions like
+    # `t => T` reach inside `x(t)`. Compound keys still match before recursion, and
+    # SU does not re-walk the replacement, so self-referential rules don't loop.
+    return SymbolicUtils.substitute(e, unwrap_rules; filterer=_ -> true)
+end
+_deep_substitute(x::Num, ur::Dict) = wrap(_deep_substitute(unwrap(x), ur))
+function _deep_substitute(eq::Equation, ur::Dict)
+    return Equation(
+        _deep_substitute(unwrap(eq.lhs), ur), _deep_substitute(unwrap(eq.rhs), ur)
+    )
+end
+_deep_substitute(x, ::Dict) = x
+
 function substitute_all(x::Subtype, rules::Dict; include_derivatives=true)
-    result = substitute(x, rules)
+    unwrap_rules = Dict(unwrap(k) => unwrap(v) for (k, v) in rules)
+    result = _deep_substitute(x, unwrap_rules)
     if include_derivatives
         result = Symbolics.substitute_in_deriv(result, rules)
     end

src/Symbolics/fourier.jl

@@ -29,16 +29,47 @@ function trig_reduce(x)
     x = simplify_exp_products(x) # simplify products of exps
     x = exp_to_trig(x)
     x = Num(simplify_complex(expand(x)))
-    return x # simplify_fractions(x)# (a*c^2 + b*c)/c^2 = (a*c + b)/c
+    return reduce_denominator(x) # collapse trig identities left in the denominator
+end
+
+"""
+    reduce_denominator(x)
+
+Reduce trigonometric identities left in the denominator of a fraction.
+
+Fraction denominators arise from combining fractions (`add_div` merges `a/b + c/d`
+into a single fraction) and from solving linear systems (`rearrange!` divides by the
+coefficient determinant). The exponential round-trip in [`trig_reduce`](@ref) only
+reaches the numerator, so identities such as the Pythagorean
+`cos(ωt)² + sin(ωt)² => 1` survive in the denominator. Reduce the denominator on
+its own (it holds no nested fraction, so the recursion terminates) and leave the
+numerator untouched. Trig-free denominators are returned unchanged.
+"""
+function reduce_denominator(x)
+    u = unwrap(x)
+    isdiv(u) || return x
+    num, den = arguments(u)
+    contains_trig(den) || return x
+    return wrap(num) / trig_reduce(wrap(den))
+end
+
+"Return true if any subterm of `x` is a sin or cos."
+contains_trig(x::Num) = contains_trig(unwrap(x))
+contains_trig(x) = false
+function contains_trig(x::BasicSymbolic)
+    is_trig(x) && return true
+    SymbolicUtils.iscall(x) || return false
+    return any(contains_trig, arguments(x))
 end
 
 "Return true if `f` is a sin or cos."
 is_trig(f::Num) = is_trig(unwrap(f))
 is_trig(f) = false
 function is_trig(f::BasicSymbolic)
     f = ispow(f) ? arguments(f)[1] : f
-    isterm(f) && operation(f) ∈ [cos, sin] && return true
-    return false
+    isterm(f) || return false
+    op = operation(f)
+    return op === cos || op === sin
 end
 
 """

src/utils.jl

@@ -37,4 +37,8 @@ end
 is_identity(A::Matrix{Num}) = (@eqsym A == Matrix{Num}(LinearAlgebra.I, size(A)...))
 hasnan(x::Matrix{Num}) = any(my_isnan, unwrap.(x))
 my_isnan(x) = isnan(x)
-my_isnan(x::BasicSymbolic) = false
+function my_isnan(x::BasicSymbolic)
+    SymbolicUtils.isconst(x) || return false
+    v = Symbolics.value(x)
+    return v isa Number && isnan(v)
+end

test/symbolics.jl

@@ -152,6 +152,51 @@ end
     @eqtest fourier_cos_term((cos(f * t)^2 + cos(f * t))^3, 0, t) == 23//16
 end
 
+@testset "trig_reduce" begin
+    using QuestBase: trig_reduce
+
+    @variables f t a ω
+
+    # the Pythagorean identity collapses in a bare expression …
+    @eqtest trig_reduce(cos(f * t)^2 + sin(f * t)^2) == 1
+
+    # … and, crucially, in the denominator of a combined fraction. The exponential
+    # round-trip only touches the numerator, so without `reduce_denominator` the
+    # identity survives in the denominator, `get_independent` sees a time-dependent
+    # denominator and discards the whole fraction, collapsing the Krylov-Bogoliubov
+    # slow-flow equations to `0 ~ d/dT`.
+    @eqtest trig_reduce(a / (cos(f * t)^2 + sin(f * t)^2)) == a
+    @eqtest trig_reduce(a / (ω * cos(f * t)^2 + ω * sin(f * t)^2)) == a / ω
+
+    # trig-free denominators pass through `reduce_denominator` untouched
+    using QuestBase: reduce_denominator
+    @eqtest reduce_denominator(a / (ω^2 + 1)) == a / (ω^2 + 1)
+end
+
+@testset "rearrange! reduces determinant denominators" begin
+    # `symbolic_linear_solve` returns nested fractions whose denominators hold the
+    # coefficient determinant; for a trigonometric ansatz it contains cos² + sin² = 1.
+    # `rearrange!` must flatten and collapse it, otherwise every consumer works with
+    # superficially time-dependent denominators (this stalled the Krylov-Bogoliubov
+    # order-2 averaging for hours).
+    using QuestBase: HarmonicEquation, HarmonicVariable, DifferentialEquation
+    using QuestBase: rearrange!, d
+    @variables ω t T u1(T) v1(T) x(t)
+
+    eqs = [
+        cos(ω * t) * d(u1, T) + sin(ω * t) * d(v1, T) ~ u1,
+        -sin(ω * t) * d(u1, T) + cos(ω * t) * d(v1, T) ~ v1,
+    ]
+    hvars = [HarmonicVariable(u1, "u", "u", ω, x), HarmonicVariable(v1, "v", "v", ω, x)]
+    natural_eq = DifferentialEquation(d(x, t, 2) + x ~ 0, x)
+    eom = HarmonicEquation(eqs, hvars, natural_eq)
+
+    rearrange!(eom, d([u1, v1], T))
+    lhs = QuestBase.Num.(getfield.(eom.equations, :lhs))
+    @eqtest lhs[1] == u1 * cos(ω * t) - v1 * sin(ω * t)
+    @eqtest lhs[2] == u1 * sin(ω * t) + v1 * cos(ω * t)
+end
+
 @testset "_apply_termwise" begin
     using QuestBase: _apply_termwise