Add conversions from MT to Hyperplane and HalfSpace (#2193)

* Add conversions from MT too Hyperplane and HalfSpace * test MT * add docstrings * Update Project.toml * Update src/Sets/HalfSpace.jl Co-authored-by: Christian Schilling <[email protected]> * Update src/Sets/Hyperplane.jl Co-authored-by: Christian Schilling <[email protected]> * Update Project.toml Co-authored-by: Christian Schilling <[email protected]> * Update src/Sets/Hyperplane.jl Co-authored-by: Christian Schilling <[email protected]> * Update src/Sets/HalfSpace.jl * Update src/Sets/Hyperplane.jl Co-authored-by: Christian Schilling <[email protected]>
JuliaReach · Jul 17, 2020 · 4dee6a2 · 4dee6a2
1 parent ed100b0
commit 4dee6a2
Show file tree

Hide file tree

Showing 9 changed files with 256 additions and 1 deletion.
diff --git a/Project.toml b/Project.toml
@@ -32,6 +32,7 @@ JuMP = "0.21"
 Makie = "0.9"
 MathProgBase = "0.7"
 MiniQhull = "0.1"
+ModelingToolkit = "1.1.3, 3.10"
 Optim = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21"
 Plots = "0.25, 0.26, 0.27, 0.28, 0.29"
 Polyhedra = "0.6"
@@ -50,11 +51,13 @@ GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MiniQhull = "978d7f02-9e05-4691-894f-ae31a51d76ca"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
 TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["CDDLib", "Distributions", "Documenter", "Expokit", "GR", "IntervalMatrices", "Makie", "Optim", "Plots", "Polyhedra", "TaylorModels", "Test"]
+test = ["CDDLib", "Distributions", "Documenter", "Expokit", "GR", "IntervalMatrices",
+        "Makie", "ModelingToolkit", "Optim", "Plots", "Polyhedra", "TaylorModels", "Test"]
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -7,6 +7,7 @@ Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"

diff --git a/src/Initialization/init_ModelingToolkit.jl b/src/Initialization/init_ModelingToolkit.jl
@@ -0,0 +1,11 @@
+eval(quote
+    using .ModelingToolkit: get_variables,
+                            gradient,
+                            simplify,
+                            to_symbolic,
+                            Operation
+
+end)
+
+eval(load_modeling_toolkit_hyperplane())
+eval(load_modeling_toolkit_halfspace())
diff --git a/src/Sets/HalfSpace.jl b/src/Sets/HalfSpace.jl
@@ -519,3 +519,119 @@ end
 # TODO: after #2032, #2041 remove use of this function
 _normal_Vector(P::LazySet) = [LinearConstraint(convert(Vector, c.a), c.b) for c in constraints_list(P)]
 _normal_Vector(c::LinearConstraint) = LinearConstraint(convert(Vector, c.a), c.b)
+
+
+# ============================================
+# Functionality that requires ModelingToolkit
+# ============================================
+function load_modeling_toolkit_halfspace()
+return quote
+
+"""
+    HalfSpace(expr::Operation, vars=get_variables(expr); N::Type{<:Real}=Float64)
+
+Return the half-space given by a symbolic expression.
+
+### Input
+
+- `expr` -- symbolic expression that describes a half-space
+- `vars` -- (optional, default: `get_variables(expr)`), if an array of variables is given,
+            use those as the ambient variables in the set with respect to which derivations
+            take place; otherwise, use only the variables which appear in the given
+            expression (but be careful because the order may change; see the examples below for details)
+- `N`    -- (optional, default: `Float64`) the numeric type of the returned half-space
+
+### Output
+
+A `HalfSpace`.
+
+### Examples
+
+```julia
+julia> using ModelingToolkit
+
+julia> vars = @variables x y
+(x, y)
+
+julia> HalfSpace(x - y <= 2)
+HalfSpace{Float64,Array{Float64,1}}([1.0, -1.0], 2.0)
+
+julia> HalfSpace(x >= y)
+HalfSpace{Float64,Array{Float64,1}}([1.0, -1.0], -0.0)
+
+julia> vars = @variables x[1:4]
+(Operation[x₁, x₂, x₃, x₄],)
+
+julia> HalfSpace(x[1] >= x[2], x)
+HalfSpace{Float64,Array{Float64,1}}([-1.0, 1.0, 0.0, 0.0], -0.0)
+```
+
+Be careful with using the default `vars` value because it may introduce a wrong
+order.
+
+```julia
+julia> vars = @variables x y
+(x, y)
+
+julia> HalfSpace(2x ≥ 5y - 1) # wrong
+HalfSpace{Float64,Array{Float64,1}}([5.0, -2.0], 1.0)
+
+julia> HalfSpace(2x ≥ 5y - 1, vars) # correct
+HalfSpace{Float64,Array{Float64,1}}([-2.0, 5.0], 1.0)
+```
+
+### Algorithm
+
+It is assumed that the expression is of the form
+`EXPR0: α*x1 + ⋯ + α*xn + γ CMP β*x1 + ⋯ + β*xn + δ`,
+where `CMP` is one among `<`, `<=`, `≤`, `>`, `>=` or `≥`.
+This expression is transformed, by rearrangement and substitution, into the
+canonical form `EXPR1 : a1 * x1 + ⋯ + an * xn ≤ b`. The method used to identify
+the coefficients is to take derivatives with respect to the ambient variables `vars`.
+Therefore, the order in which the variables appear in `vars` affects the final result.
+Note in particular that strict inequalities are relaxed as being smaller-or-equal.
+Finally, the returned set is the half-space with normal vector `[a1, …, an]` and
+displacement `b`.
+"""
+function HalfSpace(expr::Operation, vars=get_variables(expr); N::Type{<:Real}=Float64)
+
+    # find sense and normalize
+    if expr.op == <
+        a, b = expr.args
+        sexpr = simplify(a - b)
+
+    elseif expr.op == >
+        a, b = expr.args
+        sexpr = simplify(b - a)
+
+    elseif (expr.op == |) && (expr.args[1].op == <)
+        a, b = expr.args[1].args
+        sexpr = simplify(a - b)
+
+    elseif (expr.op == |) && (expr.args[2].op == <)
+        a, b = expr.args[2].args
+        sexpr = simplify(a - b)
+
+    elseif (expr.op == |) && (expr.args[1].op == >)
+        a, b = expr.args[1].args
+        sexpr = simplify(b - a)
+
+    elseif (expr.op == |) && (expr.args[2].op == >)
+        a, b = expr.args[2].args
+        sexpr = simplify(b - a)
+
+    else
+        throw(ArgumentError("expected an expression describing a half-space, got $expr"))
+    end
+
+    # compute the linear coefficients by taking first order derivatives
+    coeffs = [N(α.value) for α in gradient(sexpr, collect(vars))]
+
+    # get the constant term by expression substitution
+    dvars = Dict(to_symbolic(vi) => zero(N) for vi in vars)
+    β = -N(ModelingToolkit.SymbolicUtils.substitute(to_symbolic(sexpr), dvars, fold=true))
+
+    return HalfSpace(coeffs, β)
+end
+
+end end  # quote / load_modeling_toolkit_halfspace()
diff --git a/src/Sets/Hyperplane.jl b/src/Sets/Hyperplane.jl
@@ -478,3 +478,78 @@ function translate(hp::Hyperplane{N}, v::AbstractVector{N}; share::Bool=false
     b = hp.b + dot(hp.a, v)
     return Hyperplane(a, b)
 end
+
+# ============================================
+# Functionality that requires ModelingToolkit
+# ============================================
+function load_modeling_toolkit_hyperplane()
+return quote
+
+"""
+    Hyperplane(expr::Operation, vars=get_variables(expr); N::Type{<:Real}=Float64)
+
+Return the hyperplane given by a symbolic expression.
+
+### Input
+
+- `expr` -- symbolic expression that describes a hyperplane
+- `vars` -- (optional, default: `get_variables(expr)`), if an array of variables is given,
+            use those as the ambient variables in the set with respect to which derivations
+            take place; otherwise, use only the variables which appear in the given
+            expression (but be careful because the order may change; in doubt, pass `vars` explicitly)
+- `N`    -- (optional, default: `Float64`) the numeric type of the returned hyperplane
+
+### Output
+
+A `Hyperplane`.
+
+### Examples
+
+```julia
+julia> using ModelingToolkit
+
+julia> vars = @variables x y
+(x, y)
+
+julia> Hyperplane(x - y == 2)
+Hyperplane{Float64,Array{Float64,1}}([1.0, -1.0], 2.0)
+
+julia> Hyperplane(x == y)
+Hyperplane{Float64,Array{Float64,1}}([1.0, -1.0], -0.0)
+
+julia> vars = @variables x[1:4]
+(Operation[x₁, x₂, x₃, x₄],)
+
+julia> Hyperplane(x[1] == x[2], x)
+Hyperplane{Float64,Array{Float64,1}}([1.0, -1.0, 0.0, 0.0], -0.0)
+```
+
+### Algorithm
+
+It is assumed that the expression is of the form
+`EXPR0: α*x1 + ⋯ + α*xn + γ == β*x1 + ⋯ + β*xn + δ`.
+This expression is transformed, by rearrangement and substitution, into the
+canonical form `EXPR1 : a1 * x1 + ⋯ + an * xn == b`. The method used to identify
+the coefficients is to take derivatives with respect to the ambient variables `vars`.
+Therefore, the order in which the variables appear in `vars` affects the final result.
+Finally, the returned set is the hyperplane with normal vector `[a1, …, an]` and
+displacement `b`.
+"""
+function Hyperplane(expr::Operation, vars=get_variables(expr); N::Type{<:Real}=Float64)
+    (expr.op == ==) || throw(ArgumentError("expected an expression of the form `ax == b`, got $expr"))
+
+    # simplify to the form a*x + β == 0
+    a, b = expr.args
+    sexpr = simplify(a - b)
+
+    # compute the linear coefficients by taking first order derivatives
+    coeffs = [N(α.value) for α in gradient(sexpr, collect(vars))]
+
+    # get the constant term by expression substitution
+    dvars = Dict(to_symbolic(vi) => zero(N) for vi in vars)
+    β = -N(ModelingToolkit.SymbolicUtils.substitute(to_symbolic(sexpr), dvars, fold=true))
+
+    return Hyperplane(coeffs, β)
+end
+
+end end  # quote / load_modeling_toolkit_hyperplane()
diff --git a/src/init.jl b/src/init.jl
@@ -4,6 +4,7 @@ function __init__()
     @require Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4" include("Initialization/init_Expokit.jl")
     @require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("Initialization/init_Makie.jl")
     @require MiniQhull = "978d7f02-9e05-4691-894f-ae31a51d76ca" include("Initialization/init_MiniQhull.jl")
+    @require ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" include("Initialization/init_ModelingToolkit.jl")
     @require Optim = "429524aa-4258-5aef-a3af-852621145aeb" include("Initialization/init_Optim.jl")
     @require Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" include("Initialization/init_Polyhedra.jl")
 end

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -10,6 +10,9 @@ using IntervalArithmetic: IntervalBox
 import Distributions, Expokit, IntervalMatrices, Optim, TaylorModels
 using IntervalMatrices: ±, IntervalMatrix
 using TaylorModels: set_variables, TaylorModelN
+@static if VERSION >= v"1.3"
+    using ModelingToolkit
+end
 
 # ==============================
 # Non-exported helper functions

diff --git a/test/unit_HalfSpace.jl b/test/unit_HalfSpace.jl
@@ -132,3 +132,34 @@ for N in [Float64, Float32]
     @test normalize(hs1, N(1)) == HalfSpace(N[1//3, 2//3], N(1))
     @test normalize(hs1, N(Inf)) == HalfSpace(N[1//2, 1], N(3//2))
 end
+
+# tests that require ModelingToolkit
+@static if VERSION >= v"1.3" && isdefined(@__MODULE__, :ModelingToolkit)
+
+    vars = @variables x y
+    @test HalfSpace(2x + 3y < 5) == HalfSpace([2.0, 3.0], 5.0)
+    @test HalfSpace(2x + 3y < 5, vars) == HalfSpace([2.0, 3.0], 5.0)
+    @test HalfSpace(2x + 3y < 5, N=Int) == HalfSpace([2, 3], 5)
+
+    @test HalfSpace(2x + 3y > 5) == HalfSpace([-2.0, -3.0], -5.0)
+    @test HalfSpace(2x + 3y > 5, vars) == HalfSpace([-2.0, -3.0], -5.0)
+
+    @test HalfSpace(2x + 3y ≤ 5) == HalfSpace([2.0, 3.0], 5.0)
+    @test HalfSpace(2x + 3y ≤ 5, vars) == HalfSpace([2.0, 3.0], 5.0)
+
+    @test HalfSpace(2x <= 5y - 1) == HalfSpace([2.0, -5.0], -1.0)
+    @test HalfSpace(2x ≤ 5y - 1) == HalfSpace([2.0, -5.0], -1.0)
+
+    @test HalfSpace(2x + 3y ≥ 5) == HalfSpace([-2.0, -3.0], -5.0)
+    @test HalfSpace(2x + 3y ≥ 5, vars) == HalfSpace([-2.0, -3.0], -5.0)
+
+    # doesn't work because get_vars returns variables [y, x]
+    # => both tests below require vars to pass
+    @test HalfSpace(2x ≥ 5y - 1, vars) == HalfSpace([-2.0, 5.0], 1.0)
+    @test HalfSpace(2x >= 5y - 1, vars) == HalfSpace([-2.0, 5.0], 1.0)
+
+    # test with sparse variables
+    @variables x[1:5]
+    @test HalfSpace(2x[1] + 5x[4] <= 10., x) == HalfSpace([2.0, 0.0, 0.0, 5.0, 0.0], 10.0)
+    @test HalfSpace(2x[1] + 5x[4] >= -10. + x[3], x) == HalfSpace([-2.0, 0.0, 1.0, -5.0, 0.0], 10.0)
+end
diff --git a/test/unit_Hyperplane.jl b/test/unit_Hyperplane.jl
@@ -115,3 +115,17 @@ for N in [Float64]
     empty_intersection, v = is_intersection_empty(b, hp, true)
     @test !empty_intersection && !is_intersection_empty(b, hp) && v ∈ hp && v ∈ b
 end
+
+# tests that require ModelingToolkit
+@static if VERSION >= v"1.3" && isdefined(@__MODULE__, :ModelingToolkit)
+    vars = @variables x y
+    @test Hyperplane(2x + 3y == 5) == Hyperplane([2.0, 3.0], 5.0)
+    @test Hyperplane(2x + 3y == 5, N=Int) == Hyperplane([2, 3], 5)
+    @test Hyperplane(2x + 3y == 5, vars) == Hyperplane([2.0, 3.0,], 5.0)
+    @test Hyperplane(2x == 5y) == Hyperplane([2.0, -5.0,], 0.0)
+    @test Hyperplane(2x == 5y, vars) == Hyperplane([2.0, -5.0,], 0.0)
+
+    # test with sparse variables
+    @variables x[1:5]
+    @test Hyperplane(2x[1] + 5x[4] == 10., x) == Hyperplane([2.0, 0.0, 0.0, 5.0, 0.0], 10.0)
+end