Update to MP renamings (#31)

* Update to MP renamings * Use MP#bl/rename * Fixes * Fixes * Fix ci * Remove Sequences * Updates * Fixes * Fix doc script * Checkout v3 * Fix * Remove COI * Update symmetric_group * Fixes * Fix * Fixes * rewrite the definitions of groups in symmetry.jl (#32) * rewrite the definitions of groups in symmetry.jl * add Random * Update src/symmetry.jl Co-authored-by: Marek Kaluba <marek.kaluba@kit.edu> * Fixes * fixes after the rewrite of symmetry.jl (#33) * sort out deepcopy_internal stuff * move imports/exports away from symmetry.jl * update action.jl to the rewrite of symmetry.jl * Fixes * Use latest SOS * Update GH action --------- Co-authored-by: Marek Kaluba <marek.kaluba@kit.edu>
blegat · Sep 26, 2023 · 19bf2eb · 19bf2eb
1 parent 187c7c4
commit 19bf2eb
Show file tree

Hide file tree

Showing 16 changed files with 301 additions and 377 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -22,7 +22,7 @@ jobs:
             os: ubuntu-latest
             arch: x86
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
@@ -42,6 +42,6 @@ jobs:
         with:
           depwarn: error
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v3
         with:
           file: lcov.info
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
@@ -9,11 +9,19 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@latest
         with:
           # Build documentation on the latest Julia 1.x
           version: '1'
+      - name: Install dependencies
+        shell: julia --project=docs/ {0}
+        run: |
+          using Pkg
+          Pkg.add([
+              PackageSpec(path=pwd()),
+          ])
+          Pkg.instantiate()
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.instantiate();
                                        Pkg.develop(PackageSpec(path=pwd()))'

diff --git a/Project.toml b/Project.toml
@@ -10,17 +10,18 @@ GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 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.4"
+MultivariatePolynomials = "0.5"
 MutableArithmetics = "1"
 Reexport = "1"
-SumOfSquares = "0.6"
+SumOfSquares = "0.7"
 julia = "1.6"
 
 [extras]

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -1,6 +1,8 @@
 [deps]
 CSDP = "0a46da34-8e4b-519e-b418-48813639ff34"
+Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+CondensedMatterSOS = "676dea7e-7535-42cf-a5c4-14732a6b7aba"
 Cyclotomics = "da8f5974-afbb-4dc8-91d8-516d5257c83b"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

diff --git a/examples/Benchmark_1.jl b/examples/Benchmark_1.jl
@@ -7,26 +7,22 @@
 
 using Test #src
 using CondensedMatterSOS
-import MultivariatePolynomials
-const MP = MultivariatePolynomials
+import MultivariatePolynomials as MP
 @spin σ[1:3]
 heisenberg_hamiltonian(σ, true)
 
 ## Let's pick a solver from [this list](https://jump.dev/JuMP.jl/dev/installation/#Getting-Solvers).
 
-using CSDP
-solver = optimizer_with_attributes(
-    () -> MOIU.CachingOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()), CSDP.Optimizer()),
-    MOI.Silent() => false,
-);
+import Clarabel
+solver = Clarabel.Optimizer
 
 # We can compute a lower bound `-2√2` to the ground state energy as follow:
 
 function hamiltonian_energy(N, maxdegree, solver; symmetry=true, consecutive=false, kws...)
     @spin σ[1:N]
     H = heisenberg_hamiltonian(σ, true)
     G = Lattice1Group(N)
-    cone = NonnegPolyInnerCone{SumOfSquares.COI.HermitianPositiveSemidefiniteConeTriangle}()
+    cone = NonnegPolyInnerCone{MOI.HermitianPositiveSemidefiniteConeTriangle}()
     @assert iseven(maxdegree)
     cert = SumOfSquares.Certificate.FixedBasis(
         cone,
@@ -74,10 +70,10 @@ bound
 # The reduction is obtained by block diagonalizing with a change of polynomial
 # basis to the isotypical basis.
 
-display([M.basis.polynomials for M in ν.sub_moment_matrices])
+display([M.basis.polynomials for M in ν.blocks])
 
-@test length(ν.sub_moment_matrices) == 7 #src
-[M.basis.polynomials for M in ν.sub_moment_matrices]
+@test length(ν.blocks) == 7 #src
+[M.basis.polynomials for M in ν.blocks]
 
 # Let's try this for 3 sites. First without symmetry.
 
@@ -102,7 +98,7 @@ bound, gram, ν = hamiltonian_energy(
 
 # Let's look at the isotypical basis.
 
-display([M.basis.polynomials for M in ν.sub_moment_matrices])
+display([M.basis.polynomials for M in ν.blocks])
 
 # Now let's define a function for our common use case.
 
@@ -118,7 +114,7 @@ function f(L, d=1, consecutive=false; symmetry=true)
         symmetry=symmetry,
     )
     @show bound
-    for M in ν.sub_moment_matrices
+    for M in ν.blocks
         display(round.(M.basis.polynomials, digits=6))
     end
     println("E/N = ", bound / L)

diff --git a/examples/Glass.jl b/examples/Glass.jl
@@ -19,11 +19,8 @@ hamiltonian(σ)
 
 # Let's pick a solver from [this list](https://jump.dev/JuMP.jl/dev/installation/#Getting-Solvers).
 
-using CSDP
-solver = optimizer_with_attributes(
-    () -> MOIU.CachingOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()), CSDP.Optimizer()),
-    MOI.Silent() => true
-)
+import Clarabel
+solver = Clarabel.Optimizer
 
 # We can compute a lower bound `-0.8031` to the ground state energy as follow:
 
