From 62a4a3c4535a6d8d472d81e4f529fc4f14293f2d Mon Sep 17 00:00:00 2001 From: Alvaro de Diego Date: Mon, 20 Apr 2020 11:01:31 +0200 Subject: [PATCH 1/4] now using ModelingToolkit for Differentiation --- src/CoherentStructures.jl | 3 +++ src/streammacros.jl | 20 ++++++++++++++++++++ src/velocityfields.jl | 1 + 3 files changed, 24 insertions(+) diff --git a/src/CoherentStructures.jl b/src/CoherentStructures.jl index ef0fae44..c877b89e 100644 --- a/src/CoherentStructures.jl +++ b/src/CoherentStructures.jl @@ -42,6 +42,9 @@ import JuAFEM const JFM = JuAFEM using RecipesBase import SymEngine +import ModelingToolkit +using ModelingToolkit: Variable, Differential, simplified_expr, + expand_derivatives, Expression, Operation, Constant # contains a list of exported functions include("exports.jl") diff --git a/src/streammacros.jl b/src/streammacros.jl index 84551493..58b84c61 100644 --- a/src/streammacros.jl +++ b/src/streammacros.jl @@ -219,6 +219,8 @@ end # symbolic differentiation of expressions using SymEngine # ######################################################################################## +ModelingToolkit.derivative(::typeof(sign), x, ::Val{1}) = 0 + sgn(x) = (x > 0) ? one(x) : (x < 0) ? -one(x) : zero(x) heaviside(x) = 0 < x ? one(x) : zero(x) window(x, a, b) = (a <= x < b) ? one(x) : zero(x) @@ -232,6 +234,13 @@ diff_dict[:zero] = :zero diff_dict[:window] = :zero function expr_diff(expr::Expr, var::Symbol) + D = Differential(var) + + toolkit_expr = convert(Expression, expr) + d_expr = simplified_expr(expand_derivatives(D(toolkit_expr))) + d_expr = convert(Expression, simple_simplifier(d_expr)) + return simplified_expr(propagate_constants(d_expr)) + # not a nice way to differentiate expressions, but ReverseDiffSource # is broken. expr_sym = SymEngine.Basic(expr) @@ -300,6 +309,10 @@ function simple_simplifier(expr::Expr) if expr.head !== :call return Expr(expr.head, args...) end + if args[1] === :one + return 1 + end + if args[1] === :zero return 0 end @@ -334,6 +347,13 @@ function simple_simplifier(expr::Expr) return Expr(expr.head, args...) end +function propagate_constants(expr::ModelingToolkit.Operation) + return expr.op(propagate_constants.(expr.args)...) +end +propagate_constants(x::Constant) = x.value +propagate_constants(x) = x + + simple_simplifier(expr) = expr expr_grad(expr, coord_vars::Vector{Symbol}) = expr_diff.(expr, coord_vars) diff --git a/src/velocityfields.jl b/src/velocityfields.jl index db96e5d1..e05bafbd 100644 --- a/src/velocityfields.jl +++ b/src/velocityfields.jl @@ -26,6 +26,7 @@ bickleyJetEqVari! = ODE.ODEFunction{true}((DU, U, p, t) -> DU .= bickleyJetEqVar # rotating double gyre flow [Mosovsky & Meiss, 2011] @define_stream Ψ_rot_dgyre begin st = heaviside(t)*heaviside(1-t)*t^2*(3-2*t) + heaviside(t-1) + heaviside(x)= 0.5*(sign(x) + 1) Ψ_P = sin(2π*x)*sin(π*y) Ψ_F = sin(π*x)*sin(2π*y) Ψ_rot_dgyre = (1-st) * Ψ_P + st * Ψ_F From eab27968af5472ec611d9967e9a054babafe2606 Mon Sep 17 00:00:00 2001 From: Alvaro de Diego Date: Mon, 20 Apr 2020 11:14:46 +0200 Subject: [PATCH 2/4] fixed error with sign function --- src/streammacros.jl | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/streammacros.jl b/src/streammacros.jl index 58b84c61..ea3765f9 100644 --- a/src/streammacros.jl +++ b/src/streammacros.jl @@ -346,17 +346,25 @@ function simple_simplifier(expr::Expr) end return Expr(expr.head, args...) end +simple_simplifier(expr) = expr + function propagate_constants(expr::ModelingToolkit.Operation) - return expr.op(propagate_constants.(expr.args)...) + if all(isconstant.(expr.args)) + return expr.op(propagate_constants.(expr.args)...) + else + return Operation(expr.op, propagate_constants.(expr.args)) + end end propagate_constants(x::Constant) = x.value propagate_constants(x) = x +isconstant(x::Operation) = !(x.op isa Variable) && all(isconstant.(x.args)) +isconstant(x::Constant) = true + -simple_simplifier(expr) = expr -expr_grad(expr, coord_vars::Vector{Symbol}) = expr_diff.(expr, coord_vars) +expr_grad(expr, coord_vars::Vector{Symbol}) = expr_diff.(expr, coord_vars) function hessian(expr, coord_vars) ∇expr = expr_grad(expr, coord_vars) ∇²expr = expr_grad(expr) From 526976f2baf89bd41ea6c32baa931b7acf8c54a0 Mon Sep 17 00:00:00 2001 From: Alvaro de Diego Date: Mon, 20 Apr 2020 11:51:09 +0200 Subject: [PATCH 3/4] clean up parts of the code responsible for SymEngine --- src/CoherentStructures.jl | 1 - src/streammacros.jl | 78 +-------------------------------------- 2 files changed, 2 insertions(+), 77 deletions(-) diff --git a/src/CoherentStructures.jl b/src/CoherentStructures.jl index c877b89e..97e153a2 100644 --- a/src/CoherentStructures.jl +++ b/src/CoherentStructures.jl @@ -41,7 +41,6 @@ const VD = VoronoiDelaunay import JuAFEM const JFM = JuAFEM using RecipesBase -import SymEngine import ModelingToolkit using ModelingToolkit: Variable, Differential, simplified_expr, expand_derivatives, Expression, Operation, Constant diff --git a/src/streammacros.jl b/src/streammacros.jl index ea3765f9..cb154335 100644 --- a/src/streammacros.jl +++ b/src/streammacros.jl @@ -196,7 +196,7 @@ end function streamline_derivatives(H::Symbol, formulas::Expr) # symbols that are not supposed to be substituted # (additional to symbols defined in Base) - bound_symbols = [:x, :y, :p, :t, keys(diff_dict)...] + bound_symbols = [:x, :y, :p, :t] H = substitutions(H, formulas, bound_symbols) # symbolic gradient and hessian (note the broadcast) @@ -216,23 +216,11 @@ function streamline_derivatives(H::Symbol, formulas::Expr) end ######################################################################################## -# symbolic differentiation of expressions using SymEngine # +# symbolic differentiation of expressions using ModelingToolkit # ######################################################################################## ModelingToolkit.derivative(::typeof(sign), x, ::Val{1}) = 0 -sgn(x) = (x > 0) ? one(x) : (x < 0) ? -one(x) : zero(x) -heaviside(x) = 0 < x ? one(x) : zero(x) -window(x, a, b) = (a <= x < b) ? one(x) : zero(x) - -# manually define derivatives for functions that SymEngine cant differentiate -diff_dict = Dict() -diff_dict[:abs] = :sgn -diff_dict[:sgn] = :zero -diff_dict[:heaviside] = :zero -diff_dict[:zero] = :zero -diff_dict[:window] = :zero - function expr_diff(expr::Expr, var::Symbol) D = Differential(var) @@ -241,68 +229,9 @@ function expr_diff(expr::Expr, var::Symbol) d_expr = convert(Expression, simple_simplifier(d_expr)) return simplified_expr(propagate_constants(d_expr)) - # not a nice way to differentiate expressions, but ReverseDiffSource - # is broken. - expr_sym = SymEngine.Basic(expr) - d_expr_sym = SymEngine.diff(expr_sym, var) - d_expr = Meta.parse(SymEngine.toString.(d_expr_sym)) - - # resolve derivatives that SymEngine doesn't know using diff_dict - d_expr = additional_derivatives(d_expr) - - # clean up unresolved substitutions that result from SymEnginge treating - # unknown derivatives - d_expr = substitution_cleaner(d_expr) - - # clean up zeros - d_expr = simple_simplifier(d_expr) end expr_diff(expr, var::Symbol) = expr == var ? 1 : 0 -function additional_derivatives(expr::Expr) - # some functions like abs(x) are not treated by SymEngine and block the - # expression. For example, diff(Basic(:(abs(x^2+1)),:x)) returns a SymEngine Object - # whose string representation is parsed to :(Derivative(abs(1 + x ^ 2), x)), - # whose AST is: - # Expr - # head: Symbol call - # args: Array{Any}((3,)) - # 1: Symbol Derivative - # 2: Expr - # ... Body of expression ... - # 3: Symbol x - # typ: Any - - # detect expressions of this form - if expr.head === :call && expr.args[1] === :Derivative - f = expr.args[2].args[1] - var = expr.args[3] - f_arg = expr.args[2].args[2] - if haskey(diff_dict, f) # try if diff_dict provides a rescue - df = diff_dict[f] - inner_d = expr_diff(f_arg, var) - df_computed_manually = :($df($f_arg) * $inner_d) - return additional_derivatives(df_computed_manually) # handle nested problems - end - end - return Expr(expr.head, additional_derivatives.(expr.args)...) # call recursively on subexpressions -end -additional_derivatives(expr) = expr - -# A second thing that SymEngine does is returning expressions of the form -# Subs(ex1, symb1, ex2). Resolve these substitutions -function substitution_cleaner(expr::Expr) - if expr.head === :call && expr.args[1] === :Subs - return substitution_cleaner(sym_subst( - expr.args[2], - expr.args[3], - expr.args[4], - )) - end - return Expr(expr.head, substitution_cleaner.(expr.args)...) -end -substitution_cleaner(expr) = expr - # perform some basic simplifications like getting rid of ones and zeros function simple_simplifier(expr::Expr) args = simple_simplifier.(expr.args) @@ -348,7 +277,6 @@ function simple_simplifier(expr::Expr) end simple_simplifier(expr) = expr - function propagate_constants(expr::ModelingToolkit.Operation) if all(isconstant.(expr.args)) return expr.op(propagate_constants.(expr.args)...) @@ -362,8 +290,6 @@ propagate_constants(x) = x isconstant(x::Operation) = !(x.op isa Variable) && all(isconstant.(x.args)) isconstant(x::Constant) = true - - expr_grad(expr, coord_vars::Vector{Symbol}) = expr_diff.(expr, coord_vars) function hessian(expr, coord_vars) ∇expr = expr_grad(expr, coord_vars) From d0934cbd4ac1d4ada29d2dd4be7a8c4265ea0c48 Mon Sep 17 00:00:00 2001 From: Alvaro de Diego Date: Mon, 20 Apr 2020 12:46:05 +0200 Subject: [PATCH 4/4] added ModelingToolkit as a dependency --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 0e9ff2f1..1d174520 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,7 @@ IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" JuAFEM = "30d91d44-8115-11e8-1d28-c19a5ac16de8" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" @@ -22,7 +23,6 @@ SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" VoronoiDelaunay = "72f80fcb-8c52-57d9-aff0-40c1a3526986" @@ -35,12 +35,12 @@ GeometricalPredicates = "0.3, 0.4" Interpolations = "0.11, 0.12" IterativeSolvers = "0.7, 0.8" LinearMaps = "2" +ModelingToolkit = "3.0.1" NearestNeighbors = "0.4.4" OrdinaryDiffEq = "5.18" ProgressMeter = "1" RecipesBase = "0.7, 0.8, 1.0" StaticArrays = "0.10, 0.11, 0.12" -SymEngine = "0.6, 0.7, 0.8" Tensors = "1.0.1" VoronoiDelaunay = "0.4" julia = "1"