feat: move fast_substitute to Symbolics, implement SII.symbolic_evaluate

Project.toml

@@ -81,7 +81,7 @@ SciMLBase = "2"
 Setfield = "1"
 SpecialFunctions = "2"
 StaticArrays = "1.1"
-SymbolicIndexingInterface = "0.3"
+SymbolicIndexingInterface = "0.3.12"
 SymbolicLimits = "0.2.0"
 SymbolicUtils = "1.4"
 julia = "1.10"

src/variable.jl

@@ -410,6 +410,98 @@ end
 
 SymbolicIndexingInterface.getname(x, val=_fail) = _getname(unwrap(x), val)
 
+function SymbolicIndexingInterface.symbolic_evaluate(ex::Union{Num, Arr, Symbolic, Equation, Inequality}, d::Dict; kwargs...)
+    fixpoint_sub(ex, d; kwargs...)
+end
+
+function fixpoint_sub(x, dict; operator = Nothing)
+    y = fast_substitute(x, dict; operator)
+    while !isequal(x, y)
+        y = x
+        x = fast_substitute(y, dict; operator)
+    end
+
+    return x
+end
+
+const Eq = Union{Equation, Inequality}
+function fast_substitute(eq::Eq, subs; operator = Nothing)
+    if eq isa Inequality
+        Inequality(fast_substitute(eq.lhs, subs; operator),
+            fast_substitute(eq.rhs, subs; operator),
+            eq.relational_op)
+    else
+        Equation(fast_substitute(eq.lhs, subs; operator),
+            fast_substitute(eq.rhs, subs; operator))
+    end
+end
+function fast_substitute(eq::T, subs::Pair; operator = Nothing) where {T <: Eq}
+    T(fast_substitute(eq.lhs, subs; operator), fast_substitute(eq.rhs, subs; operator))
+end
+function fast_substitute(eqs::AbstractArray, subs; operator = Nothing)
+    fast_substitute.(eqs, (subs,); operator)
+end
+function fast_substitute(eqs::AbstractArray, subs::Pair; operator = Nothing)
+    fast_substitute.(eqs, (subs,); operator)
+end
+for (exprType, subsType) in Iterators.product((Num, Symbolics.Arr), (Any, Pair))
+    @eval function fast_substitute(expr::$exprType, subs::$subsType; operator = Nothing)
+        fast_substitute(value(expr), subs; operator)
+    end
+end
+function fast_substitute(expr, subs; operator = Nothing)
+    if (_val = get(subs, expr, nothing)) !== nothing
+        return _val
+    end
+    istree(expr) || return expr
+    op = fast_substitute(operation(expr), subs; operator)
+    args = SymbolicUtils.unsorted_arguments(expr)
+    if !(op isa operator)
+        canfold = Ref(!(op isa Symbolic))
+        args = let canfold = canfold
+            map(args) do x
+                x′ = fast_substitute(x, subs; operator)
+                canfold[] = canfold[] && !(x′ isa Symbolic)
+                x′
+            end
+        end
+        canfold[] && return op(args...)
+    end
+    similarterm(expr,
+        op,
+        args,
+        symtype(expr);
+        metadata = metadata(expr))
+end
+function fast_substitute(expr, pair::Pair; operator = Nothing)
+    a, b = pair
+    isequal(expr, a) && return b
+    if a isa AbstractArray
+        for (ai, bi) in zip(a, b)
+            expr = fast_substitute(expr, ai => bi; operator)
+        end
+    end
+    istree(expr) || return expr
+    op = fast_substitute(operation(expr), pair; operator)
+    args = SymbolicUtils.unsorted_arguments(expr)
+    if !(op isa operator)
+        canfold = Ref(!(op isa Symbolic))
+        args = let canfold = canfold
+            map(args) do x
+                x′ = fast_substitute(x, pair; operator)
+                canfold[] = canfold[] && !(x′ isa Symbolic)
+                x′
+            end
+        end
+        canfold[] && return op(args...)
+    end
+    similarterm(expr,
+        op,
+        args,
+        symtype(expr);
+        metadata = metadata(expr))
+end
+
 function getparent(x, val=_fail)
     maybe_parent = getmetadata(x, Symbolics.GetindexParent, nothing)
     if maybe_parent !== nothing

