Update to StarAlgebras

Project.toml

@@ -7,10 +7,15 @@ version = "0.0.1-DEV"
 AutoGrad = "6710c13c-97f1-543f-91c5-74e8f7d95b35"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
+MultivariateMoments = "f4abf1af-0426-5881-a0da-e2f168889b5e"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
+PolyJuMP = "ddf597a6-d67e-5340-b84c-e37d84115374"
 SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204"
+StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
+SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
 TrigPolys = "bbdedc48-cb31-4a37-9fe3-b015aecc8dd3"
 
 [compat]

src/BMSOS.jl

@@ -5,6 +5,8 @@ using LinearAlgebra
 import NLopt
 import TrigPolys
 import AutoGrad
+import StarAlgebras as SA
+import MultivariateBases as MB
 
 """
     fgradfft(U::AbstractArray, p::TrigPolys.TrigPoly; mapfn = Base.map)
@@ -157,6 +159,17 @@ function sos_decomp(
     return (; fvals, Uopt, ret, stats)
 end
 
+function sos_decomp(p::SA.AlgebraElement; args...)
+    d = div(MP.maxdegree(SA.basis(p)), 2)
+    @assert MP.nvariables(MB.explicit_basis(p)) == 1
+    a = SA.coeffs(
+        p,
+        MB.SubBasis{MB.Trigonometric}(MP.monomials(MP.variables(p), 0:2d)),
+    )
+    t = TrigPolys.TrigPoly(a[1], a[2:2:2d], a[3:2:(2d+1)])
+    return sos_decomp(t; args...)
+end
+
 function sos_opt(
     p::TrigPolys.TrigPoly;
     rank = 2,

src/MOI_wrapper.jl

@@ -14,16 +14,28 @@ const DEFAULT = Dict{String,Any}(
 )
 
 mutable struct Optimizer <: MOI.AbstractOptimizer
-    poly::Union{Nothing,MP.AbstractPolynomial}
+    poly::Union{Nothing,SA.AlgebraElement}
     options::Dict{String,Any}
     fvals::Union{Nothing,Vector{Float64}}
     Uopt::Union{Nothing,Matrix{Float64}}
     ret::Union{Nothing,Symbol}
+    solve_time::Float64
+    silent::Bool
     function Optimizer()
-        return new(nothing, copy(DEFAULT), nothing, nothing, nothing, nothing)
+        return new(
+            nothing,
+            copy(DEFAULT),
+            nothing,
+            nothing,
+            nothing,
+            NaN,
+            false,
+        )
     end
 end
 
+MOI.get(::Optimizer, ::MOI.SolverName) = "BMSOS"
+
 function MOI.supports(optimizer::Optimizer, attr::MOI.RawOptimizerAttribute)
     return haskey(optimizer.options, attr.name)
 end
@@ -41,6 +53,13 @@ function MOI.set(optimizer::Optimizer, attr::MOI.RawOptimizerAttribute, value)
     return
 end
 
+MOI.supports(::Optimizer, ::MOI.Silent) = true
+
+function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
+    optimizer.silent = value
+    return
+end
+
 function MOI.supports_constraint(
     ::Optimizer,
     ::Type{MOI.VectorAffineFunction{Float64}},
@@ -57,24 +76,38 @@ function MOI.add_constraint(
         error("Nonconstant polynomials are not supported yet!")
     end
     # FIXME don't ignore `set.certificate`
-    optimizer.poly = MP.polynomial(func.constants, set.monomials)
+    optimizer.poly = MB.algebra_element(func.constants, set.basis)
     return MOI.ConstraintIndex{typeof(func),typeof(set)}(0)
 end
 
 function MOI.empty!(optimizer::Optimizer)
-    return optimizer.poly = nothing
+    optimizer.poly = nothing
+    optimizer.solve_time = NaN
+    return
 end
+
 MOI.is_empty(optimizer::Optimizer) = optimizer.poly === nothing
+
+MOI.supports_incremental_interface(::Optimizer) = true
+
+function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
+    return MOI.Utilities.default_copy_to(dest, src)
+end
+
 function MOI.optimize!(optimizer::Optimizer)
     options =
         Dict{Symbol,Any}(Symbol(key) => val for (key, val) in optimizer.options)
+    start_time = time()
     ret = sos_decomp(optimizer.poly; options...)
+    optimizer.solve_time = time() - start_time
     optimizer.fvals = ret.fvals
     optimizer.Uopt = ret.Uopt
     optimizer.ret = ret.ret
     return
 end
 
+MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec) = optimizer.solve_time
+
 function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
     # TODO if we use the MOI API of NLopt, we can just redirect
     #      the call to it