Fix SymbolicUtils v4+ compatibility issues

This commit addresses multiple API compatibility issues with SymbolicUtils v4.9+
and Symbolics v7.2+ where literal numbers are wrapped in symbolic containers.

Changes:
- rules.jl: Use recursive helpers for variadic rules instead of map/prod/sum
  which create symbolic mapreduce expressions. Add eval_numeric helper to
  evaluate numeric expressions like max(0,1). Fix poly_deg to return Julia
  numbers instead of symbolic wrappers. Add extract_numeric_val helper for
  decompose_rational.

- roots.jl: Add extract_numeric helper function to extract numeric values
  from symbolic substitution results in solve_newton and find_roots.

- integral.jl: Use isequal() instead of == for symbolic comparison, and
  === for nothing comparison to avoid symbolic boolean context errors.

- numeric_utils.jl: Comprehensive rewrite of accept_solution to properly
  handle both Num and BasicSymbolic types. Check for Num before Number
  (since Num <: Number). Extract numeric values from symbolic zero
  representations using Symbolics.value(Num(...)).

These changes ensure that the package works correctly with the new
SymbolicUtils v4+ internal representation where even literal numbers
are wrapped in BasicSymbolic containers.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

Author: claude

src/integral.jl

@@ -81,7 +81,7 @@ function integrate(
 
     eq = expand(eq)
 
-    if x == nothing
+    if x === nothing
         vars = get_variables(eq)
         if length(vars) > 1
             error("Multiple symbolic variables detect. Please pass the independent variable to `integrate`")
@@ -356,7 +356,7 @@ end
 function try_symbolic(eq, x, has_sym_consts = false, params = []; plan = default_plan())
     y = integrate_symbolic(eq, x; plan)
 
-    if y == nothing
+    if y === nothing
         if has_sym_consts && !isempty(params)
             @info("Symbolic integration failed. Try changing constant parameters ([$(join(params, ", "))]) to numerical values.")
         end

src/numeric_utils.jl

@@ -25,24 +25,78 @@ function accept_solution(eq, x, sol; plan = default_plan())
         # Δ = substitute(diff(sol, x) - expr(eq), Dict(x => x₀))
         S = subs_symbols(eq, x; include_x = true, plan.radius)
         # Also substitute any symbols in sol that may not be in eq
+        # Use value() to ensure consistent comparison with x (which may be BasicSymbolic)
+        x_val = value(x)
         for v in get_variables(value(sol))
-            if !haskey(S, v) && !isequal(v, x)
+            if !haskey(S, v) && !isequal(v, x_val)
                 S[v] = Complex(randn())
             end
         end
         Δ = substitute(diff(sol, x) - expr(eq), S)
+
+        # Check if Δ is zero - handle both Num and BasicSymbolic cases
+        # For Num, isequal works; for BasicSymbolic, we need to extract the value
+        if isequal(Δ, 0)
+            return 0.0
+        end
+        # Try to extract numeric value for BasicSymbolic zero
+        try
+            Δ_val = Symbolics.value(Num(Δ))
+            if Δ_val isa Real || Δ_val isa Complex
+                return abs(Δ_val)
+            end
+        catch
+        end
+
         result = abs(Δ)
-        # Ensure result is numeric, not symbolic
-        if result isa Number
-            return result
-        else
-            # Try to extract numeric value
+
+        # Note: Num <: Number, so we must check for Num BEFORE Number
+        # First try to extract numeric value from Num type
+        if result isa Num
+            inner = Symbolics.unwrap(result)
+            # Check if unwrapped value is a concrete number (not symbolic)
+            if inner isa Real || inner isa Complex
+                return abs(inner)
+            end
+            # Check if inner is structurally zero (e.g., abs(0))
+            if isequal(inner, 0) || (SymbolicUtils.iscall(inner) && SymbolicUtils.operation(inner) === abs &&
+                                     let arg = SymbolicUtils.arguments(inner)[1]
+                                         isequal(arg, 0) || (arg isa Real && arg == 0)
+                                     end)
+                return 0.0
+            end
+            # Try value extraction for symbolic wrapper (SymbolicUtils v4+)
+            try
+                val = Symbolics.value(result)
+                if val isa Real || val isa Complex
+                    return abs(val)
+                end
+            catch
+            end
+            # If still symbolic, return Inf
+            return Inf
+        end
+        # For non-Num concrete numbers
+        if result isa Real || result isa Complex
+            return abs(result)
+        end
+        # Try to extract numeric value from symbolic wrapper (SymbolicUtils v4+)
+        try
+            val = Symbolics.value(Num(result))
+            if val isa Real || val isa Complex
+                return abs(val)
+            end
+        catch
+        end
+        # Fallback: try unwrap
+        try
             val = Symbolics.unwrap(result)
-            if val isa Number
+            if val isa Real || val isa Complex
                 return abs(val)
             end
-            return Inf
+        catch
         end
+        return Inf
     catch e
         #
     end

src/roots.jl

@@ -1,3 +1,15 @@
+# Helper to extract numeric value from symbolic substitution result
+function extract_numeric(T, expr)
+    try
+        # Try to get the underlying value from symbolic wrapper
+        val = Symbolics.value(Num(expr))
+        return T(val)
+    catch
+        # Fall back to direct conversion
+        return T(expr)
+    end
+end
+
 # solve_newton is a symbolic Newton-Ralphson solver
 #   f is a symbolic equation to be solved (f ~ 0)
 #   x is the variable to solve
@@ -8,8 +20,8 @@ function solve_newton(T, p, ∂p, x, x₀, zs; abstol = 1.0e-10, maxiter = 50, s
 
     for i in 1:maxiter
         d[x] = xₙ
-        f = T(substitute(p, d))
-        f′ = T(substitute(∂p, d))
+        f = extract_numeric(T, substitute(p, d))
+        f′ = extract_numeric(T, substitute(∂p, d))
         ρ = sum(1 / (xₙ - z) for z in zs; init = 0)
         xₙ₊₁ = xₙ - s * f / (f′ - s * ρ * f)
 
@@ -60,7 +72,7 @@ function find_roots(T, p, x; abstol = 1.0e-8, num_roots = 0)
                 push!(s, z)
 
                 z = conj(z)
-                if abs(Complex{T}(substitute(p, Dict(x => z)))) < abstol
+                if abs(extract_numeric(Complex{T}, substitute(p, Dict(x => z)))) < abstol
                     push!(zs, z)
                     push!(s, z)
                 end

src/rules.jl

@@ -1,5 +1,15 @@
 ########################## Predicates #########################################
 
+# Helper to extract numeric value from symbolic substitution result
+function extract_numeric_val(T, expr)
+    try
+        val = Symbolics.value(Num(expr))
+        return T(val)
+    catch
+        return T(expr)
+    end
+end
+
 is_pos_int(x) = is_proper(x) && isinteger(x) && x > 0
 is_neg_int(x) = is_proper(x) && isinteger(x) && x < 0
 is_int_gt_one(x) = is_proper(x) && isinteger(x) && x > 1
@@ -54,10 +64,15 @@ is_var(eq) = isequal(eq, var(eq))
 
 @syms Ω(x)
 
+# Helper to compute max of applying Ω to a list of expressions
+apply_max_omega(xs) = length(xs) == 1 ? Ω(xs[1]) : max(Ω(xs[1]), apply_max_omega(xs[2:end]))
+# Helper to compute sum of applying Ω to a list of expressions
+apply_sum_omega(xs) = length(xs) == 1 ? Ω(xs[1]) : Ω(xs[1]) + apply_sum_omega(xs[2:end])
+
 # rules to calculate the degree of a polynomial expression or NaN if not a poly
 p_rules = [
-    @rule Ω(+(~~xs)) => max(map(Ω, ~~xs)...)
-    @rule Ω(~x * ~y) => Ω(~x) + Ω(~y)
+    @rule Ω(+(~~xs)) => apply_max_omega(~~xs)
+    @rule Ω(*(~~xs)) => apply_sum_omega(~~xs)
     @rule Ω(^(~x, ~k::is_proper)) => is_pos_int(~k) ? ~k * Ω(~x) : NaN
     @rule Ω(~x::is_number) => 0
     @rule Ω(~x::is_var) => 1
@@ -66,16 +81,66 @@ p_rules = [
 
 function poly_deg(eq)
     ω = Prewalk(Chain(p_rules))(Ω(value(eq)))
-    return substitute(ω, Dict())
+    result = substitute(ω, Dict())
+    # Evaluate numeric expressions like max(0, 1) -> 1
+    result = eval_numeric(result)
+    # Convert symbolic numbers to Julia numbers
+    # This is needed because SymbolicUtils v4+ wraps literal numbers in symbolic containers
+    try
+        return Symbolics.value(Num(result))
+    catch
+        return result  # Return as-is if conversion fails
+    end
 end
 
-is_poly(eq) = !isnan(poly_deg(eq))
-is_linear_poly(eq) = isequal(poly_deg(eq), 1)
+# Helper function to evaluate purely numeric symbolic expressions
+function eval_numeric(expr)
+    # First try to extract as Julia number
+    if expr isa Number
+        return expr
+    end
+
+    if !SymbolicUtils.iscall(expr)
+        # Not a call - might be a symbolic literal number or symbol
+        # Try to extract value
+        try
+            return Symbolics.value(Num(expr))
+        catch
+            return expr
+        end
+    end
+
+    op = SymbolicUtils.operation(expr)
+    args = SymbolicUtils.arguments(expr)
+
+    # Recursively evaluate arguments
+    evaluated_args = [eval_numeric(a) for a in args]
+
+    # Check if all arguments are Julia numbers now
+    if all(a -> a isa Number, evaluated_args)
+        # Evaluate the operation
+        try
+            return op(evaluated_args...)
+        catch
+            return expr
+        end
+    end
+
+    return expr
+end
+
+is_poly(eq) = let deg = poly_deg(eq); is_number(deg) && !isnan(deg) end
+is_linear_poly(eq) = poly_deg(eq) == 1
+
+# Helper to distribute Ω over a list with sum
+apply_sum_omega_kernel(xs) = length(xs) == 1 ? Ω(xs[1]) : Ω(xs[1]) + apply_sum_omega_kernel(xs[2:end])
+# Helper to distribute Ω over a list with product
+apply_prod_omega_kernel(xs) = length(xs) == 1 ? Ω(xs[1]) : Ω(xs[1]) * apply_prod_omega_kernel(xs[2:end])
 
 # rules to extract the kernel (the holomorphic portion) of an expression
 s_rules = [
-    @rule Ω(+(~~xs)) => sum(map(Ω, ~~xs))
-    @rule Ω(*(~~xs)) => prod(map(Ω, ~~xs))
+    @rule Ω(+(~~xs)) => apply_sum_omega_kernel(~~xs)
+    @rule Ω(*(~~xs)) => apply_prod_omega_kernel(~~xs)
     # @rule Ω(~x / ~y) => Ω(~x) * Ω(^(~y, -1))
     @rule Ω(^(~x::is_linear_poly, ~k::is_pos_int)) => ^(~x, ~k)
     @rule Ω(~x::is_var) => ~x
@@ -264,9 +329,9 @@ function decompose_rational(eq)
     for i in 1:n
         x₀ = test_point(false, 1.0)
         d = Dict(x => x₀)
-        b[i] = 1 / substitute(eq, d)
+        b[i] = 1 / extract_numeric_val(Complex, substitute(eq, d))
         for j in 1:n
-            A[i, j] = substitute(F[j], d)
+            A[i, j] = extract_numeric_val(Complex, substitute(F[j], d))
         end
     end
 
@@ -276,9 +341,14 @@ function decompose_rational(eq)
     return p
 end
 
+# Helper for rational rules - distribute Ω over sum
+apply_sum_omega_rat(xs) = length(xs) == 1 ? Ω(xs[1]) : Ω(xs[1]) + apply_sum_omega_rat(xs[2:end])
+# Helper for rational rules - distribute Ω over product
+apply_prod_omega_rat(xs) = length(xs) == 1 ? Ω(xs[1]) : Ω(xs[1]) * apply_prod_omega_rat(xs[2:end])
+
 rational_rules = [
-    @rule Ω(+(~~xs)) => sum(map(Ω, ~~xs))
-    @rule Ω(*(~~xs)) => prod(map(Ω, ~~xs))
+    @rule Ω(+(~~xs)) => apply_sum_omega_rat(~~xs)
+    @rule Ω(*(~~xs)) => apply_prod_omega_rat(~~xs)
     @rule Ω(^(~x::is_poly, ~k::is_neg_int)) => expand(
         ^(
             decompose_rational(~x),