Merge pull request #83 from JuliaReach/lf-polynomials

update sdp and polynomial interface

Project.toml

@@ -3,14 +3,12 @@ uuid = "1b4d18b6-9e5d-11e9-236c-f792b01831f8"
 version = "0.1.1"
 
 [deps]
-DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 
 [compat]
-DynamicPolynomials = "0.3, 0.4"
 ForwardDiff = "0.10"
 IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
 IntervalOptimisation = "0.4.1"
@@ -22,6 +20,7 @@ julia = "1.6"
 
 [extras]
 IntervalOptimisation = "c7c68f13-a4a2-5b9a-b424-07d005f8d9d2"
+MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/src/lib/methods.md

@@ -21,10 +21,4 @@ enclose
 
 ```@docs
 relative_precision
-```
-
-## Specifics to sum-of-squares optimization
-
-```@docs
-new_sos
-```
+```

src/RangeEnclosures.jl

@@ -1,6 +1,6 @@
 module RangeEnclosures
 
-using DynamicPolynomials, Requires, Reexport
+using Requires, Reexport
 @reexport using ForwardDiff
 @reexport using IntervalArithmetic
 const Interval_or_IntervalBox = Union{Interval, IntervalBox}
@@ -17,6 +17,7 @@ function __init__()
     @require SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" include("sdp.jl")
     @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" include("taylormodels.jl")
     @require IntervalOptimisation = "c7c68f13-a4a2-5b9a-b424-07d005f8d9d2" include("intervaloptimisation.jl")
+    @require MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" include("polynomials.jl")
 end
 
 #================

src/enclose.jl

@@ -43,12 +43,6 @@ function enclose(f::Function, dom::Interval_or_IntervalBox,
     return _enclose(solver, f, dom; kwargs...)
 end
 
-function enclose(p::AbstractPolynomialLike, dom::Interval_or_IntervalBox,
-                 solver::AbstractEnclosureAlgorithm=NaturalEnclosure(); kwargs...)
-    f(x...) = p(variables(p) => x)
-    return _enclose(solver, f, dom; kwargs...)
-end
-
 function enclose(f::Function, dom::Interval_or_IntervalBox,
                  method::Vector; kwargs...)
    return mapreduce(ξ -> _enclose(ξ, f, dom; kwargs...), ∩, method)

src/polynomials.jl

@@ -0,0 +1,12 @@
+using .MultivariatePolynomials
+
+function enclose(p::AbstractPolynomialLike, dom::Interval_or_IntervalBox,
+                 solver::AbstractEnclosureAlgorithm=NaturalEnclosure(); kwargs...)
+    f(x...) = p(variables(p) => x)
+    return _enclose(solver, f, dom; kwargs...)
+end
+
+function enclose(p::AbstractPolynomialLike, dom::Interval_or_IntervalBox,
+                 solver::SumOfSquaresEnclosure; kwargs...)
+    return _enclose(solver, p, dom; kwargs...)
+end

src/sdp.jl

@@ -1,100 +1,31 @@
-#======================================
-Methods using semidefinite programming
-======================================#
+# #======================================
+# Methods using semidefinite programming
+# ======================================#
 using .SumOfSquares
 
-"""
-    new_sos(backend; kwargs...)
-
-Return a new (empty) sum-of-squares optimization problem for the given backend.
-
-### Input
-
-- `backend` -- the backend (also called JuMP `solver`)
-- `kwargs`  -- additional keyword arguments
-
-### Output
-
-An instance of `SOSModel` for the given backend and options.
-"""
-@inline function new_sos(backend, kwargs...)
-    𝑂 = Dict(kwargs)
-
-    if :QUIET ∈ keys(𝑂)
-        # for mosek solver
-        SOSModel(backend, QUIET=𝑂[:QUIET])
-    else
-        SOSModel(backend)
-    end
-end
-
-
-function _enclose(sose::SumOfSquaresEnclosure, f::Function, dom::Interval_or_IntervalBox;
+function _enclose(sose::SumOfSquaresEnclosure, p::AbstractPolynomialLike, dom::Interval_or_IntervalBox;
                   kwargs...)
 
-    _enclose_SumOfSquares(f, dom, sose.order, sose.backend; kwargs...)
-end
-
-function _enclose_SumOfSquares(f::Function, dom::Interval, order::Int,
-                               backend; kwargs...)
-
-    # polynomial variables
-    @polyvar x
-    p = f(x)
-
-    # box constraints
-    B = @set inf(dom) <= x && x <= sup(dom)
-
-    # ============
-    # Upper bound
-    # ============
-    model = new_sos(backend, kwargs...)
-    @variable(model, γ) # JuMP decision variable
-    @constraint(model, p <= γ, domain=B, maxdegree=order)
-    @objective(model, Min, γ)
-    optimize!(model)
-    upper_bound = objective_value(model)
-
-    # ============
-    # Lower bound
-    # ============
-    model = new_sos(backend, kwargs...)
-    @variable(model, γ) # JuMP decision variable
-    @constraint(model, p >= γ, domain=B, maxdegree=order)
-    @objective(model, Max, γ)
-    optimize!(model)
-    lower_bound = objective_value(model)
-
-    return Interval(lower_bound, upper_bound)
-end
-
-function _enclose_SumOfSquares(f::Function, dom::IntervalBox{N}, order::Int,
-                               backend; kwargs...) where {N}
-
-    # polynomial variables
-    @polyvar x[1:N]
-    p = f(x...)
+    x = variables(p)
 
-    # box constraints
-    Bi =[@set inf(dom[i]) <= x[i] && x[i] <= sup(dom[i]) for i in 1:N]
-    B = reduce(intersect, Bi)
+    B = reduce(intersect, @set inf(domi) <= xi && xi <= sup(domi) for (xi, domi) in zip(x, dom))
 
     # ============
     # Upper bound
     # ============
-    model = new_sos(backend, kwargs...)
+    model = SOSModel(sose.backend, kwargs...)
     @variable(model, γ) # JuMP decision variable
-    @constraint(model, p <= γ, domain=B, maxdegree=order)
+    @constraint(model, p <= γ, domain=B, maxdegree=sose.order)
     @objective(model, Min, γ)
     optimize!(model)
     upper_bound = objective_value(model)
 
     # ============
     # Lower bound
     # ============
-    model = new_sos(backend, kwargs...)
+    model = SOSModel(sose.backend, kwargs...)
     @variable(model, γ) # JuMP decision variable
-    @constraint(model, p >= γ, domain=B, maxdegree=order)
+    @constraint(model, p >= γ, domain=B, maxdegree=sose.order)
     @objective(model, Max, γ)
     optimize!(model)
     lower_bound = objective_value(model)

test/multivariate.jl

@@ -31,4 +31,8 @@ end
     @test rleft ≤ 1e-5 && rright ≤ 1e-5
     # Note: DynamicPolynomials automatically expands p, and evaluation using
     # interval arithmetic gives a worse left bound than the factored expression.
+
+    x = enclose(p, dom, SumOfSquaresEnclosure(backend=SDPA.Optimizer))
+    @test isapprox(inf(x), 0.0; atol=1e-3)
+    @test isapprox(sup(x), 670.612; atol=1e-3)
 end

test/univariate.jl

@@ -37,11 +37,6 @@ end
     xref = Interval(4.83299, 10.5448)
     rleft, rright = relative_precision(x, xref)
     @test rleft ≤ 1e-5 && rright ≤ 1e-5
-
-    x = enclose(f, dom, SumOfSquaresEnclosure(backend=SDPA.Optimizer))
-    xref = Interval(4.8333, 10.541)
-    rleft, rright = relative_precision(x, xref)
-    @test rleft ≤ 1e-5 && rright ≤ 1e-5
 end
 
 @testset "Test univariate polynomial input" begin
@@ -53,4 +48,9 @@ end
     xref = Interval(-5.66667, 19.8334)
     rleft, rright = relative_precision(x, xref)
     @test rleft ≤ 1e-5 && rright ≤ 1e-5
+
+    x = enclose(p, dom, SumOfSquaresEnclosure(backend=SDPA.Optimizer))
+    xref = Interval(4.8333, 10.541)
+    rleft, rright = relative_precision(x, xref)
+    @test rleft ≤ 1e-5 && rright ≤ 1e-5
 end