Skip to content

Commit

Permalink
Merge pull request #14 from blegat/bl/invariant_polytope_inplace
Browse files Browse the repository at this point in the history
Modify optimizer in-place in invariant_polytope
  • Loading branch information
blegat committed Apr 8, 2020
2 parents 490f44b + 0b78bd0 commit b767d51
Show file tree
Hide file tree
Showing 10 changed files with 246 additions and 92 deletions.
5 changes: 4 additions & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
MathematicalSystems = "d14a8603-c872-5ed3-9ece-53e0e82e39da"
MultivariateMoments = "f4abf1af-0426-5881-a0da-e2f168889b5e"
MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
Optim = "429524aa-4258-5aef-a3af-852621145aeb"
ParameterJuMP = "774612a8-9878-5177-865a-ca53ae2495f9"
Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
Expand Down Expand Up @@ -53,8 +54,10 @@ julia = "1"

[extras]
CSDP = "0a46da34-8e4b-519e-b418-48813639ff34"
ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["LinearAlgebra", "Test", "CSDP"]
test = ["LinearAlgebra", "Test", "CSDP", "ECOS", "GLPK"]
3 changes: 3 additions & 0 deletions src/SwitchOnSafety.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,9 @@ module SwitchOnSafety

using LinearAlgebra, SparseArrays, Statistics

import MutableArithmetics
const MA = MutableArithmetics

using DynamicPolynomials
using MultivariatePolynomials
using SemialgebraicSets
Expand Down
105 changes: 65 additions & 40 deletions src/balanced_complex_polytope.jl
Original file line number Diff line number Diff line change
@@ -1,69 +1,94 @@
using JuMP
using ParameterJuMP
using Polyhedra
mutable struct BalancedComplexPolytope{T, VT <: AbstractVector{T}, D<:Polyhedra.FullDim}
mutable struct BalancedComplexPolytope{T, CT, VT <: AbstractVector{CT}, D<:Polyhedra.FullDim}
d::D
points::Vector{VT}
optimizer_constructor
model::Union{Nothing, JuMP.Model}
real_z::Union{Nothing, Vector{ParameterJuMP.ParameterRef}}
imag_z::Union{Nothing, Vector{ParameterJuMP.ParameterRef}}
t_0::Union{Nothing, JuMP.VariableRef}
function BalancedComplexPolytope{T, VT, D}(
model::Union{Nothing, MOI.ModelLike}
sum_con::Union{Nothing, MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
real_z::Union{Nothing, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}}
imag_z::Union{Nothing, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}}
t_0::Union{Nothing, MOI.VariableIndex}
function BalancedComplexPolytope{T, CT, VT, D}(
d::Polyhedra.FullDim, points::Polyhedra.PointIt,
optimizer_constructor=nothing) where {T, VT, D}
new{T, VT, D}(Polyhedra.FullDim_convert(D, d),
Polyhedra.lazy_collect(points), optimizer_constructor, nothing, nothing, nothing, nothing)
optimizer_constructor=nothing) where {T, CT, VT, D}
new{T, CT, VT, D}(Polyhedra.FullDim_convert(D, d),
Polyhedra.lazy_collect(points), optimizer_constructor, nothing, nothing, nothing, nothing)
end
end
function BalancedComplexPolytope(d::Polyhedra.FullDim, points::Polyhedra.PointIt, args...)
return BalancedComplexPolytope{Polyhedra.coefficient_type(eltype(points)), eltype(points),
typeof(d)}(d, points, args...)
VT = eltype(points)
CT = Polyhedra.coefficient_type(VT)
return BalancedComplexPolytope{real(CT), CT, VT, typeof(d)}(d, points, args...)
end

Polyhedra.FullDim(v::BalancedComplexPolytope) = v.d

function _reset_model(p::BalancedComplexPolytope)
p.model = nothing
p.real_z = nothing
p.imag_z = nothing
p.t_0 = nothing
end

