Merge pull request #3532 from JuliaReach/schillic/symengine

`SymEngine` code to create `HalfSpace`/`Hyperplane`

docs/Project.toml

@@ -14,6 +14,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a"
 
@@ -34,4 +35,5 @@ Polyhedra = "0.7"
 RecipesBase = "1"
 StaticArrays = "1"
 Symbolics = "5"
+SymEngine = "0.7 - 0.11"
 TaylorModels = "0.6 - 0.7"

docs/make.jl

@@ -1,6 +1,6 @@
 using Documenter, LazySets
 import Polyhedra, Optim, ExponentialUtilities, TaylorModels, Distributions,
-       MiniQhull, Symbolics
+       MiniQhull, Symbolics, SymEngine
 
 include("init.jl")
 

docs/src/lib/utils.md

@@ -57,6 +57,17 @@ LazySets.PolynomialZonotopeSampler
 LazySets._vec
 ```
 
+## SymEngine
+
+```@docs
+LazySets.free_symbols
+LazySets._is_linearcombination
+LazySets._is_halfspace
+LazySets._is_hyperplane
+LazySets.convert(::Type{HalfSpace{N}}, ::Expr; vars::Vector{Basic}=Basic[]) where {N}
+LazySets.convert(::Type{Hyperplane{N}}, ::Expr; vars::Vector{Basic}=Basic[]) where {N}
+```
+
 ## Functions for numbers
 
 ```@docs

src/Initialization/init_SymEngine.jl

@@ -0,0 +1,84 @@
+using .SymEngine: Basic, subs
+
+import .SymEngine: free_symbols
+
+"""
+   _is_linearcombination(L::Basic)
+
+Determine whether the expression `L` is a linear combination of its symbols.
+
+### Input
+
+- `L` -- expression
+
+### Output
+
+`true` if `L` is a linear combination or false otherwise.
+
+### Examples
+
+```jldoctest
+julia> using LazySets: _is_linearcombination
+
+julia> _is_linearcombination(:(2*x1 - 4))
+true
+
+julia> _is_linearcombination(:(6.1 - 5.3*f - 0.1*g))
+true
+
+julia> _is_linearcombination(:(2*x1^2))
+false
+
+julia> _is_linearcombination(:(x1^2 - 4*x2 + x3 + 2))
+false
+```
+"""
+function _is_linearcombination(L::Basic)
+    return all(isempty.(free_symbols.(diff.(L, free_symbols(L)))))
+end
+
+_is_linearcombination(L::Expr) = _is_linearcombination(convert(Basic, L))
+
+"""
+    free_symbols(expr::Expr, set_type::Type{LazySet})
+
+Return the free symbols in an expression that represents a given set type.
+
+### Input
+
+- `expr` -- symbolic expression
+
+### Output
+
+A list of symbols, in the form of SymEngine `Basic` objects.
+
+### Examples
+
+```jldoctest
+julia> using LazySets: free_symbols
+
+julia> free_symbols(:(x1 <= -0.03), HalfSpace)
+1-element Vector{SymEngine.Basic}:
+ x1
+
+julia> free_symbols(:(x1 + x2 <= 2*x4 + 6), HalfSpace)
+3-element Vector{SymEngine.Basic}:
+ x2
+ x1
+ x4
+```
+"""
+function free_symbols(::Expr, ::Type{<:LazySet}) end
+
+function free_symbols(expr::Expr)
+    if _is_hyperplane(expr)
+        return free_symbols(expr, HyperPlane)
+    elseif _is_halfspace(expr)
+        return free_symbols(expr, Halfspace)
+    else
+        error("the free symbols for the expression $(expr) is not implemented")
+    end
+end
+
+eval(load_symengine_hyperplane())
+eval(load_symengine_halfspace())

src/Sets/HalfSpace.jl

@@ -641,6 +641,154 @@ function load_symbolics_halfspace()
     end
 end  # quote / load_symbolics_halfspace()
 