@@ -49,8 +46,8 @@ bound
 
 # But with a smaller basis:
 
-@test length(ν.sub_moment_matrices) == 2 #src
-[M.basis.monomials for M in ν.sub_moment_matrices]
+@test length(ν.blocks) == 2 #src
+[M.basis.monomials for M in ν.blocks]
 
 # Using term sparsity with chordal completion, we get a smaller bound:
 
@@ -60,8 +57,8 @@ bound
 
 # But with an even smaller basis:
 
-@test length(ν.sub_moment_matrices) == 5 #src
-[M.basis.monomials for M in ν.sub_moment_matrices]
+@test length(ν.blocks) == 5 #src
+[M.basis.monomials for M in ν.blocks]
 
 # ## 2D
 
@@ -94,16 +91,16 @@ bound
 
 # But with a smaller basis:
 
-@test length(ν.sub_moment_matrices) == 2 #src
-[M.basis.monomials for M in ν.sub_moment_matrices]
+@test length(ν.blocks) == 2 #src
+[M.basis.monomials for M in ν.blocks]
 
 # Using term sparsity with chordal completion, we get a smaller bound:
 
 bound, gram, ν = hamiltonian_energy(2, 2, 2, solver, sparsity = MonomialSparsity(ChordalCompletion()))
-@test bound ≈ -5.1878 rtol=1e-3 #src
+@test bound ≈ -5.4659 rtol=1e-3 #src
 bound
 
 # But with an even smaller basis:
 
-@test length(ν.sub_moment_matrices) == 9 #src
-[M.basis.monomials for M in ν.sub_moment_matrices]
+@test length(ν.blocks) == 9 #src
+[M.basis.monomials for M in ν.blocks]
diff --git a/examples/Ising.jl b/examples/Ising.jl
@@ -42,8 +42,8 @@ bound
 
 # But with a smaller basis:
 
-@test length(ν.sub_moment_matrices) == 2 #src
-[M.basis.monomials for M in ν.sub_moment_matrices]
+@test length(ν.blocks) == 2 #src
+[M.basis.monomials for M in ν.blocks]
 
 # Using term sparsity with chordal completion, we get a smaller bound:
 
@@ -53,5 +53,5 @@ bound
 
 # But with an even smaller basis:
 
-@test length(ν.sub_moment_matrices) == 5 #src
-[M.basis.monomials for M in ν.sub_moment_matrices]
+@test length(ν.blocks) == 5 #src
+[M.basis.monomials for M in ν.blocks]
diff --git a/src/CondensedMatterSOS.jl b/src/CondensedMatterSOS.jl
@@ -35,15 +35,14 @@ function MP.name_base_indices(var::SpinVariable) # Used to print variable
     end
 end
 
-include("sequences.jl")
-import .Sequences
-
 include("operators.jl")
 include("monom-set.jl")
 
 include("ising.jl")
 
 include("symmetry.jl")
+export Lattice1Group
+
 include("action.jl")
 include("sos.jl")
 

diff --git a/src/action.jl b/src/action.jl
@@ -1,3 +1,4 @@
+import SumOfSquares.Symmetry
 export Action
 
 struct Action <: Symmetry.OnMonomials
