feat: update to symbolics V7 and SymbolicUtils v4

Author: oameye

.github/workflows/Tests.yml

@@ -40,7 +40,6 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - 'pre'
           - 'lts'
           - '1'
         os:

Project.toml

@@ -21,7 +21,7 @@ Test = "1.10"
 OrderedCollections = "1.8"
 Aqua = "0.8.11"
 ExplicitImports = "1.11"
-JET = "0.9.18, 0.10.0"
+JET = "0.9.18, 0.10.0, 0.11"
 CheckConcreteStructs = "0.1.0"
 
 [extras]

src/DifferentialEquation.jl

@@ -163,7 +163,10 @@ corresponding to second-order differential equations.
 function is_rearranged_standard(eom::DifferentialEquation, degree=2)
     tvar = get_independent_variables(eom)[1]
     D = Differential(tvar)^degree
-    return isequal(getfield.(values(eom.equations), :lhs), D.(get_variables(eom)))
+    lhs = getfield.(values(eom.equations), :lhs)
+    rhs = D.(get_variables(eom))
+    diffs = Symbolics.simplify.(lhs .- rhs)
+    return all(is_literal_zero, diffs)
 end
 
 """
@@ -195,7 +198,17 @@ function rearrange!(eom::DifferentialEquation, new_lhs::Vector{Num})
     return nothing
 end
 function get_variables_nums(vars::Vector{Num})
-    unique(flatten([Num.(get_variables(x)) for x in vars]))
+    # Symbolics v7: `get_variables(Differential(t, n)(x(t)))` returns the derivative term
+    # itself, so we must explicitly strip derivatives to recover the dependent variable.
+    out = Num[]
+    for expr in vars
+        sym = Symbolics.unwrap(expr)
+        while Symbolics.is_derivative(sym)
+            sym = first(Symbolics.arguments(sym))
+        end
+        push!(out, Num(sym))
+    end
+    return out
 end # TODO: remove this function or at least better names
 
 """

src/HarmonicEquation.jl

@@ -62,11 +62,23 @@ end
 
 "Get the parameters (not time nor variables) of a HarmonicEquation"
 function _parameters(eom::HarmonicEquation)
