Update to SumOfSquares v0.6

.github/workflows/ci.yml

@@ -12,7 +12,7 @@ jobs:
       fail-fast: false
       matrix:
         include:
-          - version: '1.0'
+          - version: '1.6'
             os: ubuntu-latest
             arch: x64
           - version: '1'
@@ -21,9 +21,6 @@ jobs:
           - version: '1'
             os: ubuntu-latest
             arch: x86
-          - version: 'nightly'
-            os: ubuntu-latest
-            arch: x64
     steps:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1

Project.toml

@@ -6,19 +6,22 @@ version = "0.1.0"
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 Combinatorics = "1.0.2"
 DataStructures = "0.18.7"
-MultivariatePolynomials = "0.3.11"
-MutableArithmetics = "0.2.10, 0.3"
-Reexport = "0.2, 1"
-SumOfSquares = "0.4.3, 0.5"
+MultivariatePolynomials = "0.4"
+MutableArithmetics = "1"
+Reexport = "1"
+SumOfSquares = "0.6"
+julia = "1.6"
 
 [extras]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/Project.toml

@@ -14,4 +14,4 @@ SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-Documenter = "~0.26"
+Documenter = "~0.27"

examples/Benchmark_1.jl

@@ -7,6 +7,8 @@
 
 using Test #src
 using CondensedMatterSOS
+import MultivariatePolynomials
+const MP = MultivariatePolynomials
 @spin σ[1:3]
 heisenberg_hamiltonian(σ, true)
 
@@ -16,36 +18,12 @@ using CSDP
 solver = optimizer_with_attributes(
     () -> MOIU.CachingOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()), CSDP.Optimizer()),
     MOI.Silent() => false,
-)
+);
 
 # We can compute a lower bound `-2√2` to the ground state energy as follow:
 
-include(joinpath(dirname(dirname(pathof(CondensedMatterSOS))), "examples", "symmetry.jl"))
 function hamiltonian_energy(N, maxdegree, solver; symmetry=true, consecutive=false, kws...)
     @spin σ[1:N]
-    function action(mono::CondensedMatterSOS.SpinMonomial, el::DirectSum)
-        isempty(mono.variables) && return 1 * mono
-        sign = 1
-        vars = map(values(mono.variables)) do var
-            rel_id = var.id - σ[1][1].id
-            rel_index = var.index + 1
-            @assert σ[rel_index][rel_id + 1] == var
-            id = ((rel_id + el.c.id) % el.c.n) + σ[1][1].id
-            index = (rel_index^el.k.p) - 1
-            new_var = CondensedMatterSOS.SpinVariable(id, index)
-            if el.k.k.id != 0 && el.k.k.id != index + 1
-                sign *= -1
-            end
-            return new_var
-        end
-        return sign * CondensedMatterSOS.SpinMonomial(vars)
-    end
-    function action(term::CondensedMatterSOS.SpinTerm, el::DirectSum)
-        return MP.coefficient(term) * action(MP.monomial(term), el)
-    end
-    function action(poly::CondensedMatterSOS.SpinPolynomial, el::DirectSum)
-        return MP.polynomial([action(term, el) for term in MP.terms(poly)])
-    end
     H = heisenberg_hamiltonian(σ, true)
     G = Lattice1Group(N)
     cone = NonnegPolyInnerCone{SumOfSquares.COI.HermitianPositiveSemidefiniteConeTriangle}()
@@ -54,11 +32,9 @@ function hamiltonian_energy(N, maxdegree, solver; symmetry=true, consecutive=fal
         cone,
         MonomialBasis(MP.monomials(σ[1], 0:div(maxdegree, 2), consecutive=consecutive)),
     )
