Merge pull request #1025 from ChrisRackauckas-Claude/bump-symbolics-mtk-compat

Bump Symbolics v7, SymbolicUtils v4, ModelingToolkit v11 compat

Author: ChrisRackauckas

Project.toml

@@ -85,8 +85,7 @@ LuxCore = "1.0.1"
 LuxLib = "1.3"
 MCMCChains = "7"
 MLDataDevices = "1.2.0"
-MethodOfLines = "0.11.6"
-ModelingToolkit = "9.46, 10"
+ModelingToolkit = "11"
 MonteCarloMeasurements = "1.1"
 NeuralOperators = "0.5, 0.6"
 Optimisers = "0.3.3, 0.4"
@@ -105,8 +104,8 @@ SciMLBase = "2.56"
 Statistics = "1.10"
 StochasticDiffEq = "6.69.1"
 SymbolicIndexingInterface = "0.3.31"
-SymbolicUtils = "3.7.2"
-Symbolics = "6.14"
+SymbolicUtils = "4.12"
+Symbolics = "7.8"
 TensorBoardLogger = "0.1.24"
 Test = "1.10"
 WeightInitializers = "1.0.3"
@@ -127,7 +126,6 @@ LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LuxCUDA = "d0bbae9a-e099-4d5b-a835-1c6931763bda"
 LuxCore = "bb33d45b-7691-41d6-9220-0943567d0623"
 LuxLib = "82251201-b29d-42c6-8e01-566dec8acb11"
-MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
@@ -136,4 +134,4 @@ TensorBoardLogger = "899adc3e-224a-11e9-021f-63837185c80f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Boltz", "CUDA", "DiffEqNoiseProcess", "ExplicitImports", "FastGaussQuadrature", "Flux", "Hwloc", "InteractiveUtils", "LineSearches", "LuxCUDA", "LuxCore", "LuxLib", "MethodOfLines", "OptimizationOptimJL", "OrdinaryDiffEq", "ReTestItems", "StochasticDiffEq", "TensorBoardLogger", "Test"]
+test = ["Aqua", "Boltz", "CUDA", "DiffEqNoiseProcess", "ExplicitImports", "FastGaussQuadrature", "Flux", "Hwloc", "InteractiveUtils", "LineSearches", "LuxCUDA", "LuxCore", "LuxLib", "OptimizationOptimJL", "OrdinaryDiffEq", "ReTestItems", "StochasticDiffEq", "TensorBoardLogger", "Test"]

docs/src/tutorials/constraints.md