+# =====================================
+# Functionality that requires SymEngine
+# =====================================
+
+function load_symengine_halfspace()
+    return quote
+        """
+            _is_halfspace(expr::Expr)
+
+        Determine whether the given expression corresponds to a half-space.
+
+        ### Input
+
+        - `expr` -- a symbolic expression
+
+        ### Output
+
+        `true` if `expr` corresponds to a half-space or `false` otherwise.
+
+        ### Examples
+
+        ```jldoctest
+        julia> using LazySets: _is_halfspace
+
+        julia> all(_is_halfspace.([:(x1 <= 0), :(x1 < 0), :(x1 > 0), :(x1 >= 0)]))
+        true
+
+        julia> _is_halfspace(:(2*x1 <= 4))
+        true
+
+        julia> _is_halfspace(:(6.1 <= 5.3*f - 0.1*g))
+        true
+
+        julia> _is_halfspace(:(2*x1^2 <= 4))
+        false
+
+        julia> _is_halfspace(:(x1^2 > 4*x2 - x3))
+        false
+
+        julia> _is_halfspace(:(x1 > 4*x2 - x3))
+        true
+        ```
+        """
+        function _is_halfspace(expr::Expr)::Bool
+
+            # check that there are three arguments
+            # these are the comparison symbol, the left hand side and the right hand side
+            if (length(expr.args) != 3) || !(expr.head == :call)
+                return false
+            end
+
+            # check that this is an inequality
+            if !(expr.args[1] in [:(<=), :(<), :(>=), :(>)])
+                return false
+            end
+
+            # convert to symengine expressions
+            lhs, rhs = convert(Basic, expr.args[2]), convert(Basic, expr.args[3])
+
+            # check if the expression defines a half-space
+            return _is_linearcombination(lhs) && _is_linearcombination(rhs)
+        end
+
+        """
+            convert(::Type{HalfSpace{N}}, expr::Expr; vars=nothing) where {N}
+
+        Return a `LazySet.HalfSpace` given a symbolic expression that represents a half-space.
+
+        ### Input
+
+        - `expr` -- a symbolic expression
+        - `vars` -- (optional, default: `nothing`): set of variables with respect to which
+                    the gradient is taken; if nothing, it takes the free symbols in the given expression
+
+        ### Output
+
+        A `HalfSpace`, in the form `ax <= b`.
+
+        ### Examples
+
+        ```jldoctest convert_halfspace
+        julia> convert(HalfSpace, :(x1 <= -0.03))
+        HalfSpace{Float64, Vector{Float64}}([1.0], -0.03)
+
+        julia> convert(HalfSpace, :(x1 < -0.03))
+        HalfSpace{Float64, Vector{Float64}}([1.0], -0.03)
+
+        julia> convert(HalfSpace, :(x1 > -0.03))
+        HalfSpace{Float64, Vector{Float64}}([-1.0], 0.03)
+
+        julia> convert(HalfSpace, :(x1 >= -0.03))
+        HalfSpace{Float64, Vector{Float64}}([-1.0], 0.03)
+
+        julia> convert(HalfSpace, :(x1 + x2 <= 2*x4 + 6))
+        HalfSpace{Float64, Vector{Float64}}([1.0, 1.0, -2.0], 6.0)
+        ```
+
+        You can also specify the set of "ambient" variables, even if not
+        all of them appear:
+
+        ```jldoctest convert_halfspace
+        julia> using SymEngine: Basic
+
+        julia> convert(HalfSpace, :(x1 + x2 <= 2*x4 + 6), vars=Basic[:x1, :x2, :x3, :x4])
+        HalfSpace{Float64, Vector{Float64}}([1.0, 1.0, 0.0, -2.0], 6.0)
+        ```
+        """
+        function convert(::Type{HalfSpace{N}}, expr::Expr; vars::Vector{Basic}=Basic[]) where {N}
+            @assert _is_halfspace(expr) "the expression :(expr) does not correspond to a half-space"
+
+            # check sense of the inequality, assuming < or <= by default
+            got_geq = expr.args[1] in [:(>=), :(>)]
+
+            # get sides of the inequality
+            lhs, rhs = convert(Basic, expr.args[2]), convert(Basic, expr.args[3])
+
+            # a1 x1 + ... + an xn + K [cmp] 0 for cmp in <, <=, >, >=
+            eq = lhs - rhs
+            if isempty(vars)
+                vars = free_symbols(eq)
+            end
+            K = subs(eq, [vi => zero(N) for vi in vars]...)
+            a = convert(Basic, eq - K)
+
+            # convert to numeric types
+            K = convert(N, K)
+            a = convert(Vector{N}, diff.(a, vars))
+
+            if got_geq
+                return HalfSpace(-a, K)
+            else
+                return HalfSpace(a, -K)
+            end
+        end
+
+        # type-less default half-space conversion
+        function convert(::Type{HalfSpace}, expr::Expr; vars::Vector{Basic}=Basic[])
+            return convert(HalfSpace{Float64}, expr; vars=vars)
+        end
+
+        function free_symbols(expr::Expr, ::Type{<:HalfSpace})
+            # get sides of the inequality
+            lhs, rhs = convert(Basic, expr.args[2]), convert(Basic, expr.args[3])
+            return free_symbols(lhs - rhs)
+        end
+    end
+end  # quote / load_symengine_halfspace()
+
 """
     complement(H::HalfSpace)