diff --git a/src/sos_polynomial.jl b/src/sos_polynomial.jl index 7e95de754..35eec6144 100644 --- a/src/sos_polynomial.jl +++ b/src/sos_polynomial.jl @@ -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, @@ -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]) ] @@ -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 @@ -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, diff --git a/test/issue_344.jl b/test/issue_344.jl new file mode 100644 index 000000000..8d79721ed --- /dev/null +++ b/test/issue_344.jl @@ -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)