test/runtests.jl

@@ -63,6 +63,9 @@ if GROUP == "All" || GROUP == "Core" || GROUP == "SymbolicIndexingInterface"
     @safetestset "SymbolicIndexingInterface Parameter Indexing Test" begin
         include("symbolic_indexing_interface_parameter_indexing.jl")
     end
+    @safetestset "SymbolicIndexingInterface Symbolic Evaluate Test" begin
+        include("symbolic_indexing_interface_symbolic_evaluate.jl")
+    end
 end
 
 if GROUP == "Downstream"

test/symbolic_indexing_interface_symbolic_evaluate.jl

@@ -0,0 +1,36 @@
+using Symbolics
+using SymbolicIndexingInterface
+using Symbolics: Differential, Operator
+
+@variables t x(t) y(t)
+@variables p[1:3, 1:3] q[1:3]
+
+bar(x, p) = p * x
+@register_array_symbolic bar(x::AbstractVector, p::AbstractMatrix) begin
+    size = size(x)
+    eltype = promote_type(eltype(x), eltype(p))
+end
+
+D = Differential(t)
+
+expr1 = x + y + D(x)
+@test isequal(symbolic_evaluate(expr1, Dict(x => 3)), 3 + y + D(3))
+@test isequal(symbolic_evaluate(expr1, Dict(x => 3); operator = Operator), 3 + y + D(x))
+@test isequal(symbolic_evaluate(expr1, Dict(x => 1, D(x) => 2)), y + 3)
+@test symbolic_evaluate(expr1, Dict(x => 1, D(x) => 2, y => 3)) == 6
+@test isequal(symbolic_evaluate(expr1, Dict(x => 3, y => 3x), operator = Operator), 12 + D(x))
+@test symbolic_evaluate(expr1, Dict(x => 3, y => 3x, D(x) => 2)) == 14
+
+expr2 = bar(q, p)
+@test isequal(symbolic_evaluate(expr2, Dict(p => ones(3, 3))), bar(q, ones(3, 3)))
+@test symbolic_evaluate(expr2, Dict(p => ones(3, 3), q => ones(3))) == 3ones(3)
+
+expr3 = bar(3q, 3p)
+@test isequal(symbolic_evaluate(expr3, Dict(p => ones(3, 3))), bar(3q, 3ones(3, 3)))
+@test symbolic_evaluate(expr3, Dict(p => ones(3, 3), q => ones(3))) == 27ones(3)
+
+expr4 = D(x) ~ 3x + y
+@test isequal(symbolic_evaluate(expr4, Dict(x => 3)), D(3) ~ 9 + y)
+@test isequal(symbolic_evaluate(expr4, Dict(x => 3); operator = Operator), D(x) ~ y + 9)
+@test isequal(symbolic_evaluate(expr4, Dict(x => 1, D(x) => 2)), 2 ~ 3 + y)
+@test isequal(symbolic_evaluate(expr4, Dict(x => 1, D(x) => 2, y => 3)), 2 ~ 6)

test/symbolic_indexing_interface_trait.jl

@@ -10,3 +10,11 @@ using SymbolicIndexingInterface
 @variables y[1:3]
 @test symbolic_type(y) == ArraySymbolic()
 @test all(symbolic_type.(collect(y)) .== (ScalarSymbolic(),))
+
+@variables x y z
+subs = Dict(x => 0.1, y => 2z)
+subs2 = merge(subs, Dict(z => 2x+3))
+
+@test symbolic_evaluate(x, subs) == 0.1
+@test isequal(symbolic_evaluate(y, subs), 2z)
+@test symbolic_evaluate(y, subs2) == 6.4