-    all_symbols = flatten([
-        cat(get_variables(eq.lhs), get_variables(eq.rhs); dims=1) for eq in eom.equations
-    ])
-    # subtract the set of independent variables (i.e., time) from all free symbols
-    return setdiff(all_symbols, get_variables(eom), get_independent_variables(eom))
+    symbols = Num[]
+    for eq in eom.equations
+        vars = union(Symbolics.get_variables(eq.lhs), Symbolics.get_variables(eq.rhs))
+        vars = sort!(collect(vars); by=string)
+        for sym in vars
+            if Symbolics.is_derivative(sym)
+                sym = first(Symbolics.arguments(sym))
+            end
+            push!(symbols, Num(sym))
+        end
+    end
+    vars = Set(get_variables(eom))
+    indep = Set(get_independent_variables(eom))
+    params = filter(s -> !(s in vars || s in indep), symbols)
+    params = unique(params)
+    sort!(params; by=string)
+    return params
 end
 
 """

src/Symbolics/Symbolics_utils.jl

@@ -2,6 +2,13 @@
 expand_all(x::Num) = Num(expand_all(x.val))
 _apply_termwise(f, x::Num) = wrap(_apply_termwise(f, unwrap(x)))
 
+# Symbolics v7 note: `unwrap` now always returns `BasicSymbolic`; use `Symbolics.value`
+# when you want the underlying literal number (via `Const`).
+_literal_number(x::Num) = (v = Symbolics.value(x); v isa Number ? v : nothing)
+_literal_number(x::BasicSymbolic) =
+    (v = SymbolicUtils.unwrap_const(x); v isa Number ? v : nothing)
+is_literal_zero(x) = (v = _literal_number(x); v !== nothing && iszero(v))
+
 "Expands using SymbolicUtils.expand and expand_exp_power (changes exp(x)^n to exp(x*n)"
 function expand_all(x)
     result = Postwalk(expand_exp_power)(SymbolicUtils.expand(x))
@@ -15,7 +22,14 @@ function expand_fraction(x::BasicSymbolic)
     elseif ismul(x)
         return _apply_termwise(expand_fraction, x)
     elseif isdiv(x)
-        return sum([arg / x.den for arg in arguments(x.num)])
+        # Only distribute division over addition in the numerator.
+        # In SymbolicUtils v4, `arguments(x.num)` for a multiplication returns factors,
+        # so splitting unconditionally would produce incorrect results (e.g. d*(a+b)/c).
+        num, den = x.num, x.den
+        if isadd(num)
+            return sum([expand_fraction(arg) / den for arg in arguments(num)])
+        end
+        return expand_fraction(num) / expand_fraction(den)
     else
         return x
     end
@@ -131,7 +145,7 @@ function is_harmonic(x::Num, t::Num)::Bool
     isempty(t_terms) && return true
     trigs = is_trig.(t_terms)
 
-    if !prod(trigs)
+        if !all(trigs)
         return false
     else
         powers = [max_power(first(arguments(term.val)), t) for term in t_terms[trigs]]

src/Symbolics/drop_powers.jl

@@ -17,10 +17,14 @@ function drop_powers(expr::Num, vars::Vector{Num}, deg::Int)
     Symbolics.@variables ϵ
     subs_expr = deepcopy(expr)
     rules = Dict([var => ϵ * var for var in unique(vars)])
-    subs_expr = Symbolics.expand(substitute_all(subs_expr, rules))
+    subs_expr = Symbolics.expand(substitute_all(subs_expr, rules; include_derivatives=false))
     max_deg = max_power(subs_expr, ϵ)
     removal = Dict([ϵ^d => Num(0) for d in deg:max_deg])
-    res = substitute_all(substitute_all(subs_expr, removal), Dict(ϵ => Num(1)))
+    res = substitute_all(
+        substitute_all(subs_expr, removal; include_derivatives=false),
+        Dict(ϵ => Num(1));
+        include_derivatives=false,
+    )
     return Symbolics.expand(res)
 end
 
@@ -30,7 +34,7 @@ end
 
 # calls the above for various types of the first argument
 function drop_powers(eq::Equation, var::Vector{Num}, deg::Int)
-    return drop_powers(eq.lhs, var, deg) .~ drop_powers(eq.lhs, var, deg)
+    return drop_powers(eq.lhs, var, deg) .~ drop_powers(eq.rhs, var, deg)
 end
 function drop_powers(eqs::Vector{Equation}, var::Vector{Num}, deg::Int)
     return [
@@ -45,9 +49,12 @@ drop_powers(x, vars, deg::Int) = drop_powers(wrap(x), vars, deg)
 function max_power(x::Num, y::Num)
     terms = get_all_terms(x)
     powers = power_of.(terms, y)
-    literal_powers = Int[
-        Int(SymbolicUtils.unwrap_const(p)) for p in powers if SymbolicUtils.is_literal_number(p)
-    ]
+    literal_powers = Int[]
+    for p in powers
+        pv = SymbolicUtils.unwrap_const(p)
+        pv isa Number || continue
+        push!(literal_powers, Int(pv))
+    end
     isempty(literal_powers) && return 0
     return maximum(literal_powers)
 end

src/Symbolics/exponentials.jl

@@ -43,8 +43,9 @@ function simplify_exp_products_mul(expr)
     rest_ind = setdiff(1:length(arguments(expr)), ind)
     rest = isempty(rest_ind) ? 1 : prod(arguments(expr)[rest_ind])
     total = isempty(ind) ? 0 : sum(getindex.(arguments.(arguments(expr)[ind]), 1))
-    if SymbolicUtils.is_literal_number(total)
-        return iszero(SymbolicUtils.unwrap_const(total)) ? rest : rest * exp(total)
+    total_val = SymbolicUtils.unwrap_const(total)
+    if total_val isa Number
+        return iszero(total_val) ? rest : rest * exp(total)
     end
     return rest * exp(total)
 end