function _build_model(brp::BalancedComplexPolytope)
function _build_model(brp::BalancedComplexPolytope{T}) where T
@assert brp.model === nothing
n = length(brp.points)
brp.model = ParameterJuMP.ModelWithParams(brp.optimizer_constructor)
brp.model = MOI.instantiate(brp.optimizer_constructor, with_bridge_type = T)
# Real part of `t`
α = @variable(brp.model, [1:n])
α = Vector{MOI.VariableIndex}(undef, n)
# Imaginary part of `t`
β = @variable(brp.model, [1:n])
q = @variable(brp.model, [1:n])
β = Vector{MOI.VariableIndex}(undef, n)
q = Vector{MOI.VariableIndex}(undef, n)
# |t| ≤ q
for i in 1:n
@constraint(brp.model, [q[i], α[i], β[i]] in SecondOrderCone())
vars, con = MOI.add_constrained_variables(brp.model, MOI.SecondOrderCone(3))
q[i], α[i], β[i] = vars
end
# sum |t| ≤ sum q ≤ 1
@constraint(brp.model, sum(q) 1)
return sum(α[i] .* real.(brp.points[i]) for i in 1:n) - sum(β[i] .* imag.(brp.points[i]) for i in 1:n),
sum(α[i] .* imag.(brp.points[i]) for i in 1:n) + sum(β[i] .* real.(brp.points[i]) for i in 1:n)
end

function _fix(p::BalancedComplexPolytope, v)
JuMP.fix.(brp.real_z, real.(v))
JuMP.fix.(brp.imag_z, imag.(v))
brp.sum_con = MOI.add_constraint(brp.model, _sum(q, T), MOI.LessThan(one(T)))
rp = [real.(p) for p in brp.points]
ip = [imag.(p) for p in brp.points]
real_hull = MOI.Utilities.operate!.(
-, T,
_convex_combination(α, rp, brp.d),
_convex_combination(β, ip, brp.d))
imag_hull = MOI.Utilities.operate!.(
+, T,
_convex_combination(α, ip, brp.d),
_convex_combination(β, rp, brp.d))
return real_hull, imag_hull
end

function _parametrized_model(brp::BalancedComplexPolytope, v)
real_hull, imag_hull = _build_model(brp)
brp.real_z = ParameterJuMP.add_parameters(brp.model, collect(real.(v)))
brp.imag_z = ParameterJuMP.add_parameters(brp.model, collect(imag.(v)))
@constraint(brp.model, brp.real_z .== real_hull)
@constraint(brp.model, brp.imag_z .== imag_hull)
brp.real_z = [MOI.add_constraint(brp.model, real_hull[i], MOI.EqualTo(real(v[i]))) for i in 1:brp.d]
brp.imag_z = [MOI.add_constraint(brp.model, imag_hull[i], MOI.EqualTo(imag(v[i]))) for i in 1:brp.d]
end

function _ratio_model(brp::BalancedComplexPolytope, v)
function _ratio_model(brp::BalancedComplexPolytope{T}, v) where T
real_hull, imag_hull = _build_model(brp)
brp.t_0 = @variable(brp.model)
@constraint(brp.model, brp.t_0 .* real.(v) .== real_hull)
@constraint(brp.model, brp.t_0 .* imag.(v) .== imag_hull)
brp.t_0 = MOI.add_variable(brp.model)
for i in 1:brp.d
push!(real_hull[i].terms, MOI.ScalarAffineTerm(-real(v[i]), brp.t_0))
push!(imag_hull[i].terms, MOI.ScalarAffineTerm(-imag(v[i]), brp.t_0))
end
brp.real_z = [MOI.add_constraint(brp.model, real_hull[i], MOI.EqualTo(zero(T))) for i in 1:brp.d]
brp.imag_z = [MOI.add_constraint(brp.model, imag_hull[i], MOI.EqualTo(zero(T))) for i in 1:brp.d]
end

function _update_model(p::BalancedComplexPolytope{T}, v) where T
vars, con = MOI.add_constrained_variables(p.model, MOI.SecondOrderCone(3))
q, α, β = vars
MOI.modify(p.model, p.sum_con, MOI.ScalarCoefficientChange(q, one(T)))
for i in 1:p.d
MOI.modify(p.model, p.real_z[i], MOI.ScalarCoefficientChange(α, real(v[i])))
MOI.modify(p.model, p.imag_z[i], MOI.ScalarCoefficientChange(β, imag(v[i])))
end
end