-    display(cert.basis.monomials)
-    certificate = SymmetricIdeal(
+    certificate = Symmetry.Ideal(
+        Symmetry.Pattern(G, Action(σ)),
         cert,
-        G,
-        action,
     )
     if symmetry
         energy(H, maxdegree, solver; certificate = certificate, kws...)
@@ -71,7 +47,7 @@ bound, gram, ν = hamiltonian_energy(
     2,
     solver,
     symmetry = false,
-    sparsity = NoSparsity(),
+    sparsity = SumOfSquares.Sparsity.NoPattern(),
 )
 @test bound ≈ -6 rtol=1e-6 #src
 bound

examples/Glass.jl

@@ -25,15 +25,15 @@ solver = optimizer_with_attributes(
     MOI.Silent() => true
 )
 
-# We can compute a lower bound `-0.853452` to the ground state energy as follow:
+# We can compute a lower bound `-0.8031` to the ground state energy as follow:
 
 function hamiltonian_energy(N, maxdegree, solver; kws...)
     @spin σ[1:N]
     H = hamiltonian(σ)
     energy(H, maxdegree, solver; kws...)
 end
 bound, gram, ν = hamiltonian_energy(2, 2, solver, sparsity = NoSparsity())
-@test bound ≈ -0.853452 rtol=1e-5 #src
+@test bound ≈ -0.8031 rtol=1e-3 #src
 bound
 
 # We can see that the moment matrix uses all monomials:
@@ -44,7 +44,7 @@ bound
 # Using term sparsity with block/cluster completion, we get the same bound:
 
 bound, gram, ν = hamiltonian_energy(2, 2, solver)
-@test bound ≈ -0.853452 rtol=1e-5 #src
+@test bound ≈ -0.8031 rtol=1e-3 #src
 bound
 
 # But with a smaller basis:
@@ -55,7 +55,7 @@ bound
 # Using term sparsity with chordal completion, we get a smaller bound:
 
 bound, gram, ν = hamiltonian_energy(2, 2, solver, sparsity = MonomialSparsity(ChordalCompletion()))
-@test bound ≈ -1.097288 rtol=1e-5 #src
+@test bound ≈ -0.87058 rtol=1e-3 #src
 bound
 
 # But with an even smaller basis:
@@ -70,15 +70,15 @@ bound
 @spin σ[1:2, 1:2]
 hamiltonian(σ)
 
-# We can compute a lower bound `-1.990513` to the ground state energy as follow:
+# We can compute a lower bound `-4.15244` to the ground state energy as follow:
 
 function hamiltonian_energy(N, M, maxdegree, solver; kws...)
     @spin σ[1:N, 1:M]
     H = hamiltonian(σ)
     energy(H, maxdegree, solver; kws...)
 end
 bound, gram, ν = hamiltonian_energy(2, 2, 2, solver, sparsity = NoSparsity())
-@test bound ≈ -1.990513 rtol=1e-5 #src
+@test bound ≈ -4.15244 rtol=1e-3 #src
 bound
 
 # We can see that the moment matrix uses all monomials:
@@ -89,7 +89,7 @@ bound
 # Using term sparsity with block/cluster completion, we get the same bound:
 
 bound, gram, ν = hamiltonian_energy(2, 2, 2, solver)
-@test bound ≈ -1.990513 rtol=1e-5 #src
+@test bound ≈ -4.15244 rtol=1e-3 #src
 bound
 
 # But with a smaller basis:
@@ -100,7 +100,7 @@ bound
 # Using term sparsity with chordal completion, we get a smaller bound:
 
 bound, gram, ν = hamiltonian_energy(2, 2, 2, solver, sparsity = MonomialSparsity(ChordalCompletion()))
-@test bound ≈ -2.655704 rtol=1e-5 #src
+@test bound ≈ -5.1878 rtol=1e-3 #src
 bound
 
 # But with an even smaller basis:

src/CondensedMatterSOS.jl

@@ -43,6 +43,8 @@ include("monom-set.jl")
 
 include("ising.jl")
 
+include("symmetry.jl")
+include("action.jl")
 include("sos.jl")
 
 export @spin

src/action.jl

@@ -0,0 +1,23 @@
+export Action
+
+struct Action <: Symmetry.OnMonomials
+    σ
+end
+Symmetry.SymbolicWedderburn.coeff_type(::Action) = Float64
+function Symmetry.SymbolicWedderburn.action(a::Action, el::DirectSum, mono::CondensedMatterSOS.SpinMonomial)
+    isempty(mono.variables) && return 1 * mono
+    sign = 1
+    vars = map(values(mono.variables)) do var
+        rel_id = var.id - a.σ[1][1].id
+        rel_index = var.index + 1
+        @assert a.σ[rel_index][rel_id + 1] == var
+        id = ((rel_id + el.c.id) % el.c.n) + a.σ[1][1].id
+        index = (rel_index^el.k.p) - 1
+        new_var = CondensedMatterSOS.SpinVariable(id, index)
+        if el.k.k.id != 0 && el.k.k.id != index + 1
+            sign *= -1
+        end
+        return new_var
+    end
+    return sign * CondensedMatterSOS.SpinMonomial(vars)
+end

src/operators.jl

@@ -61,6 +61,20 @@ function Base.:*(a::SpinMonomial, b::SpinMonomial)
     end
     return coef * c
 end
+
+function Base.:*(a::SpinTerm, b::Union{SpinMonomial, SpinVariable})
+    t = MP.monomial(a) * b
+    return MP.term(MP.coefficient(a) * MP.coefficient(t), MP.monomial(t))
+end
+function Base.:*(a::Union{SpinMonomial, SpinVariable}, b::SpinTerm)
+    t = a * MP.monomial(b)
+    return MP.term(MP.coefficient(t) * MP.coefficient(b), MP.monomial(t))
+end
+function Base.:*(a::SpinTerm, b::SpinTerm)
+    t = a * MP.monomial(b)
+    return MP.term(MP.coefficient(t) * MP.coefficient(b), MP.monomial(t))
+end
+
 MP.multconstant(α, m::SpinMonomial) = SpinTerm(α, m)
 MP.multconstant(m::SpinMonomial, α) = SpinTerm(α, m)
 function Base.:(==)(a::SpinVariable, b::SpinVariable)

src/sos.jl

@@ -1,4 +1,5 @@
-using SumOfSquares
+using Reexport
+@reexport using SumOfSquares
 export optimizer_with_attributes, MOI, MOIU, NoSparsity, MonomialSparsity, ChordalCompletion
 
 export energy
@@ -22,9 +23,9 @@ of the SumOfSquares package.
 """
 function energy(H, maxdegree, solver;
     cone=NonnegPolyInnerCone{SumOfSquares.COI.HermitianPositiveSemidefiniteConeTriangle}(),
-    sparsity=MonomialSparsity(),
+    sparsity=Sparsity.Monomial(),
     non_sparse=SumOfSquares.Certificate.MaxDegree(cone, MonomialBasis, maxdegree),
-    certificate=sparsity isa NoSparsity ? non_sparse : SumOfSquares.Certificate.SparseIdeal(sparsity, non_sparse),
+    certificate=sparsity isa Sparsity.NoPattern ? non_sparse : SumOfSquares.Certificate.Sparsity.Ideal(sparsity, non_sparse),
     kws...
 )
     model = MOI.instantiate(solver, with_bridge_type=Float64)
@@ -34,11 +35,11 @@ function energy(H, maxdegree, solver;
     MOI.Bridges.add_bridge(model, PolyJuMP.ZeroPolynomialBridge{Complex{Float64}})
     MOI.Bridges.add_bridge(model, SumOfSquares.COI.Bridges.Variable.HermitianToSymmetricPSDBridge{Float64})
     MOI.Bridges.add_bridge(model, SumOfSquares.COI.Bridges.Constraint.SplitZeroBridge{Float64})
-    γ = MOI.SingleVariable(MOI.add_variable(model))
+    γ = MOI.add_variable(model)
     poly = convert(SpinPolynomial{Complex{Float64}}, H) - (1.0 + 0.0im) * γ
     c = SumOfSquares.add_constraint(model, poly, SOSCone(); ideal_certificate=certificate, kws...)
     MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), γ)
+    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), γ)
     MOI.optimize!(model)
     if MOI.get(model, MOI.TerminationStatus()) != MOI.OPTIMAL
         @warn("Termination status: $(MOI.get(model, MOI.TerminationStatus())), $(MOI.get(model, MOI.RawStatusString()))")

examples/symmetry.jl → src/symmetry.jl

@@ -1,6 +1,8 @@
 using SumOfSquares
-include(joinpath(dirname(dirname(pathof(SumOfSquares))), "examples", "symmetry.jl"))
-include(joinpath(dirname(dirname(pathof(SumOfSquares))), "examples", "scaled_perm.jl"))
+using GroupsCore
+using PermutationGroups
+
+export Lattice1Group
 
 struct KleinElement <: GroupElement
     id::Int
@@ -180,6 +182,7 @@ Base.:^(a::DirectSum, b::DirectSum) = conj(a, b)
 struct Lattice1Group <: Group
     n::Int
 end
+Base.eltype(::Lattice1Group) = DirectSum
 Base.one(el::DirectSum) = DirectSum(one(el.c), one(el.k))
 Base.one(L::Lattice1Group) = DirectSum(CyclicElem(L.n, 0), one(KleinPermGroup()))
 function PermutationGroups.gens(L::Lattice1Group)

src/types.jl

@@ -53,6 +53,10 @@ struct SpinTerm{T} <: MP.AbstractTerm{T}
     monomial::SpinMonomial
 end
 
+function MP.term(coefficient, monomial::SpinMonomial)
+    return SpinTerm(coefficient, monomial)
+end
+
 # TODO move to MP
 MP.convertconstant(::Type{SpinTerm{T}}, α) where {T} = convert(T, α) * MP.constantmonomial(SpinTerm{T})
 Base.copy(t::SpinTerm) = SpinTerm(MP.coefficient(t), copy(MP.monomial(t)))

test/runtests.jl

@@ -1,4 +1,5 @@
 using MultivariatePolynomials
+const MP = MultivariatePolynomials
 using Test
 using CondensedMatterSOS
 
@@ -154,10 +155,10 @@ end
     @test monovec([1, sigmax[1], sigmax[2]]) == [sigmax[1], sigmax[2], 1]
     @test CMS.monomials([sigmax[1], sigmax[2]], 0) == [1]
     for consecutive in [false, true]
-        @test monomials([sigmax[1], sigmax[2]], 0, consecutive=consecutive) == [1]
-        @test all(monomials([sigmax[1], sigmax[2]], 2, consecutive=consecutive) .== [CMS.SpinVariable(sigmax[1].id, i) * CMS.SpinVariable(sigmax[2].id, j) for i in 0:2 for j in 0:2])
+        @test MP.monomials([sigmax[1], sigmax[2]], 0, consecutive=consecutive) == [1]
+        @test all(MP.monomials([sigmax[1], sigmax[2]], 2, consecutive=consecutive) .== [CMS.SpinVariable(sigmax[1].id, i) * CMS.SpinVariable(sigmax[2].id, j) for i in 0:2 for j in 0:2])
     end
-    @test all(monomials([sigmax[1], sigmax[2], sigmax[4]], 2, consecutive=true) .== append!(
+    @test all(MP.monomials([sigmax[1], sigmax[2], sigmax[4]], 2, consecutive=true) .== append!(
         append!(
             [CMS.SpinVariable(sigmax[1].id, i) * CMS.SpinVariable(sigmax[2].id, j) for i in 0:2 for j in 0:2],
             [CMS.SpinVariable(sigmax[2].id, i) * CMS.SpinVariable(sigmax[4].id, j) for i in 0:2 for j in 0:2]