Merge pull request #69 from kalmarek/mk/schreier_hom

add Schreier evaluation of extension homomorphisms

Project.toml

@@ -20,7 +20,7 @@ AbstractPermutations = "0.3"
 Cyclotomics = "0.3"
 GroupsCore = "0.5"
 DynamicPolynomials = "0.6"
-PermutationGroups = "0.6.2"
+PermutationGroups = "0.6.3"
 PrecompileTools = "1"
 PrettyTables = "2"
 Primes = "0.4, 0.5"

examples/action_polynomials.jl

@@ -1,8 +1,5 @@
 using DynamicPolynomials
-const DP = DynamicPolynomials
-using GroupsCore
-
-const SW = SymbolicWedderburn
+import DynamicPolynomials as DP
 
 # Defining action on polynomials by acting on terms and monomials:
 function SW.action(a::SW.Action, el::GroupElement, poly::DP.AbstractPolynomial)

examples/dihedral.jl

@@ -88,3 +88,6 @@ function GroupsCore.order(T::Type, el::DihedralElement)
     iszero(el.id) && return T(1)
     return T(div(el.n, gcd(el.n, el.id)))
 end
+
+# this is needed for using them in StarAlgebra:
+SA.comparable(::Type{DihedralElement}) = (a, b) -> hash(a) < hash(b)

examples/ex_C2_linear.jl

@@ -40,10 +40,7 @@ using Test
 @testset "Decompose in basis" begin
     g = gens(G, 1)
 
-    basis = StarAlgebras.FixedBasis(
-        monomials([x, y], 0:4),
-        StarAlgebras.DiracMStructure(*),
-    )
+    basis = SA.FixedBasis(monomials([x, y], 0:4), SA.DiracMStructure(*))
     k = SymbolicWedderburn.action(By90Rotation(), g, basis[2])
     ehom = SymbolicWedderburn.ExtensionHomomorphism(By90Rotation(), basis)
     idcs, vals = SymbolicWedderburn.decompose(k, ehom)
@@ -60,10 +57,8 @@ end
 
 @testset "induced Matrix Representation" begin
     g = gens(G, 1)
-    monomial_basis = StarAlgebras.FixedBasis(
-        monomials([x, y], 0:4),
-        StarAlgebras.DiracMStructure(*),
-    )
+    monomial_basis =
+        SA.FixedBasis(monomials([x, y], 0:4), SA.DiracMStructure(*))
     ehom =
         SymbolicWedderburn.ExtensionHomomorphism(By90Rotation(), monomial_basis)
     m = droptol!(SymbolicWedderburn.induce(By90Rotation(), ehom, g), 1e-15)

examples/ex_S4.jl

@@ -1,12 +1,3 @@
-using SymbolicWedderburn
-using PermutationGroups
-
-using DynamicPolynomials
-
-include(joinpath(@__DIR__, "action_polynomials.jl"))
-include(joinpath(@__DIR__, "sos_problem.jl"))
-include(joinpath(@__DIR__, "solver.jl"))
-
 const N = 4
 @polyvar x[1:N]
 

examples/ex_robinson_form.jl

@@ -11,7 +11,7 @@ using LinearAlgebra
 using SparseArrays
 
 using SymbolicWedderburn
-using SymbolicWedderburn.StarAlgebras
+import SymbolicWedderburn.SA as SA
 using DynamicPolynomials
 
 # definitions of general actions on polynomials:

examples/run_examples.jl

@@ -1,7 +1,24 @@
 using Pkg
 Pkg.activate(@__DIR__)
 Pkg.instantiate()
+
+module Examples
+using LinearAlgebra
+using SparseArrays
+
+using GroupsCore
+import SymbolicWedderburn as SW
+import SymbolicWedderburn.StarAlgebras as SA
+using PermutationGroups
+import PermutationGroups.AP as AP
+
+include(joinpath(@__DIR__, "action_polynomials.jl"))
+include(joinpath(@__DIR__, "sos_problem.jl"))
+include(joinpath(@__DIR__, "solver.jl"))
+
 include(joinpath(@__DIR__, "ex_C2_linear.jl"))
 include(joinpath(@__DIR__, "ex_S4.jl"))
 include(joinpath(@__DIR__, "ex_motzkin.jl"))
 include(joinpath(@__DIR__, "ex_robinson_form.jl"))
+
+end

examples/solver.jl