function _fix(bcp::BalancedComplexPolytope, v)
for i in 1:bcp.d
if bcp.t_0 === nothing
MOI.set(bcp.model, MOI.ConstraintSet(), bcp.real_z[i], MOI.EqualTo(real(v[i])))
MOI.set(bcp.model, MOI.ConstraintSet(), bcp.imag_z[i], MOI.EqualTo(imag(v[i])))
else
MOI.modify(bcp.model, bcp.real_z[i], MOI.ScalarCoefficientChange(bcp.t_0, -real(v[i])))
MOI.modify(bcp.model, bcp.imag_z[i], MOI.ScalarCoefficientChange(bcp.t_0, -imag(v[i])))
end
end
end
66 changes: 48 additions & 18 deletions src/balanced_real_polytope.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,14 @@ mutable struct BalancedRealPolytope{T, VT <: AbstractVector{T}, D<:Polyhedra.Ful
d::D
points::Vector{VT}
optimizer_constructor
model::Union{Nothing, JuMP.Model}
z::Union{Nothing, Vector{ParameterJuMP.ParameterRef}}
t_0::Union{Nothing, JuMP.VariableRef}
model::Union{Nothing, MOI.ModelLike}
sum_con::Union{Nothing, MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
z::Union{Nothing, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}}
t_0::Union{Nothing, MOI.VariableIndex}
function BalancedRealPolytope{T, VT, D}(
d::Polyhedra.FullDim, points::Polyhedra.PointIt,
optimizer_constructor=nothing) where {T, VT, D}
new{T, VT, D}(Polyhedra.FullDim_convert(D, d), Polyhedra.lazy_collect(points), optimizer_constructor, nothing, nothing, nothing)
new{T, VT, D}(Polyhedra.FullDim_convert(D, d), Polyhedra.lazy_collect(points), optimizer_constructor, nothing, nothing, nothing, nothing)
end
end
function BalancedRealPolytope(d::Polyhedra.FullDim, points::Polyhedra.PointIt, args...)
Expand Down Expand Up @@ -56,27 +57,56 @@ function Polyhedra.nextindex(rep::BalancedRealPolytope{T}, idx::Polyhedra.PointI
end
end

function _build_model(brp::BalancedRealPolytope)
function _abs_constraint(model::MOI.ModelLike, t::MOI.VariableIndex, q::MOI.VariableIndex, T::Type)
ft = MOI.SingleVariable(t)
fq = MOI.SingleVariable(q)
MOI.add_constraint(model, MOI.Utilities.operate(-, T, fq, ft), MOI.GreaterThan(zero(T)))
MOI.add_constraint(model, MOI.Utilities.operate(+, T, fq, ft), MOI.GreaterThan(zero(T)))
end

function _build_model(brp::BalancedRealPolytope{T}) where T
@assert brp.model === nothing
n = length(brp.points)
brp.model = ParameterJuMP.ModelWithParams(brp.optimizer_constructor)
t = @variable(brp.model, [1:n])
q = @variable(brp.model, [1:n])
brp.model = MOI.instantiate(brp.optimizer_constructor, with_bridge_type = T)
t = MOI.add_variables(brp.model, n)
q = MOI.add_variables(brp.model, n)
# |t| ≤ q
@constraint(brp.model, t .≤ q)
@constraint(brp.model, -t .≤ q)
for i in 1:n
_abs_constraint(brp.model, t[i], q[i], T)
end
# sum |t| ≤ sum q ≤ 1
@constraint(brp.model, sum(q) 1)
return sum(t[i] .* brp.points[i] for i in 1:n)
brp.sum_con = MOI.add_constraint(brp.model, _sum(q, T), MOI.LessThan(one(T)))
return _convex_combination(t, brp.points, brp.d)
end

function _parametrized_model(brp::BalancedRealPolytope, v)
hull = _build_model(brp)
brp.z = ParameterJuMP.add_parameters(brp.model, collect(v))
@constraint(brp.model, brp.z .== hull)
brp.z = [MOI.add_constraint(brp.model, hull[i], MOI.EqualTo(v[i])) for i in 1:brp.d]
end

function _ratio_model(brp::BalancedRealPolytope, v)
function _ratio_model(brp::BalancedRealPolytope{T}, v) where T
hull = _build_model(brp)
brp.t_0 = @variable(brp.model)
@constraint(brp.model, brp.t_0 .* v .== hull)
brp.t_0 = MOI.add_variable(brp.model)
for i in 1:brp.d
push!(hull[i].terms, MOI.ScalarAffineTerm(-v[i], brp.t_0))
end
brp.z = [MOI.add_constraint(brp.model, hull[i], MOI.EqualTo(zero(T))) for i in 1:brp.d]
end
function _update_model(p::BalancedRealPolytope{T}, v) where T
t = MOI.add_variable(p.model)
q = MOI.add_variable(p.model)
_abs_constraint(p.model, t, q, T)
MOI.modify(p.model, p.sum_con, MOI.ScalarCoefficientChange(q, one(T)))
for i in 1:p.d
MOI.modify(p.model, p.z[i], MOI.ScalarCoefficientChange(t, v[i]))
end
end

function _fix(brp::BalancedRealPolytope, v)
for i in 1:brp.d
if brp.t_0 === nothing
MOI.set(brp.model, MOI.ConstraintSet(), brp.z[i], MOI.EqualTo(v[i]))
else
MOI.modify(brp.model, brp.z[i], MOI.ScalarCoefficientChange(brp.t_0, -v[i]))
end
end
end
58 changes: 45 additions & 13 deletions src/conitope.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,10 @@ mutable struct Conitope{T, VT <: AbstractVector{T}, D<:Polyhedra.FullDim}
d::D
points::Vector{VT}
optimizer_constructor
model::Union{Nothing, JuMP.Model}
z::Union{Nothing, Vector{ParameterJuMP.ParameterRef}}
t_0::Union{Nothing, JuMP.VariableRef}
model::Union{Nothing, MOI.ModelLike}
sum_con::Union{Nothing, MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
z::Union{Nothing, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}}}
t_0::Union{Nothing, MOI.VariableIndex}
function Conitope{T, VT, D}(
d::Polyhedra.FullDim, points::Polyhedra.PointIt,
optimizer_constructor=nothing) where {T, VT, D}
Expand All @@ -23,22 +24,53 @@ end

Polyhedra.FullDim(v::Conitope) = v.d

function _build_model(brp::Conitope)
function _sum(x::Vector{MOI.VariableIndex}, T::Type)
MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(one(T), xi) for xi in x],
zero(T))
end
function _convex_combination(t::Vector{MOI.VariableIndex}, points::Vector{Vector{T}}, d::Int) where T
return [MOI.ScalarAffineFunction([
MOI.ScalarAffineTerm(points[i][j], t[i]) for i in eachindex(t)
], zero(T)) for j in 1:d]
end
function _build_model(brp::Conitope{T}) where T
@assert brp.model === nothing
n = length(brp.points)
brp.model = ParameterJuMP.ModelWithParams(brp.optimizer_constructor)
@variable(brp.model, t[1:n] 0)
@constraint(brp.model, sum(t) 1)
return sum(t[i] .* brp.points[i] for i in 1:n)
brp.model = MOI.instantiate(brp.optimizer_constructor, with_bridge_type = T)
t, ct = MOI.add_constrained_variables(brp.model, MOI.Nonnegatives(n))
# sum t_i ≤ 1
brp.sum_con = MOI.add_constraint(brp.model, _sum(t, T), MOI.LessThan(one(T)))
return _convex_combination(t, brp.points, brp.d)
end
function _parametrized_model(brp::Conitope, v)
hull = _build_model(brp)
brp.z = ParameterJuMP.add_parameters(brp.model, collect(v))
@constraint(brp.model, brp.z .≤ hull)
brp.z = [MOI.add_constraint(brp.model, hull[i], MOI.GreaterThan(v[i])) for i in 1:brp.d]
end
function _ratio_model(brp::Conitope, v)
function _ratio_model(brp::Conitope{T}, v) where T
hull = _build_model(brp)
brp.t_0 = @variable(brp.model)
@constraint(brp.model, brp.t_0 .* v .≤ hull)
brp.t_0 = MOI.add_variable(brp.model)
for i in 1:brp.d
push!(hull[i].terms, MOI.ScalarAffineTerm(-v[i], brp.t_0))
end
brp.z = [MOI.add_constraint(brp.model, hull[i], MOI.GreaterThan(zero(T))) for i in 1:brp.d]
end
function _update_model(p::Conitope{T}, v) where T
ts, ct = MOI.add_constrained_variables(p.model, MOI.Nonnegatives(1))
t = ts[1]
MOI.modify(p.model, p.sum_con, MOI.ScalarCoefficientChange(t, one(T)))
for i in 1:p.d
MOI.modify(p.model, p.z[i], MOI.ScalarCoefficientChange(t, v[i]))
end
end

