substitute using SymbolicUtils

Fixes #339

MWE:

```julia
using ModelingToolkit
@parameters α β γ δ
s = α => γ
expr = (((1 / β - 1) + δ) / α) ^ (1 / (α - 1))
ModelingToolkit.substitute(expr, s)
new  = (((1 / β - 1) + δ) / γ) ^ (1 / (γ - 1))
isequal(expr, new)
```

Author: ChrisRackauckas

src/ModelingToolkit.jl

@@ -13,6 +13,9 @@ import GeneralizedGenerated
 using DocStringExtensions
 using Base: RefValue
 
+import SymbolicUtils
+import SymbolicUtils: to_symbolic, FnType, simplify
+
 import TreeViews
 
 """
@@ -112,17 +115,18 @@ export Reaction, ReactionSystem
 export Differential, expand_derivatives, @derivatives
 export IntervalDomain, ProductDomain, ⊗, CircleDomain
 export Equation, ConstrainedEquation
-export simplify_constants
-
 export Operation, Expression, Variable
+export independent_variable, states, parameters, equations
+
 export calculate_jacobian, generate_jacobian, generate_function
 export calculate_tgrad, generate_tgrad
 export calculate_gradient, generate_gradient
 export calculate_factorized_W, generate_factorized_W
 export calculate_hessian, generate_hessian
 export calculate_massmatrix, generate_diffusion_function
-export independent_variable, states, parameters, equations
-export simplified_expr, rename, get_variables, substitute_expr!
+
+export simplified_expr, rename, get_variables
+export simplify, substitute
 export build_function
 export @register
 export modelingtoolkitize

src/differentials.jl

@@ -49,10 +49,10 @@ function expand_derivatives(O::Operation)
 
         return sum(1:length(o.args)) do i
             derivative(o, i) * expand_derivatives(D(o.args[i]))
-        end |> simplify_constants
+        end |> simplify
     end
 
-    return simplify_constants(O)
+    return simplify(O)
 end
 expand_derivatives(x) = x
 

src/direct.jl

@@ -8,7 +8,7 @@ an array of variable expressions.
 """
 function gradient(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
     out = [expand_derivatives(Differential(v)(O)) for v in vars]
-    simplify ? simplify_constants.(out) : out
+    simplify ? simplify.(out) : out
 end
 
 """
@@ -21,7 +21,7 @@ an array of variable expressions.
 """
 function jacobian(ops::AbstractVector{<:Expression}, vars::AbstractVector{<:Expression}; simplify = true)
     out = [expand_derivatives(Differential(v)(O)) for O in ops, v in vars]
-    simplify ? simplify_constants.(out) : out
+    simplify ? simplify.(out) : out
 end
 
 """
@@ -34,7 +34,7 @@ an array of variable expressions.
 """
 function hessian(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
     out = [expand_derivatives(Differential(v2)(Differential(v1)(O))) for v1 in vars, v2 in vars]
-    simplify ? simplify_constants.(out) : out
+    simplify ? simplify.(out) : out
 end
 
 function simplified_expr(O::Operation)

src/simplify.jl

@@ -1,6 +1,3 @@
-import SymbolicUtils
-import SymbolicUtils: FnType
-
 # ModelingToolkit -> SymbolicUtils
 SymbolicUtils.istree(x::Operation) = true
 function SymbolicUtils.operation(x::Operation)
@@ -38,10 +35,12 @@ SymbolicUtils.symtype(x::Expression) = Number
 
 # SymbolicUtils -> ModelingToolkit
 
-function simplify_constants(expr)
+function SymbolicUtils.simplify(expr::Expression)
     SymbolicUtils.simplify(expr) |> to_mtk
 end
 
+@deprecate simplify_constants(ex) simplify(ex)
+
 to_mtk(x) = x
 to_mtk(x::Number) = Constant(x)
 to_mtk(v::SymbolicUtils.Sym{T}) where {T} = Variable{T}(nameof(v))

src/systems/diffeqs/abstractodesystem.jl

@@ -69,13 +69,13 @@ function calculate_factorized_W(sys::AbstractODESystem, simplify=true)
     Wfact = lu(W, Val(false), check=false).factors
 
     if simplify
-        Wfact = simplify_constants.(Wfact)
+        Wfact = simplify.(Wfact)
     end
 
     W_t = - LinearAlgebra.I/gam + jac
     Wfact_t = lu(W_t, Val(false), check=false).factors
     if simplify
-        Wfact_t = simplify_constants.(Wfact_t)
+        Wfact_t = simplify.(Wfact_t)
     end
     sys.Wfact[] = Wfact
     sys.Wfact_t[] = Wfact_t
@@ -113,7 +113,7 @@ function calculate_massmatrix(sys::AbstractODESystem, simplify=true)
             error("Only semi-explicit constant mass matrices are currently supported")
         end
     end
-    M = simplify ? simplify_constants.(M) : M
+    M = simplify ? simplify.(M) : M
     # M should only contain concrete numbers
     M = map(x->x isa Constant ? x.value : x, M)
     M == I ? I : M

src/utils.jl

@@ -100,17 +100,32 @@ end
 
 # variable substitution
 """
-substitute_expr!(expr::Operation, s::Pair{Operation, Operation})
+substitute(expr::Operation, s::Pair{Operation, Operation})
+substitute(expr::Operation, s::Dict{Operation, Operation})
+substitute(expr::Operation, s::Vector{Pair{Operation, Operation}})
 
 Performs the substitution `Operation => Operation` on the `expr` Operation.
 """
-substitute_expr!(expr::Constant, s::Pair{Operation, Operation}) = nothing
-function substitute_expr!(expr::Operation, s::Pair{Operation, Operation})
-    if !is_singleton(expr)
-        expr.args .= replace(expr.args, s)
-        for arg in expr.args
-            substitute_expr!(arg, s)
-        end
-    end
-    return nothing
+substitute(expr::Constant, s) = nothing
+
+function substitute(expr::Operation, s::Pair{Operation, Operation})
+    substitute(expr,Dict(to_symbolic(s[1]) => to_symbolic(s[2])))
+end
+
+function substitute(expr::Operation, dict::Dict{Operation,Operation})
+    substitute(expr,Dict(map(Pair,
+                         map(to_symbolic, keys(dict)),
+                         map(to_symbolic, values(dict)))))
 end
+
+function substitute(expr::Operation, s::Vector)
+    substitute(expr,Dict(map(Pair,
+                         map(to_symbolic, first.(s)),
+                         map(to_symbolic, last.(s)))))
+end
+
+function substitute(expr::Operation, s::Dict{<:SymbolicUtils.Term,<:SymbolicUtils.Term})
+    SymbolicUtils.RuleSet([SymbolicUtils.@rule(~x::(x->haskey(s, x)) => s[~x])])(expr) |> to_mtk
+end
+
+@deprecate substitute_expr!(expr,s) substitute(expr,s)

test/variable_utils.jl

@@ -7,6 +7,7 @@ sol = ModelingToolkit.get_variables(expr)
 
 @parameters γ
 s = α => γ
+expr = (((1 / β - 1) + δ) / α) ^ (1 / (α - 1))
+ModelingToolkit.substitute(expr, s)
 new = (((1 / β - 1) + δ) / γ) ^ (1 / (γ - 1))
-ModelingToolkit.substitute_expr!(expr, s)
 @test isequal(expr, new)