@@ -1,4 +1,3 @@
-using JuMP
 import SCS
 
 function scs_optimizer(;

examples/sos_problem.jl

@@ -1,16 +1,7 @@
-using LinearAlgebra
-using SparseArrays
-
-using GroupsCore
-
-using DynamicPolynomials
-using SymbolicWedderburn
-import SymbolicWedderburn.StarAlgebras
-
 using JuMP
 
-function DynamicPolynomials.coefficients(p, b::StarAlgebras.FixedBasis)
-    return DynamicPolynomials.coefficients(p, b.elts)
+function DP.coefficients(p, b::SA.FixedBasis)
+    return DP.coefficients(p, b.elts)
 end
 
 function invariant_constraint!(
@@ -28,19 +19,13 @@ function invariant_constraint!(
 end
 
 function sos_problem(poly::AbstractPolynomial)
-    vars = DynamicPolynomials.variables(poly)
+    vars = DP.variables(poly)
 
-    basis_psd = DynamicPolynomials.monomials(
-        vars,
-        0:DynamicPolynomials.maxdegree(poly)÷2,
-    )
+    basis_psd = DP.monomials(vars, 0:DP.maxdegree(poly)÷2)
 
-    basis_constraints = StarAlgebras.FixedBasis(
-        DynamicPolynomials.monomials(
-            vars,
-            0:DynamicPolynomials.maxdegree(poly),
-        ),
-        StarAlgebras.DiracMStructure(*),
+    basis_constraints = SA.FixedBasis(
+        DP.monomials(vars, 0:DP.maxdegree(poly)),
+        SA.DiracMStructure(*),
     )
 
     M = [basis_constraints[x*y] for x in basis_psd, y in basis_psd]
@@ -54,7 +39,7 @@ function sos_problem(poly::AbstractPolynomial)
 
     objective = poly - t
     for (idx, b) in enumerate(basis_constraints)
-        c = DynamicPolynomials.coefficient(objective, b)
+        c = DP.coefficient(objective, b)
         JuMP.@constraint sos_model LinearAlgebra.dot(P, M .== idx) == c
     end
 
@@ -64,7 +49,7 @@ end
 function sos_problem(
     poly::AbstractPolynomial,
     invariant_vs::AbstractVector,
-    basis_constraints::StarAlgebras.AbstractBasis,
+    basis_constraints::SA.AbstractBasis,
     basis_psd,
     T=Float64,
 )
@@ -80,7 +65,7 @@ function sos_problem(
     # preallocating
     M_orb = similar(M, T)
 
-    C = DynamicPolynomials.coefficients(poly - t, basis_constraints)
+    C = DP.coefficients(poly - t, basis_constraints)
 
     for iv in invariant_vs
         c = dot(C, iv)
@@ -94,18 +79,18 @@ end
 
 function sos_problem(
     poly::AbstractPolynomial,
-    wedderburn::SymbolicWedderburn.WedderburnDecomposition,
+    wedderburn::SW.WedderburnDecomposition,
     basis_psd;
 )
     m = JuMP.Model()
 
-    M = let basis_constraints = SymbolicWedderburn.basis(wedderburn)
+    M = let basis_constraints = SA.basis(wedderburn)
         [basis_constraints[x*y] for x in basis_psd, y in basis_psd]
     end
 
     JuMP.@variable m t
     JuMP.@objective m Max t
-    psds = map(SymbolicWedderburn.direct_summands(wedderburn)) do ds
+    psds = map(SW.direct_summands(wedderburn)) do ds
         dim = size(ds, 1)
         P = JuMP.@variable m [1:dim, 1:dim] Symmetric
         JuMP.@constraint m P in PSDCone()
@@ -116,15 +101,13 @@ function sos_problem(
     # Mπs = zeros.(eltype(wedderburn), size.(psds))
     M_orb = similar(M, eltype(wedderburn))
 
-    C = DynamicPolynomials.coefficients(
-        poly - t,
-        SymbolicWedderburn.basis(wedderburn),
-    )
+    C = DP.coefficients(poly - t, SW.basis(wedderburn))
+
     for iv in invariant_vectors(wedderburn)
         c = dot(C, iv)
         M_orb = invariant_constraint!(M_orb, M, iv)
-        # Mπs = SymbolicWedderburn.diagonalize!(Mπs, M_orb, wedderburn)
-        Mπs = SymbolicWedderburn.diagonalize(M_orb, wedderburn)
+        # Mπs = SW.diagonalize!(Mπs, M_orb, wedderburn)
+        Mπs = SW.diagonalize(M_orb, wedderburn)
 
         JuMP.@constraint m sum(
             dot(Mπ, Pπ) for (Mπ, Pπ) in zip(Mπs, psds) if !iszero(Mπ)
@@ -136,18 +119,19 @@ end
 function sos_problem(
     poly::AbstractPolynomial,
     G::Group,
-    action::SymbolicWedderburn.Action,
+    action::SW.Action,
     T=Float64;
     decompose_psd=true,
     semisimple=false
 )
-    max_deg = DynamicPolynomials.maxdegree(poly)
-    vars = DynamicPolynomials.variables(poly)
-    basis_psd = DynamicPolynomials.monomials(vars, 0:max_deg÷2)
-    basis_constraints = DynamicPolynomials.monomials(vars, 0:max_deg)
+    max_deg = DP.maxdegree(poly)
+    vars = DP.variables(poly)
+    basis_psd = DP.monomials(vars, 0:max_deg÷2)
+    basis_constraints = DP.monomials(vars, 0:max_deg)
 
     if decompose_psd == true
-        wedderburn, symmetry_adaptation_time = @timed WedderburnDecomposition(
+        wedderburn, symmetry_adaptation_time =
+            @timed SW.WedderburnDecomposition(
             T,
             G,
             action,
@@ -160,15 +144,9 @@ function sos_problem(
             @timed sos_problem(poly, wedderburn, basis_psd)
     else
         (invariant_vs, basis_cnstr), symmetry_adaptation_time = @timed let G = G
-            basis = StarAlgebras.FixedBasis(
-                basis_constraints,
-                StarAlgebras.DiracMStructure(*),
-            )
-
-            tblG = SymbolicWedderburn.Characters.CharacterTable(
-                Rational{Int},
-                G,
-            )
+            basis = SA.FixedBasis(basis_constraints, SA.DiracMStructure(*))
+
+            tblG = SW.Characters.CharacterTable(Rational{Int}, G)
             iv = invariant_vectors(tblG, action, basis)
             iv, basis
         end

src/SymbolicWedderburn.jl

@@ -10,6 +10,7 @@ import AbstractPermutations as AP
 import AbstractPermutations: degree
 import PermutationGroups as PG
 using StarAlgebras
+import StarAlgebras as SA
 
 export symmetry_adapted_basis, WedderburnDecomposition
 export basis,
@@ -25,6 +26,7 @@ import .Characters: row_echelon_form!
 import .Characters.FiniteFields
 
 include("ext_homomorphisms.jl")
+include("ext_hom_schreier.jl")
 include("actions.jl")
 include("group_action_error.jl")
 include("action_characters.jl")