function _fix(brp::Conitope, v)
for i in 1:brp.d
if brp.t_0 === nothing
MOI.set(brp.model, MOI.ConstraintSet(), brp.z[i], MOI.GreaterThan(v[i]))
else
MOI.modify(brp.model, brp.z[i], MOI.ScalarCoefficientChange(brp.t_0, -v[i]))
end
end
end

function _ei(i::Integer, d::Integer, T::Type)
Expand Down
37 changes: 18 additions & 19 deletions src/polytopes.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,18 +11,11 @@ include("balanced_complex_polytope.jl")

const PolytopeLike = Union{Conitope, BalancedRealPolytope, BalancedComplexPolytope}

function _reset_model(p::PolytopeLike)
p.model = nothing
p.z = nothing
p.t_0 = nothing
end
function Polyhedra.convexhull!(p::PolytopeLike, v::SymPoint)
push!(p.points, v.point)
# Invalidates the model
_reset_model(p)
end
function _fix(p::PolytopeLike, v)
JuMP.fix.(p.z, v)
if p.model !== nothing
_update_model(p, v.point)
end
end
function Base.in(v, brp::PolytopeLike)
if isempty(brp.points)
Expand All @@ -33,11 +26,13 @@ function Base.in(v, brp::PolytopeLike)
else
_fix(brp, v)
end
JuMP.optimize!(brp.model)
status = termination_status(brp.model)
MOI.optimize!(brp.model)
status = MOI.get(brp.model, MOI.TerminationStatus())
if status == MOI.OPTIMAL
return true
elseif status == MOI.INFEASIBLE || status == MOI.INFEASIBLE_OR_UNBOUNDED
# There is no objective so it cannot be unbounded.
# Therefore, `MOI.INFEASIBLE_OR_UNBOUNDED` means infeasible.
return false
else
error("Solver returned status $status.")
Expand All @@ -49,16 +44,20 @@ function in_ratio(v, p::PolytopeLike)
# multiplied by 0.0 to be mapped to the polytope
return 0.0
end
_ratio_model(p, v)
@objective(p.model, Max, p.t_0)
JuMP.optimize!(p.model)
status = termination_status(p.model)
if p.model === nothing
_ratio_model(p, v)
MOI.set(p.model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
MOI.set(p.model, MOI.ObjectiveFunction{MOI.SingleVariable}(),
MOI.SingleVariable(p.t_0))
else
_fix(p, v)
end
MOI.optimize!(p.model)
status = MOI.get(p.model, MOI.TerminationStatus())
if status == MOI.OPTIMAL
r = JuMP.value(p.t_0)
_reset_model(p)
r = MOI.get(p.model, MOI.VariablePrimal(), p.t_0)
return r
else
_reset_model(p)
error("Solver returned status $status.")
end
end
Expand Down

2 comments on commit b767d51

@blegat
Copy link
Owner Author

@blegat blegat commented on b767d51 Apr 8, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Error while trying to register: "Tag with name 0.1.0 already exists and points to a different commit"

Please sign in to comment.