diff --git a/docs/Project.toml b/docs/Project.toml index c48a410c18..204adae797 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/make.jl b/docs/make.jl index 383ed247b5..0834dbc43c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,6 +1,6 @@ using Documenter, LazySets import Polyhedra, Optim, ExponentialUtilities, TaylorModels, Distributions, - MiniQhull, Symbolics + MiniQhull, Symbolics, SymEngine include("init.jl") diff --git a/docs/src/lib/utils.md b/docs/src/lib/utils.md index 587f7300eb..aa0a921e6c 100644 --- a/docs/src/lib/utils.md +++ b/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 diff --git a/src/Initialization/init_SymEngine.jl b/src/Initialization/init_SymEngine.jl new file mode 100644 index 0000000000..5152de0572 --- /dev/null +++ b/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()) diff --git a/src/Sets/HalfSpace.jl b/src/Sets/HalfSpace.jl index 489f40a249..b0d6b6c614 100644 --- a/src/Sets/HalfSpace.jl +++ b/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) diff --git a/src/Sets/Hyperplane.jl b/src/Sets/Hyperplane.jl index 27b313adb4..aaf4c87096 100644 --- a/src/Sets/Hyperplane.jl +++ b/src/Sets/Hyperplane.jl @@ -612,6 +612,151 @@ function load_symbolics_hyperplane() end end # quote / load_symbolics_hyperplane() +# ===================================== +# Functionality that requires SymEngine +# ===================================== + +function load_symengine_hyperplane() + return quote + """ + _is_hyperplane(expr::Expr) + + Determine whether the given expression corresponds to a hyperplane. + + ### Input + + - `expr` -- a symbolic expression + + ### Output + + `true` if `expr` corresponds to a half-space or `false` otherwise. + + ### Examples + + ```jldoctest + julia> using LazySets: _is_hyperplane + + julia> _is_hyperplane(:(x1 = 0)) + true + + julia> _is_hyperplane(:(2*x1 = 4)) + true + + julia> _is_hyperplane(:(6.1 = 5.3*f - 0.1*g)) + true + + julia> _is_hyperplane(:(2*x1^2 = 4)) + false + + julia> _is_hyperplane(:(x1^2 = 4*x2 - x3)) + false + + julia> _is_hyperplane(:(x1 = 4*x2 - x3)) + true + ``` + """ + function _is_hyperplane(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) != 2) || !(expr.head == :(=)) + return false + end + + # convert to symengine expressions + lhs = convert(Basic, expr.args[1]) + + if :args in fieldnames(typeof(expr.args[2])) + # treats the 4 in :(2*x1 = 4) + rhs = convert(Basic, expr.args[2].args[2]) + else + rhs = convert(Basic, expr.args[2]) + end + + # check if the expression defines a hyperplane + return _is_linearcombination(lhs) && _is_linearcombination(rhs) + end + + """ + convert(::Type{Hyperplane{N}}, expr::Expr; vars=nothing) where {N} + + Return a `LazySet.Hyperplane` given a symbolic expression that represents a hyperplane. + + ### 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 `Hyperplane`, in the form `ax = b`. + + ### Examples + + ```jldoctest convert_hyperplane + julia> convert(Hyperplane, :(x1 = -0.03)) + Hyperplane{Float64, Vector{Float64}}([1.0], -0.03) + + julia> convert(Hyperplane, :(x1 + 0.03 = 0)) + Hyperplane{Float64, Vector{Float64}}([1.0], -0.03) + + julia> convert(Hyperplane, :(x1 + x2 = 2*x4 + 6)) + Hyperplane{Float64, Vector{Float64}}([1.0, 1.0, -2.0], 6.0) + ``` + + You can also specify the set of "ambient" variables in the hyperplane, even if not + all of them appear: + + ```jldoctest convert_hyperplane + julia> using SymEngine: Basic + + julia> convert(Hyperplane, :(x1 + x2 = 2*x4 + 6), vars=Basic[:x1, :x2, :x3, :x4]) + Hyperplane{Float64, Vector{Float64}}([1.0, 1.0, 0.0, -2.0], 6.0) + ``` + """ + function convert(::Type{Hyperplane{N}}, expr::Expr; vars::Vector{Basic}=Basic[]) where {N} + @assert _is_hyperplane(expr) "the expression :(expr) does not correspond to a Hyperplane" + + # get sides of the inequality + lhs = convert(Basic, expr.args[1]) + + # treats the 4 in :(2*x1 = 4) + rhs = :args in fieldnames(typeof(expr.args[2])) ? convert(Basic, expr.args[2].args[2]) : + convert(Basic, expr.args[2]) + + # a1 x1 + ... + an xn + K = 0 + 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)) + + return Hyperplane(a, -K) + end + + # type-less default Hyperplane conversion + function convert(::Type{Hyperplane}, expr::Expr; vars::Vector{Basic}=Basic[]) + return convert(Hyperplane{Float64}, expr; vars=vars) + end + + function free_symbols(expr::Expr, ::Type{<:Hyperplane}) + # get sides of the equality + lhs = convert(Basic, expr.args[1]) + + # treats the 4 in :(2*x1 = 4) + rhs = :args in fieldnames(typeof(expr.args[2])) ? convert(Basic, expr.args[2].args[2]) : + convert(Basic, expr.args[2]) + return free_symbols(lhs - rhs) + end + end +end # quote / load_symengine_hyperplane() + """ distance(x::AbstractVector, H::Hyperplane{N}) where {N} diff --git a/src/init.jl b/src/init.jl index f42f17a174..aa7c16a32e 100644 --- a/src/init.jl +++ b/src/init.jl @@ -12,6 +12,7 @@ function __init__() @require SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" include("Initialization/init_SCS.jl") # @require StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" include("Initialization/init_StaticArraysCore.jl") @require RangeEnclosures = "1b4d18b6-9e5d-11e9-236c-f792b01831f8" include("Initialization/init_RangeEnclosures.jl") + @require SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" include("Initialization/init_SymEngine.jl") @require Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" include("Initialization/init_Symbolics.jl") @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" include("Initialization/init_TaylorModels.jl") @require WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" include("Initialization/init_WriteVTK.jl")