Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 22 additions & 10 deletions src/sos_polynomial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,10 @@ function vectorized_matrix(
)
end

function gram_matrix_type(B::Type, matrix_cone_type, T::Type)
return GramMatrix{T,B,T}
end

# Need these two methods to avoid ambiguity
function build_gram_matrix(
q::Vector,
Expand Down Expand Up @@ -70,21 +74,28 @@ end
# return build_gram_matrix(C, basis, T, SymmetricVectorized())
#end

# Need to annotate the element type in case `bases`
# is empty
function _build_matrix_eltype(::Type{T}, bases, f) where {T}
return Base._return_type(f, Tuple{Vector{T},eltype(bases)})
end

function build_matrix(
Q::Function,
bases::Vector{<:AbstractPolynomialBasis},
f::Function,
)
return map(eachindex(bases)) do i
return f(Q(i), bases[i])
end
::Type{T},
) where {T}
@show T
return T[f(Q(i), bases[i]) for i in eachindex(bases)]
end
function build_matrix(
Q::Function,
bases::Vector{<:Vector{<:AbstractPolynomialBasis}},
f::Function,
)
return [
::Type{T},
) where {T}
return T[
f(Q(i), bases[i][j]) for i in eachindex(bases) for
j in eachindex(bases[i])
]
Expand All @@ -96,13 +107,14 @@ function build_gram_matrix(Q::Function, bases::Vector, matrix_cone_type, T)
Q,
bases,
(Q, b) -> build_gram_matrix(Q, b, matrix_cone_type, T),
gram_matrix_type(eltype(bases), matrix_cone_type, T),
),
)
end

function build_moment_matrix(Q::Function, bases::Vector)
return MultivariateMoments.BlockDiagonalMomentMatrix(
build_matrix(Q, bases, build_moment_matrix),
build_matrix(Q, bases, build_moment_matrix, T),
)
end

Expand Down Expand Up @@ -132,10 +144,10 @@ _first(b::AbstractPolynomialBasis) = b
_first(b::Vector) = first(b)
function add_gram_matrix(
model::MOI.ModelLike,
matrix_cone_type::Type,
::Type{matrix_cone_type},
bases::Vector,
T::Type,
)
::Type{T},
) where {matrix_cone_type,T}
Qs = Vector{Vector{MOI.VariableIndex}}(undef, length(bases))
cQs = Vector{union_constraint_indices_types(matrix_cone_type)}(
undef,
Expand Down
38 changes: 38 additions & 0 deletions test/issue_344.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,38 @@
using JuMP
using SumOfSquares
using DynamicPolynomials
function sos_min(sparsity, d, obj, dom)
model =
Model(() -> MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}()))
@variable(model, t)
@objective(model, Min, t)
con_ref = @constraint(
model,
obj - t in SOSCone(),
sparsity = sparsity,
maxdegree = d,
domain = dom
)
return optimize!(model)
#return moment_matrix(con_ref)
end
@polyvar x[1:3] l[1:3]
my_game_obj = x[1]^2 + x[2]^2 + 1.0 * x[3]^2
my_eq_list = [
-1.0 + l[3] - 3l[2] - 3x[2]^2 - x[1] * x[2],
3x[2] * l[2] - x[2]^2 * l[2] - x[1]^2 * l[2],
l[3] - x[3] * l[3],
-2 + x[3] + 2x[1] + 2x[1] * l[1],
-2 - x[3] + 2x[2] + 2x[2] * l[1],
2l[1] - x[3] * l[1] - x[2]^2 * l[1] - x[1]^2 * l[1],
]
my_ineq_list = [
3x[3] - x[2]^2 - x[1]^2,
1 - x[3],
l[2],
l[3],
2 - x[3] - x[2]^2 - x[1]^2,
l[1],
]
my_Dom = basicsemialgebraicset(algebraicset(my_eq_list), my_ineq_list)
ν = sos_min(Sparsity.Monomial(ChordalCompletion()), 4, my_game_obj, my_Dom)