@@ -22,7 +22,6 @@ with Physics-Informed Neural Networks.
 
 ```@example fokkerplank
 using NeuralPDE, Lux, ModelingToolkit, Optimization, OptimizationOptimJL, LineSearches
-using Integrals, Cubature
 using DomainSets: Interval
 using IntervalSets: leftendpoint, rightendpoint
 # the example is taken from this article https://arxiv.org/abs/1910.10503
@@ -52,15 +51,19 @@ chain = Lux.Chain(Dense(1, inn, Lux.σ),
     Dense(inn, inn, Lux.σ),
     Dense(inn, 1))
 
-lb = [x_0]
-ub = [x_end]
+lb = x_0
+ub = x_end
+# Use a simple trapezoidal rule for the normalization constraint.
+# This avoids AD issues with Integrals.jl's C-based quadrature solvers.
+norm_xs = collect(range(lb, ub, length = 200))
+norm_dx = Float64(norm_xs[2] - norm_xs[1])
 function norm_loss_function(phi, θ, p)
-    function inner_f(x, θ)
-        0.01 * phi(x, θ) .- 1
+    # Evaluate phi at quadrature points (each point as a 1-element vector)
+    s = sum(1:length(norm_xs)) do i
+        first(phi([norm_xs[i]], θ))
     end
-    prob = IntegralProblem(inner_f, lb, ub, θ)
-    norm2 = solve(prob, HCubatureJL(), reltol = 1e-8, abstol = 1e-8, maxiters = 10)
-    abs(norm2[1])
+    norm_val = 0.01 * s * norm_dx
+    abs(norm_val - 1)
 end
 
 discretization = PhysicsInformedNN(chain,

src/NeuralPDE.jl

@@ -36,8 +36,8 @@ using WeightInitializers: glorot_uniform, zeros32
 using Zygote: Zygote
 
 # Symbolic Stuff
-using ModelingToolkit: ModelingToolkit, PDESystem, Differential, toexpr, defaults
-using Symbolics: Symbolics, unwrap, arguments, operation, build_expr, Num,
+using ModelingToolkit: ModelingToolkit, PDESystem, Differential, toexpr
+using Symbolics: Symbolics, unwrap, arguments, operation, Num,
     expand_derivatives
 using SymbolicUtils: SymbolicUtils
 using SymbolicIndexingInterface: SymbolicIndexingInterface

src/PDE_BPINN.jl

@@ -85,19 +85,30 @@ function get_symbols(dataset, depvars, eqs)
     # order of pinnrep.depvars, depvar_vals, BayesianPINN.dataset must be same
     to_subs = Dict(depvars .=> depvar_vals)
 
-    numform_vars = Symbolics.get_variables.(eqs)
-    Eq_vars = unique(reduce(vcat, numform_vars))
-    # got equation's depvar num format {x(t)} for use in substitute()
-
+    # Find callable depvar terms (e.g., u(x,t)) by searching the equation tree.
+    # In SymbolicUtils v4, get_variables may decompose u(x,t) into bare u, x, t,
+    # so we search directly for callable terms whose operation matches a depvar name.
+    depvar_set = Set(depvars)
     tobe_subs = Dict()
-    for a in depvars
-        for i in Eq_vars
-            expr = toexpr(i)
-            if (expr isa Expr) && (expr.args[1] == a)
-                tobe_subs[a] = i
+    function _search(term)
+        t = SymbolicUtils.unwrap(term)
+        return if SymbolicUtils.iscall(t)
+            op = SymbolicUtils.operation(t)
+            if !SymbolicUtils.iscall(op) && !(op isa Differential)
+                name = toexpr(op)
+                if name isa Symbol && name in depvar_set
+                    tobe_subs[name] = t
+                end
+            end
+            for arg in SymbolicUtils.arguments(t)
+                _search(arg)
             end
         end
     end
+    for eq in eqs
+        _search(eq.lhs)
+        _search(eq.rhs)
+    end
     # depvar symbolic and num format got, tobe_subs : Dict{Any, Any}(:y => y(t), :x => x(t))
 
     return to_subs, tobe_subs

src/discretize.jl

@@ -403,7 +403,7 @@ For more information, see `discretize` and `PINNRepresentation`.
 function SciMLBase.symbolic_discretize(pde_system::PDESystem, discretization::AbstractPINN)
     (; eqs, bcs, domain) = pde_system
     eq_params = pde_system.ps
-    defaults = pde_system.defaults
+    defaults = pde_system.initial_conditions
     (;
         chain, param_estim, additional_loss, multioutput, init_params, phi,
         derivative, strategy, logger, iteration, self_increment,
@@ -412,7 +412,7 @@ function SciMLBase.symbolic_discretize(pde_system::PDESystem, discretization::Ab
     adaloss = discretization.adaptive_loss
 
     default_p = eq_params isa SciMLBase.NullParameters ? nothing :
-        [defaults[ep] for ep in eq_params]
+        [Symbolics.value(defaults[ep]) for ep in eq_params]
 
     depvars, indvars, dict_indvars,
         dict_depvars, dict_depvar_input = get_vars(

src/symbolic_utilities.jl

@@ -1,5 +1,12 @@
 using Base.Broadcast
 
+# build_expr was removed from Symbolics.jl v7; define locally
+function build_expr(head::Symbol, args)
+    ex = Expr(head)
+    append!(ex.args, args)
+    return ex
+end
+
 """
 Override `Broadcast.__dot__` with `Broadcast.dottable(x::Function) = true`
 
@@ -155,8 +162,11 @@ function _transform_expression(
                 derivative_variables = Symbol[]
                 order = 0
                 while (_args[1] isa Differential)
-                    order += 1
-                    push!(derivative_variables, toexpr(_args[1].x))
+                    d_order = _args[1].order
+                    order += d_order
+                    for _ in 1:d_order
+                        push!(derivative_variables, toexpr(_args[1].x))
+                    end
                     _args = _args[2].args
                 end
                 depvar = _args[1]
@@ -197,7 +207,19 @@ function _transform_expression(
                     integrating_var_id = [dict_indvars[i] for i in integrating_variable]
                 else
                     integrating_variable = toexpr(_args[1].domain.variables)
-                    integrating_var_id = [dict_indvars[integrating_variable]]
+                    # In SymbolicUtils v4, a symbolic tuple may produce
+                    # Expr(:call, :tuple, :x, :y) instead of a Julia Tuple
+                    if integrating_variable isa Expr &&
+                            integrating_variable.head == :call &&
+                            integrating_variable.args[1] === tuple
+                        integrating_variable = integrating_variable.args[2:end]
+                        integrating_var_id = [
+                            dict_indvars[i]
+                                for i in integrating_variable
+                        ]
+                    else
+                        integrating_var_id = [dict_indvars[integrating_variable]]
+                    end
                 end
 
                 integrating_depvars = []
@@ -336,8 +358,10 @@ Parse ModelingToolkit equation form to the inner representation.
        (derivative(phi2, u2, [x, y], [[ε,0]], 1, θ2) + 9 * derivative(phi1, u, [x, y], [[0,ε]], 1, θ1)) - 0]
 """
 function parse_equation(pinnrep::PINNRepresentation, eq)
-    eq_lhs = isequal(expand_derivatives(eq.lhs), 0) ? eq.lhs : expand_derivatives(eq.lhs)
-    eq_rhs = isequal(expand_derivatives(eq.rhs), 0) ? eq.rhs : expand_derivatives(eq.rhs)
+    eq_lhs = SymbolicUtils._iszero(expand_derivatives(eq.lhs)) ? eq.lhs :
+        expand_derivatives(eq.lhs)
+    eq_rhs = SymbolicUtils._iszero(expand_derivatives(eq.rhs)) ? eq.rhs :
+        expand_derivatives(eq.rhs)
     left_expr = transform_expression(pinnrep, toexpr(eq_lhs))
     right_expr = transform_expression(pinnrep, toexpr(eq_rhs))
     left_expr = _dot_(left_expr)
@@ -379,7 +403,7 @@ function get_vars(indvars_, depvars_)
     depvars = Symbol[]
     dict_depvar_input = Dict{Symbol, Vector{Symbol}}()
     for d in depvars_
-        if unwrap(d) isa SymbolicUtils.BasicSymbolic
+        if SymbolicUtils.iscall(unwrap(d))
             dname = SymbolicIndexingInterface.getname(d)
             push!(depvars, dname)
             push!(

src/training_strategies.jl

@@ -422,7 +422,7 @@ function get_loss_function(
         return solve(
             prob, strategy.quadrature_alg; strategy.reltol, strategy.abstol,
             strategy.maxiters
-        )[1]
+        ).u
     end
     return (θ) -> f_(lb, ub, loss_function, θ) / area
 end

src/transform_inf_integral.jl

@@ -43,7 +43,7 @@ function v_inf(t)
 end
 
 function v_semiinf(t, a, upto_inf)
-    if a isa Num
+    if a isa Num || a isa SymbolicUtils.BasicSymbolic
         if upto_inf == true
             return :($t ./ (1 .- $t))
         else
@@ -84,6 +84,16 @@ function transform_inf_integral(
     )
     lb_ = Symbolics.tosymbol.(lb)
     ub_ = Symbolics.tosymbol.(ub)
+    # Convert bounds to plain Julia numbers where possible to avoid Symbolics
+    # arithmetic in boolean mask operations (e.g., false * Num(Inf) gives NaN)
+    lb = map(
+        (l, ls) -> ls === -Inf ? -Inf : ls === Inf ? Inf : (l isa Num ? SymbolicUtils.unwrap(l) : l),
+        lb, lb_
+    )
+    ub = map(
+        (u, us) -> us === Inf ? Inf : us === -Inf ? -Inf : (u isa Num ? SymbolicUtils.unwrap(u) : u),
+        ub, ub_
+    )
 
     if -Inf in lb_ || Inf in ub_
         if !(integrating_variable isa Array)
@@ -118,10 +128,42 @@ function transform_inf_integral(
 
         ϵ = 1 / 20 #cbrt(eps(eltypeθ))
 
-        lb = 0.0 .* _semiup + (-1.0 + ϵ) .* _inf + (-1.0 + ϵ) .* _semilw + _none .* lb +
-            lb ./ (1 .+ lb) .* _num_semiup + (-1.0 + ϵ) .* _num_semilw
-        ub = (1.0 - ϵ) .* _semiup + (1.0 - ϵ) .* _inf + 0.0 .* _semilw + _none .* ub +
-            (1.0 - ϵ) .* _num_semiup + ub ./ (1 .+ ub) .* _num_semilw
+        # Use conditional logic instead of boolean mask multiplication to avoid
+        # NaN from 0.0 * Inf (IEEE754 behavior) with symbolic infinity bounds
+        lb = map(eachindex(lb)) do i
+            if _none[i]
+                lb[i]
+            elseif _inf[i]
+                -1.0 + ϵ
+            elseif _semiup[i]
+                0.0
+            elseif _semilw[i]
+                -1.0 + ϵ
+            elseif _num_semiup[i]
+                lb[i] / (1 + lb[i])
+            elseif _num_semilw[i]
+                -1.0 + ϵ
+            else
+                lb[i]
+            end
+        end
+        ub = map(eachindex(ub)) do i
+            if _none[i]
+                ub[i]
+            elseif _inf[i]
+                1.0 - ϵ
+            elseif _semiup[i]
+                1.0 - ϵ
+            elseif _semilw[i]
+                0.0
+            elseif _num_semiup[i]
+                1.0 - ϵ
+            elseif _num_semilw[i]
+                ub[i] / (1 + ub[i])
+            else
+                ub[i]
+            end
+        end
 
         j = get_inf_transformation_jacobian(
             integrating_var_transformation, _inf, _semiup,

test/BPINN_PDE_tests.jl

@@ -262,7 +262,7 @@ end
         [t],
         [u(t)],
         [p],
-        defaults = Dict([p => 4.0])
+        initial_conditions = Dict([p => 4.0])
     )
 
     analytic_sol_func1(u0, t) = u0 + sinpi(2t) / (2π)
@@ -355,7 +355,7 @@ end
 
     @named pde_system = PDESystem(
         eqs, bcs, domains,
-        [t], [x(t), y(t), z(t)], [σ_], defaults = Dict([p => 1.0 for p in [σ_]])
+        [t], [x(t), y(t), z(t)], [σ_], initial_conditions = Dict([p => 1.0 for p in [σ_]])
     )
 
     sol1 = ahmc_bayesian_pinn_pde(
@@ -491,7 +491,7 @@ end
         [x, t],
         [u(x, t)],
         [α],
-        defaults = Dict([α => 2.0])
+        initial_conditions = Dict([α => 2.0])
     )
 
     # neccesarry for loss function construction (involves Operator masking)

test/NNPDE_tests.jl

@@ -170,6 +170,8 @@ end
 @testitem "PDE II: 2D Poisson" tags = [:nnpde1] setup = [NNPDE1TestSetup] begin
     using Lux, Random, Optimisers, DomainSets, Cubature, QuasiMonteCarlo, Integrals
     import DomainSets: Interval, infimum, supremum
+    using OptimizationOptimJL: BFGS
+    using LineSearches: BackTracking
 
     function test_2d_poisson_equation(chain, strategy)
         @parameters x y
@@ -196,7 +198,9 @@ end
         discretization = PhysicsInformedNN(chain, strategy; init_params = ps)
         @named pde_system = PDESystem(eq, bcs, domains, [x, y], [u(x, y)])
         prob = discretize(pde_system, discretization)
-        res = solve(prob, Adam(0.1); maxiters = 500, callback)
+        res = solve(prob, Adam(0.01); maxiters = 1000, callback)
+        prob = remake(prob, u0 = res.u)
+        res = solve(prob, BFGS(linesearch = BackTracking()); maxiters = 1000)
         phi = discretization.phi
 
         xs, ys = [infimum(d.domain):0.01:supremum(d.domain) for d in domains]
@@ -404,9 +408,9 @@ end
     end
 
     # Adam warmup for robustness, then BFGS for convergence
-    res = solve(prob, Adam(0.01); maxiters = 500)
+    res = solve(prob, Adam(0.01); maxiters = 1000)
     prob = remake(prob, u0 = res.u)
-    res = solve(prob, BFGS(linesearch = BackTracking()); maxiters = 500)
+    res = solve(prob, BFGS(linesearch = BackTracking()); maxiters = 1000)
 
     dx = 0.1
     xs, ts = [infimum(d.domain):dx:supremum(d.domain) for d in domains]