@@ -11,7 +12,7 @@ function Symmetry.SymbolicWedderburn.action(a::Action, el::DirectSum, mono::Cond
         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
+        id = ((rel_id + el.h.id) % el.h.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

diff --git a/src/monom-set.jl b/src/monom-set.jl
@@ -31,12 +31,12 @@ end
 function _monomials(sites::Vector{Int}, deg::Int, filter::Function, monos, consecutive)
     if consecutive
         if deg == 0
-            _add_monomials!(monos, Int[], MP.constantmonomial(SpinMonomial), 1, filter)
+            _add_monomials!(monos, Int[], MP.constant_monomial(SpinMonomial), 1, filter)
             return monos
         end
         if deg >= length(sites)
             if deg == length(sites)
-                _add_monomials!(monos, sites, MP.constantmonomial(SpinMonomial), 1, filter)
+                _add_monomials!(monos, sites, MP.constant_monomial(SpinMonomial), 1, filter)
             end
             return monos
         end
@@ -48,11 +48,11 @@ function _monomials(sites::Vector{Int}, deg::Int, filter::Function, monos, conse
                 append!(active_sites, sites[1:(n-m)])
             end
             @assert length(active_sites) == deg
-            _add_monomials!(monos, active_sites, MP.constantmonomial(SpinMonomial), 1, filter)
+            _add_monomials!(monos, active_sites, MP.constant_monomial(SpinMonomial), 1, filter)
         end
     else
         for active_sites in combinations(sites, deg)
-            _add_monomials!(monos, active_sites, MP.constantmonomial(SpinMonomial), 1, filter)
+            _add_monomials!(monos, active_sites, MP.constant_monomial(SpinMonomial), 1, filter)
         end
     end
     return monos

diff --git a/src/operators.jl b/src/operators.jl
@@ -1,3 +1,7 @@
+function MA.promote_operation(::typeof(*), ::Type{<:SpinLike}, ::Type{<:SpinLike})
+    return MP.term_type(SpinMonomial, Complex{Int})
+end
+
 function _coef_variable(a::SpinVariable, b::SpinVariable)
     i = a.index
     j = b.index
@@ -13,7 +17,7 @@ end
 function Base.:*(a::SpinVariable, b::SpinVariable)
     if a.id == b.id
         if a.index == b.index
-            return (1 + 0im) * MP.constantmonomial(a)
+            return (1 + 0im) * MP.constant_monomial(a)
             # We want to return `1` but in which type ?
             # We could use `Bool` type as it the type compatible with the most other types in Julia
             # but currently the convention is `Int` in MP, i.e. variables and monomials are `AbstractTerm{Int}`.
@@ -62,21 +66,6 @@ function Base.:*(a::SpinMonomial, b::SpinMonomial)
     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)
     return (a.id==b.id) && (a.index==b.index)
 end
@@ -122,28 +111,12 @@ function Base.isless(a::SpinMonomial, b::SpinMonomial)
     end
     return false
 end
-
-combine_plus(t1::SpinTerm, t2::SpinTerm) = (MP.coefficient(t1) + MP.coefficient(t2)) * MP.monomial(t1)
-combine_plus(t::SpinTerm, ::MA.Zero) = MA.mutable_copy(t)
-combine_plus(::MA.Zero, t::SpinTerm) = MA.mutable_copy(t)
-combine_minus(t1::SpinTerm, t2::SpinTerm) = (MP.coefficient(t1) - MP.coefficient(t2)) * MP.monomial(t1)
-combine_minus(t::SpinTerm, ::MA.Zero) = MA.mutable_copy(t)
-combine_minus(::MA.Zero, t::SpinTerm) = -t
-function MA.promote_operation(::typeof(combine_plus), ::Type{SpinTerm{S}}, ::Type{SpinTerm{T}}) where {S, T}
-    return SpinTerm{MA.promote_operation(+, S, T)}
-end
-function MA.promote_operation(::typeof(combine_minus), ::Type{SpinTerm{S}}, ::Type{SpinTerm{T}}) where {S, T}
-    return SpinTerm{MA.promote_operation(-, S, T)}
+function MP.compare(a::SpinMonomial, b::SpinMonomial)
+    if a == b
+        return 0
+    elseif a < b
+        return -1
+    else
+        return 1
+    end
 end
-compare(t1::SpinTerm, t2::SpinTerm) = monomial(t1) > monomial(t2)
-
-# TODO taken from TypedPolynomials
-join_terms_plus(terms1::AbstractArray{<:SpinTerm}, terms2::AbstractArray{<:SpinTerm}) = Sequences.mergesorted(terms1, terms2, compare, combine_plus)
-join_terms_minus(terms1::AbstractArray{<:SpinTerm}, terms2::AbstractArray{<:SpinTerm}) = Sequences.mergesorted(terms1, terms2, compare, combine_minus)
-
-Base.:+(p1::SpinPolynomial, p2::SpinPolynomial) = SpinPolynomial(join_terms_plus(MP.terms(p1), MP.terms(p2)))
-Base.:+(p::SpinPolynomial, t::SpinTerm) = SpinPolynomial(join_terms_plus(MP.terms(p), [t]))
-Base.:+(t::SpinTerm, p::SpinPolynomial) = p + t
-Base.:-(p1::SpinPolynomial, p2::SpinPolynomial) = SpinPolynomial(join_terms_minus(MP.terms(p1), MP.terms(p2)))
-Base.:-(p::SpinPolynomial, t::SpinTerm) = SpinPolynomial(join_terms_minus(MP.terms(p), [t]))
-Base.:-(t::SpinTerm, p::SpinPolynomial) = SpinPolynomial(join_terms_minus([t], MP.terms(p)))