diff --git a/Project.toml b/Project.toml index 095d1ab..1dbcf6d 100644 --- a/Project.toml +++ b/Project.toml @@ -10,19 +10,23 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420" Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1" [compat] -Cyclotomics = "0.2.2" +Cyclotomics = "0.2.3" GroupsCore = "0.3.1" PermutationGroups = "0.3" Primes = "0.4, 0.5" +StarAlgebras = "0.1.5" julia = "1" [extras] DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "DynamicPolynomials", "Random", "Pkg"] +test = ["Test", "DynamicPolynomials", "JuMP", "Pkg", "Random", "SCS"] diff --git a/examples/action_polynomials.jl b/examples/action_polynomials.jl index 0992c68..85a00be 100644 --- a/examples/action_polynomials.jl +++ b/examples/action_polynomials.jl @@ -1,10 +1,44 @@ -struct PermutingVariables <: SymbolicWedderburn.ByPermutations end +using DynamicPolynomials +MP = DynamicPolynomials.MultivariatePolynomials +using GroupsCore +using PermutationGroups -function SymbolicWedderburn.action(a::PermutingVariables, g::PermutationGroups.AbstractPerm, m::Monomial) +struct VariablePermutation <: SymbolicWedderburn.ByPermutations end + +function SymbolicWedderburn.action(a::VariablePermutation, g::PermutationGroups.AbstractPerm, m::Monomial) v = variables(m) return m(v => SymbolicWedderburn.action(a, g, v)) end -function SymbolicWedderburn.action(::PermutingVariables, g::PermutationGroups.AbstractPerm, v::AbstractVector) +function SymbolicWedderburn.action(::VariablePermutation, g::PermutationGroups.AbstractPerm, v::AbstractVector) return map(i -> v[i^g], eachindex(v)) end + +abstract type OnMonomials <: SymbolicWedderburn.ByLinearTransformation end + +function SymbolicWedderburn.action( + a::Union{VariablePermutation,OnMonomials}, + el::GroupElement, + term::MP.AbstractTerm, +) + return MP.coefficient(term) * SymbolicWedderburn.action(a, el, MP.monomial(term)) +end +function SymbolicWedderburn.action( + a::Union{VariablePermutation,OnMonomials}, + el::GroupElement, + poly::MP.AbstractPolynomial, +) + return sum([SymbolicWedderburn.action(a, el, term) for term in MP.terms(poly)]) +end + +function SymbolicWedderburn.decompose( + k::MP.AbstractPolynomialLike, + hom::SymbolicWedderburn.InducedActionHomomorphism, +) + # correct only if features(hom) == monomials + + indcs = [hom[mono] for mono in MP.monomials(k)] + coeffs = MP.coefficients(k) + + return indcs, coeffs +end diff --git a/examples/dihedral.jl b/examples/dihedral.jl new file mode 100644 index 0000000..1813024 --- /dev/null +++ b/examples/dihedral.jl @@ -0,0 +1,66 @@ +import GroupsCore + +struct DihedralGroup <: GroupsCore.Group + n::Int +end + +struct DihedralElement <: GroupsCore.GroupElement + n::Int + reflection::Bool + id::Int +end + +Base.one(G::DihedralGroup) = DihedralElement(G.n, false, 0) + +Base.eltype(::DihedralGroup) = DihedralElement +function Base.iterate(G::DihedralGroup, prev::DihedralElement=DihedralElement(G.n, false, -1)) + if prev.id + 1 >= G.n + if prev.reflection + return nothing + else + next = DihedralElement(G.n, true, 0) + end + else + next = DihedralElement(G.n, prev.reflection, prev.id + 1) + end + return next, next +end +Base.IteratorSize(::Type{DihedralGroup}) = Base.HasLength() + +GroupsCore.order(::Type{T}, G::DihedralGroup) where {T} = convert(T, 2G.n) +GroupsCore.gens(G::DihedralGroup) = [DihedralElement(G.n, false, 1), DihedralElement(G.n, true, 0)] + +# Base.rand not needed for our purposes here + +Base.parent(g::DihedralElement) = DihedralGroup(g.n) +function Base.:(==)(g::DihedralElement, h::DihedralElement) + return g.n == h.n && g.reflection == h.reflection && g.id == h.id +end + +function Base.inv(el::DihedralElement) + if el.reflection || iszero(el.id) + return el + else + return DihedralElement(el.n, false, el.n - el.id) + end +end +function Base.:*(a::DihedralElement, b::DihedralElement) + a.n == b.n || error("Cannot multiply elements from different Dihedral groups") + id = mod(a.reflection ? a.id - b.id : a.id + b.id, a.n) + return DihedralElement(a.n, a.reflection != b.reflection, id) +end + +Base.copy(a::DihedralElement) = DihedralElement(a.n, a.reflection, a.id) + +# optional functions: +function GroupsCore.order(T::Type, el::DihedralElement) + if el.reflection + return T(2) + else + if iszero(el.id) + return T(1) + else + return T(div(el.n, gcd(el.n, el.id))) + end + end +end diff --git a/examples/ex_C2_linear.jl b/examples/ex_C2_linear.jl index 0f2c537..ab17aff 100644 --- a/examples/ex_C2_linear.jl +++ b/examples/ex_C2_linear.jl @@ -9,14 +9,14 @@ G = CyclicGroup(2) x,y = @polyvar x y # (x, y) → (x+y, x-y) -struct MyAction <: SymbolicWedderburn.ByLinearTransformation end +struct By90Rotation <: SymbolicWedderburn.ByLinearTransformation end # x -> linear combination of basis elements; -SymbolicWedderburn.coeff_type(m::MyAction) = Float64 -function SymbolicWedderburn.action(a::MyAction, g::CyclicGroupElement, m::Monomial) +SymbolicWedderburn.coeff_type(::By90Rotation) = Float64 +function SymbolicWedderburn.action(::By90Rotation, g::CyclicGroupElement, m::Monomial) isone(g) && return m x, y = variables(m) - return m(x => (x+y)/sqrt(2), y => (x-y)/sqrt(2)) # action must be orthogonal + return m(x => (x+y)/sqrt(2), y => (x-y)/sqrt(2)) # we need the action to be orthogonal end function SymbolicWedderburn.decompose(k::AbstractPolynomial, hom::SymbolicWedderburn.InducedActionHomomorphism) @@ -36,13 +36,13 @@ using Test g = gens(G, 1) basis = monomials([x,y], 0:4) - k = SymbolicWedderburn.action(MyAction(), g, basis[2]) - ehom = SymbolicWedderburn.ExtensionHomomorphism(MyAction(), basis) + k = SymbolicWedderburn.action(By90Rotation(), g, basis[2]) + ehom = SymbolicWedderburn.ExtensionHomomorphism(By90Rotation(), basis) idcs, vals = SymbolicWedderburn.decompose(k, ehom) @test sum(c*basis[i] for (c,i) in zip(vals, idcs)) == k for b in basis - k = SymbolicWedderburn.action(MyAction(), g, b) + k = SymbolicWedderburn.action(By90Rotation(), g, b) idcs, vals = SymbolicWedderburn.decompose(k, ehom) @test sum(c*basis[i] for (c,i) in zip(vals, idcs)) == k end @@ -52,17 +52,17 @@ end @testset "induced Matrix Representation" begin g = gens(G, 1) - basis = monomials([x,y], 0:4) - ehom = SymbolicWedderburn.ExtensionHomomorphism(MyAction(), basis) - m = droptol!(SymbolicWedderburn.induce(MyAction(), ehom, g), 1e-15) + monomial_basis = monomials([x,y], 0:4) + ehom = SymbolicWedderburn.ExtensionHomomorphism(By90Rotation(), monomial_basis) + m = droptol!(SymbolicWedderburn.induce(By90Rotation(), ehom, g), 1e-15) @test eltype(m) == Float64 @test !(m ≈ one(m)) @test m^2 ≈ one(m) - ssimple_basis = SymbolicWedderburn.symmetry_adapted_basis(G, basis, MyAction()); + ssimple_basis = SymbolicWedderburn.symmetry_adapted_basis(Float64, G, By90Rotation(), monomial_basis); degs = SymbolicWedderburn.degree.(ssimple_basis) - mlts = SymbolicWedderburn.multiplicity.(ssimple_basis) - @test sum(d*m for (d,m) in zip(degs, mlts)) == length(basis) - @test sum(first∘size∘SymbolicWedderburn.basis, ssimple_basis) == length(basis) + mlts = multiplicity.(ssimple_basis) + @test sum(d*m for (d,m) in zip(degs, mlts)) == length(monomial_basis) + @test sum(first∘size, ssimple_basis) == length(monomial_basis) end diff --git a/examples/ex_C4.jl b/examples/ex_C4.jl deleted file mode 100644 index 7a4f2dc..0000000 --- a/examples/ex_C4.jl +++ /dev/null @@ -1,69 +0,0 @@ -using SymbolicWedderburn -using PermutationGroups -using Cyclotomics - -using SparseArrays - -using DynamicPolynomials -using SumOfSquares -using SCS - -OPTIMIZER = optimizer_with_attributes( - SCS.Optimizer, - "acceleration_lookback" => 10, - "max_iters" => 10_000, - "eps" => 1e-6, - "linear_solver" => SCS.DirectSolver, -) - -include(joinpath(@__DIR__, "action_polynomials.jl")) - -N = 4 - -@polyvar x[1:N] - -f = sum(x) + - sum(x .^ 2) + - (sum((x .+ 1) .^ 2 .* (x .+ x') .^ 2))^2 * (1 + sum(x .^ 2)) - -basis = monomials(x, 0:(DynamicPolynomials.maxdegree(f)÷2)) - -ts, obj, st = let f = f, basis = basis, m = SOSModel(OPTIMIZER) - @variable m t - @objective m Max t - @variable m sos SOSPoly(basis) - @constraint m f - t == sos - optimize!(m) - @info (m,) termination_status(m) objective_value(m) solve_time(m) - termination_status(m), objective_value(m), solve_time(m) -end - -ts_sa, obj_sa, st_sa = let f = f, - basis = basis, - m = SOSModel(OPTIMIZER), - # G = PermGroup(Perm([2:N; 1])), - G = PermGroup([perm"(1,2)", Perm([2:N; 1])]) - - t = @timed let - ssimple_basis = SymbolicWedderburn.symmetry_adapted_basis(Float64, G, basis, PermutingVariables()) - SparseMatrixCSC{Float64,Int}.(SymbolicWedderburn.basis.(ssimple_basis)) - end - - sa_basis, symmetry_adaptation_time = t.value, t.time - - let m = m, basis = basis, R = sa_basis - @variable m t - @objective m Max t - - soses = @variable m [r in R] SOSPoly(FixedPolynomialBasis(r * basis)) - @constraint m f - t == sum(soses) - - optimize!(m) - @info (m,) termination_status(m) objective_value(m) solve_time(m) symmetry_adaptation_time - termination_status(m), objective_value(m), solve_time(m) - end -end - -@assert ts == ts_sa == SumOfSquares.MOI.OPTIMAL -@assert isapprox(obj, obj_sa, atol=1e-4) -@assert st_sa < st diff --git a/examples/ex_S4.jl b/examples/ex_S4.jl new file mode 100644 index 0000000..f52fa7f --- /dev/null +++ b/examples/ex_S4.jl @@ -0,0 +1,161 @@ +using SymbolicWedderburn +using PermutationGroups +using Cyclotomics + +using SparseArrays + +using DynamicPolynomials +using SumOfSquares +using SCS + +const OPTIMIZER = optimizer_with_attributes( + SCS.Optimizer, + "acceleration_lookback" => 10, + "max_iters" => 10_000, + "eps" => 1e-6, + "linear_solver" => SCS.DirectSolver, +) + +include(joinpath(@__DIR__, "action_polynomials.jl")) + +N = 4 + +@polyvar x[1:N] + +f = sum(x) + + sum(x .^ 2) + + (sum((x .+ 1) .^ 2 .* (x .+ x') .^ 2))^2 * (1 + sum(x .^ 2)) + +no_symmetry = let f = f, basis = monomials(x, 0:(DynamicPolynomials.maxdegree(f)÷2)) + stats = @timed let m = SOSModel(OPTIMIZER) + @variable m t + @objective m Max t + @variable m sos SOSPoly(basis) + @constraint m f - t == sos + m + end + m = stats.value + optimize!(m) + @info (m,) termination_status(m) objective_value(m) creation_time=(stats.time) solve_time(m) + ( + status=termination_status(m), + objective=objective_value(m), + creation_t=stats.time, + solve_t=solve_time(m) + ) +end + +semisimple_dec = let f = f, G = PermGroup([perm"(1,2)", Perm([2:N; 1])]) + basis = monomials(x, 0:(DynamicPolynomials.maxdegree(f)÷2)) + + stats = @timed begin + sa_basis, symmetry_adaptation_time, = @timed let + ssimple_basis = SymbolicWedderburn.symmetry_adapted_basis( + Rational{Int}, + G, + VariablePermutation(), + basis, + semisimple=true + ) + sp = sparse.(ssimple_basis) + droptol!.(sp, 1e-12) + end + + m = SOSModel(OPTIMIZER) + + let basis = basis, R = sa_basis + @variable m t + @objective m Max t + + soses = @variable m [r in R] SOSPoly(FixedPolynomialBasis(r * basis)) + @constraint m f - t == sum(soses) + end + m + end + m = stats.value + optimize!(m) + @info (m,) termination_status(m) objective_value(m) solve_time(m) creation_time=(stats.time) symmetry_adaptation_time + ( + status=termination_status(m), + objective=objective_value(m), + creation_t=stats.time, + solve_t=solve_time(m) + ) +end + +@assert no_symmetry.status == MOI.OPTIMAL +@assert semisimple_dec.status == MOI.OPTIMAL +@assert isapprox(no_symmetry.objective, semisimple_dec.objective, atol=1e-4) +@assert semisimple_dec.solve_t < no_symmetry.solve_t + +# This is faster, but not quite the speedup we'd like to see. +# Let us try now to use the decomposition into simple summands +include(joinpath(@__DIR__, "util.jl")) + +wedderburn_dec = let f=f, G = PermGroup([perm"(1,2)", Perm([2:N; 1])]), T = Float64 + + basis_half = monomials(x, 0:(DynamicPolynomials.maxdegree(f)÷2)) + + wedderburn, symmetry_adaptation_time, = @timed WedderburnDecomposition( + T, + G, + VariablePermutation(), + monomials(x, 0:DynamicPolynomials.maxdegree(f)), + basis_half + ) + + M = [SymbolicWedderburn.basis(wedderburn)[x*y] for x in basis_half, y in basis_half] + + stats = @timed begin + m = let m = JuMP.Model(OPTIMIZER) + JuMP.@variable m t + JuMP.@objective m Max t + psds = [ + JuMP.@variable(m, [1:d, 1:d] in PSDCone()) + for d in size.(direct_summands(wedderburn), 1) + ] + + # preallocating + Mπs = zeros.(T, size.(psds)) + M_orb = similar(M, T) + tmps = zeros.(T, reverse.(size.(direct_summands(wedderburn)))) + + C = DynamicPolynomials.coefficients(f-t, SymbolicWedderburn.basis(wedderburn)) + + for iv in invariant_vectors(wedderburn) + c = dot(C, iv) + # C needs to be invariant under G, + # which is eqivalent to being constant on the orbits + let cc = C[findfirst(!iszero, iv)] + @assert all(C[m] == cc for m in findall(!iszero, iv)) + end + + # invariant_constraint! is defined locally in util.jl + # it essentially amounts to averaging constraint matrices defined by M + # with weights provided by iv and storing the result in (dense) matrix M_orb + M_orb = invariant_constraint!(M_orb, M, iv) + Mπs = SymbolicWedderburn.diagonalize!(Mπs, M_orb, wedderburn, tmps) + + JuMP.@constraint m sum( + dot(Mπ, Pπ) for (Mπ, Pπ) in zip(Mπs, psds) if !iszero(Mπ) + ) == c + end + m + end + end + + m = stats.value + optimize!(m) + @info (m,) termination_status(m) objective_value(m) solve_time(m) creation_time=(stats.time) symmetry_adaptation_time + ( + status=termination_status(m), + objective=objective_value(m), + creation_t=stats.time, + solve_t=solve_time(m) + ) +end + +@assert wedderburn_dec.status == MOI.OPTIMAL +@assert isapprox(wedderburn_dec.objective, semisimple_dec.objective, atol=1e-4) +@assert semisimple_dec.solve_t > wedderburn_dec.solve_t +@assert wedderburn_dec.solve_t < 2.0 diff --git a/examples/ex_motzkin.jl b/examples/ex_motzkin.jl index 302fe1a..cfd007d 100644 --- a/examples/ex_motzkin.jl +++ b/examples/ex_motzkin.jl @@ -12,7 +12,7 @@ include(joinpath(@__DIR__, "action_polynomials.jl")) OPTIMIZER = optimizer_with_attributes( SCS.Optimizer, "acceleration_lookback" => 0, - "max_iters" => 20_000, + "max_iters" => 10_000, "eps" => 2e-6, "linear_solver" => SCS.DirectSolver, ) @@ -21,9 +21,9 @@ OPTIMIZER = optimizer_with_attributes( motzkin = x^4 * y^2 + y^4 * x^2 - 3 * x^2 * y^2 + 1 g = (x^2 + y^2 + 1) -basis = monomials([x, y], 0:7) +monomial_basis = monomials([x, y], 0:7) -ts, st = let f = motzkin, basis = basis, m = SOSModel(OPTIMIZER) +ts, st = let f = motzkin, basis = monomial_basis, m = SOSModel(OPTIMIZER) @variable m t >=0 # @objective m Max t @variable m sos SOSPoly(basis) @@ -34,18 +34,22 @@ ts, st = let f = motzkin, basis = basis, m = SOSModel(OPTIMIZER) end ts_sa, st_sa = let f = motzkin, - basis = basis, - m = SOSModel(OPTIMIZER), + basis = monomial_basis, G = PermGroup(perm"(1,2)") - t = @timed let - ssimple_basis = SymbolicWedderburn.symmetry_adapted_basis(Float64, G, exponents.(basis), PermutingVariables()) - SparseMatrixCSC{Float64,Int}.(SymbolicWedderburn.basis.(ssimple_basis)) + sa_basis, symmetry_adaptation_time, = @timed let + ssimple_basis = SymbolicWedderburn.symmetry_adapted_basis( + Float64, + G, + VariablePermutation(), + exponents.(basis), + semisimple=true + ) + sp = sparse.(ssimple_basis) + droptol!.(sp, 1e-12) end - sa_basis, symmetry_adaptation_time = t.value, t.time - - let m = m, basis = basis, R = sa_basis + let m = SOSModel(OPTIMIZER), basis = basis, R = sa_basis @variable m t >= 0 # @objective m Max t diff --git a/examples/ex_robinson_form.jl b/examples/ex_robinson_form.jl new file mode 100644 index 0000000..8724a7b --- /dev/null +++ b/examples/ex_robinson_form.jl @@ -0,0 +1,269 @@ +# # Dihedral symmetry of the Robinson form + +# Based on version from SumOfSquares by Benoît Legat +# **Adapted from**: Example 5.4 of [GP04] +# +# [GP04] Gatermann, Karin and Parrilo, Pablo A. +# *Symmetry groups, semidefinite programs, and sums of squares*. +# Journal of Pure and Applied Algebra 192.1-3 (2004): 95-128. + +using LinearAlgebra +using SparseArrays + +using SymbolicWedderburn +using SymbolicWedderburn.StarAlgebras +using DynamicPolynomials +using SumOfSquares +using SCS + +const OPTIMIZER = optimizer_with_attributes( + SCS.Optimizer, + "acceleration_lookback" => 10, + "max_iters" => 3_000, + "alpha" => 1.95, + "eps" => 1e-6, + "linear_solver" => SCS.DirectSolver, +) + +@polyvar x y +const robinson_form = x^6 + y^6 - x^4 * y^2 - y^4 * x^2 - x^4 - y^4 - x^2 - y^2 + 3x^2 * y^2 + 1 + +# The Robinson form is invariant under the following action of the Dihedral group on monomials: +# The action of each element of the groups is to map the variables `x, y` to: +# +# | id | rotation | reflection | +# |----|----------|------------| +# | 0 | x, y | y, x | +# | 1 | -y, x | -x, y | +# | 2 | -x, -y | -y, -x | +# | 3 | y, -x | x, -y | + +# implementation of the dihedral group: +include(joinpath(@__DIR__, "dihedral.jl")) +# definitions of general actions on polynomials: +include(joinpath(@__DIR__, "action_polynomials.jl")) + +struct DihedralAction <: OnMonomials end + +SymbolicWedderburn.coeff_type(::DihedralAction) = Float64 +function SymbolicWedderburn.action(::DihedralAction, el::DihedralElement, mono::AbstractMonomial) + if iseven(el.reflection + el.id) + var_x, var_y = x, y + else + var_x, var_y = y, x + end + sign_x = 1 <= el.id <= 2 ? -1 : 1 + sign_y = 2 <= el.id ? -1 : 1 + return mono([x, y] => [sign_x * var_x, sign_y * var_y]) +end + +G = DihedralGroup(4) +for g in G + @assert SymbolicWedderburn.action(DihedralAction(), g, robinson_form) == robinson_form +end + +# A simple Gramm-matrix like formulation of SOS problem: +# * a single large pdd constraint and +# * one linear constraint for each monomial in the monomial basis of f. +no_symmetry = let f = robinson_form + stats = @timed let m = Model(OPTIMIZER) + @variable m t + @objective m Max t + @constraint m f - t in SOSCone() + m + end + m = stats.value + optimize!(m) + @info (m,) termination_status(m) objective_value(m) creation_time = (stats.time) solve_time(m) + + @assert isapprox(value(m[:t]), -3825 / 4096, rtol = 1e-2) + + ( + status = termination_status(m), + objective = objective_value(m), + creation_t = stats.time, + solve_t = solve_time(m) + ) +end + +# Let's decompose basis into semi-simple summands and prepare the approximately +# reduced optimization problem. Since this is an isomorphism we reduce the +# number of variables from N² to n₁² + … + nₖ², where Σᵢ nᵢ = N. That is instead +# of one psd constraint of size N we have k constraints of sizes n₁,… nₖ. +# The number of linear constraints remains the same as before. +semisimple_dec = let f = robinson_form, G = DihedralGroup(4) + basis = monomials([x, y], 0:(DynamicPolynomials.maxdegree(f) ÷ 2)) + + stats = @timed begin + sa_basis, symmetry_adaptation_time, = @timed let + ssimple_basis = SymbolicWedderburn.symmetry_adapted_basis( + Rational{Int}, + G, + DihedralAction(), + basis, + semisimple = true + ) + sp = sparse.(ssimple_basis) + droptol!.(sp, 1e-12) + end + + m = SOSModel(OPTIMIZER) + + let basis = basis, R = sa_basis + @variable m t + @objective m Max t + + soses = @variable m [r in R] SOSPoly(FixedPolynomialBasis(r * basis)) + @constraint m f - t == sum(soses) + end + m + end + m = stats.value + optimize!(m) + + @assert isapprox(value(m[:t]), -3825 / 4096, rtol = 1e-2) + + @info (m,) termination_status(m) objective_value(m) solve_time(m) creation_time = (stats.time) symmetry_adaptation_time + ( + status = termination_status(m), + objective = objective_value(m), + creation_t = stats.time, + solve_t = solve_time(m) + ) +end + +# Now given that f must be constant on orbits of monomials under the action of G +# let's exploit that to reduce the number of cnstraints (one for each orbit). +# Here we keep the single large psd constraint of size N×N. +include(joinpath(@__DIR__, "util.jl")) + +orbit_dec = let f = robinson_form, G = DihedralGroup(4), T = Float64 + basis_monoms = monomials([x,y], 0:DynamicPolynomials.maxdegree(f)) + + invariant_vs, symmetry_adaptation_time, = @timed invariant_vectors(G, DihedralAction(), basis_monoms) + + M = let basis_full = StarAlgebras.Basis{UInt32}(basis_monoms), basis_half = monomials([x,y], 0:(DynamicPolynomials.maxdegree(f) ÷ 2)) + [basis_full[x * y] for x in basis_half, y in basis_half] + end + + stats = @timed let m = JuMP.Model(OPTIMIZER) + JuMP.@variable m t + JuMP.@objective m Max t + P = JuMP.@variable(m, [1:size(M, 1), 1:size(M, 1)], PSD) + + # preallocating + M_orb = similar(M, T) + + C = DynamicPolynomials.coefficients(f - t, basis_monoms) + + for iv in invariant_vs + c = dot(C, iv) + + # invariant_constraint! is defined locally in util.jl + # it essentially amounts to averaging constraint matrices defined by M + # with weights provided by iv and storing the result in (dense) matrix M_orb + M_orb = invariant_constraint!(M_orb, M, iv) + + #= + if !iszero(M_orb) + monoms = basis_monoms[unique!(M[map(!iszero, M_orb)])] + @info "Coefficient of the $monoms: $([DynamicPolynomials.coefficient(f-t, m) for m in monoms])" dot_C_v=c + @show sparse(M_orb) + else + monoms = basis_monoms[unique!(M[SparseArrays.nonzeroinds(iv)])] + @info "Coefficient of the $monoms: $([DynamicPolynomials.coefficient(f-t, m) for m in monoms])" dot_C_v=c + @info "M_orb is zero!" + end + =# + + JuMP.@constraint m dot(M_orb, P) == c + end + m + end + + m = stats.value + optimize!(m) + + @assert isapprox(value(m[:t]), -3825 / 4096, rtol = 1e-2) + + @info (m,) termination_status(m) objective_value(m) solve_time(m) creation_time = (stats.time) symmetry_adaptation_time + ( + status = termination_status(m), + objective = objective_value(m), + creation_t = stats.time, + solve_t = solve_time(m) + ) +end + + +wedderburn_dec = let f = robinson_form, G = DihedralGroup(4), T = Float64 + + basis_half = monomials([x,y], 0:(DynamicPolynomials.maxdegree(f) ÷ 2)) + + wedderburn, symmetry_adaptation_time, = @timed WedderburnDecomposition( + T, + G, + DihedralAction(), + monomials([x,y], 0:DynamicPolynomials.maxdegree(f)), + basis_half + ) + + M = [SymbolicWedderburn.basis(wedderburn)[x * y] for x in basis_half, y in basis_half] + + stats = @timed let m = JuMP.Model(OPTIMIZER) + JuMP.@variable m t + JuMP.@objective m Max t + psds = [ + JuMP.@variable(m, [1:d, 1:d] in PSDCone()) + for d in size.(direct_summands(wedderburn), 1) + ] + + # preallocating + Mπs = zeros.(T, size.(psds)) + M_orb = similar(M, T) + tmps = zeros.(T, reverse.(size.(direct_summands(wedderburn)))) + + C = DynamicPolynomials.coefficients(f - t, SymbolicWedderburn.basis(wedderburn)) + + for iv in invariant_vectors(wedderburn) + c = dot(C, iv) + + # invariant_constraint! is defined locally in util.jl + # it essentially amounts to averaging constraint matrices defined by M + # with weights provided by iv and storing the result in (dense) matrix M_orb + M_orb = invariant_constraint!(M_orb, M, iv) + + #= + if !iszero(M_orb) + monoms = SymbolicWedderburn.basis(wedderburn)[unique!(M[map(!iszero, M_orb)])] + @info "Coefficient of the $monoms: $([DynamicPolynomials.coefficient(f-t, m) for m in monoms])" dot_C_v=c + @show sparse(M_orb) + else + monoms = SymbolicWedderburn.basis(wedderburn)[unique!(M[SparseArrays.nonzeroinds(iv)])] + @info "Coefficient of the $monoms: $([DynamicPolynomials.coefficient(f-t, m) for m in monoms])" dot_C_v=c + @info "M_orb is zero!" + end + =# + + Mπs = SymbolicWedderburn.diagonalize!(Mπs, M_orb, wedderburn, tmps) + + JuMP.@constraint m sum( + dot(Mπ, Pπ) for (Mπ, Pπ) in zip(Mπs, psds) if !iszero(Mπ) + ) == c + end + m + end + m = stats.value + optimize!(m) + + @info value(m[:t]) + @assert isapprox(value(m[:t]), -3825 / 4096, rtol = 1e-2) + + @info (m,) termination_status(m) objective_value(m) solve_time(m) creation_time = (stats.time) + ( + status = termination_status(m), + objective = objective_value(m), + creation_t = stats.time, + solve_t = solve_time(m) + ) +end diff --git a/examples/util.jl b/examples/util.jl new file mode 100644 index 0000000..ee98280 --- /dev/null +++ b/examples/util.jl @@ -0,0 +1,16 @@ +using SparseArrays + +function invariant_constraint!( + M_orb::AbstractMatrix{<:AbstractFloat}, + M::Matrix{<:Integer}, + invariant_vec::SparseVector, +) + M_orb .= zero(eltype(M_orb)) + for i in eachindex(M) + if M[i] ∈ SparseArrays.nonzeroinds(invariant_vec) + M_orb[i] += invariant_vec[M[i]] + end + end + return M_orb +end + diff --git a/src/Characters/Characters.jl b/src/Characters/Characters.jl new file mode 100644 index 0000000..9f0ba65 --- /dev/null +++ b/src/Characters/Characters.jl @@ -0,0 +1,30 @@ +module Characters +using LinearAlgebra +using SparseArrays +using Primes +using GroupsCore +using Cyclotomics + +import PermutationGroups: AbstractOrbit, degree +import PermutationGroups + +export AbstractClassFunction, Character, CharacterTable + +export conjugacy_classes, + constituents, + degree, + irreducible_characters, + isirreducible, + table + +include("gf.jl") + +include("echelon_form.jl") +include("eigenspacedecomposition.jl") +include("cmmatrix.jl") + +include("powermap.jl") +include("character_tables.jl") +include("class_functions.jl") +include("dixon.jl") +end diff --git a/src/Characters/character_tables.jl b/src/Characters/character_tables.jl new file mode 100644 index 0000000..79d0927 --- /dev/null +++ b/src/Characters/character_tables.jl @@ -0,0 +1,200 @@ +struct CharacterTable{Gr,T,O} <: AbstractMatrix{T} + group::Gr + conjugacy_classes::Vector{O} + inv_of::Vector{Int} + pmap::PowerMap + values::Matrix{T} + # rows indexed by irreducible characters, + # columns indexed by conjugacy classes +end + +## Array Interface: +Base.size(chtbl::CharacterTable) = size(chtbl.values) +Base.getindex(chtbl::CharacterTable, i::Integer, j::Integer) = chtbl.values[i, j] + +## Accessors +powermap(chtbl::CharacterTable) = chtbl.pmap +Base.parent(chtbl::CharacterTable) = chtbl.group + +conjugacy_classes(chtbl::CharacterTable) = chtbl.conjugacy_classes +nconjugacy_classes(chtbl::CharacterTable) = size(chtbl.values, 2) +nirreps(chtbl::CharacterTable) = size(chtbl.values, 1) + +## irreps +irreducible_characters(chtbl::CharacterTable) = + [Character(chtbl, i) for i in 1:size(chtbl, 1)] +irreducible_characters(T::Type, chtbl::CharacterTable) = + [Character{T}(chtbl, i) for i in 1:size(chtbl, 1)] +irreducible_characters(G::Group, cclasses = conjugacy_classes(G)) = + irreducible_characters(Rational{Int}, G, cclasses) + +function irreducible_characters( + R::Type{<:Rational}, + G::Group, + cclasses=conjugacy_classes(G) +) + return irreducible_characters(CharacterTable(R, G, cclasses)) +end + +## construcing tables + +function CharacterTable( + Fp::Type{<:FiniteFields.GF}, + G::Group, + cclasses = conjugacy_classes(G), +) + Ns = [CMMatrix(cclasses, i) for i in 1:length(cclasses)] + esd = common_esd(Ns, Fp) + @assert isdiag(esd) + + tbl = CharacterTable(G, cclasses, _inv_of(cclasses), PowerMap(cclasses), esd.basis) + + tbl = normalize!(tbl) + return tbl +end + +function CharacterTable(R::Type{<:Rational}, G::Group, cclasses=conjugacy_classes(G)) + Fp = FiniteFields.GF{dixon_prime(cclasses)} + tblFp = CharacterTable(Fp, G, cclasses) + return complex_character_table(R, tblFp) +end + +function complex_character_table( + R::Type{<:Rational}, + tblFp::CharacterTable{<:Group, <:FiniteFields.GF}, +) + charsFp = irreducible_characters(tblFp) + mult_c = _multiplicities(charsFp) + + e = size(mult_c, 3) # the exponent + + C = Cyclotomics.Cyclotomic{R,Cyclotomics.SparseVector{R,Int}} + values = Matrix{C}(undef, size(tblFp)) + + Es = [E(e, k) for k in 0:e - 1] + Threads.@threads for j in 1:size(tblFp, 2) # conjugacy_classes + for i in 1:size(tblFp, 1) # characters + # reduced_embedding may prevent overflow sometimes + values[i, j] = Cyclotomics.reduced_embedding( + sum(mult_c[i, j, k + 1] * Es[k + 1] for k in 0:e - 1) + ) + end + end + + return CharacterTable( + parent(tblFp), + conjugacy_classes(tblFp), + tblFp.inv_of, + powermap(tblFp), + values, + ) +end + +function _inv_of(cc::AbstractVector{<:AbstractOrbit}) + inv_of = zeros(Int, size(cc)) + for (i, c) in enumerate(cc) + g = inv(first(c)) + inv_of[i] = something(findfirst(k -> g in k, cc), 0) + end + any(iszero, inv_of) && throw( + ArgumentError( + "Could not find the conjugacy class for inverse of $(first(cc[findfirst(iszero, inv_of)])).", + ), + ) + return inv_of +end + +function normalize!(chtbl::CharacterTable{<:Group,<:FiniteFields.GF}) + id = one(parent(chtbl)) + for (i, χ) in enumerate(irreducible_characters(chtbl)) + k = χ(id) + if !isone(k) + chtbl.values[i, :] .*= inv(k) + end + # ⟨χ, χ⟩ = 1/d² + + deg = sqrt(inv(dot(χ, χ))) + @debug "normalizing with" dot(χ, χ) χ(id) χ + + # normalizing χ + chtbl.values[i, :] .*= deg + end + return chtbl +end + +## "fancy" io + +function Base.show(io::IO, ::MIME"text/plain", chtbl::CharacterTable) + println(io, "Character table of ", parent(chtbl), " over $(eltype(chtbl))") + + if !(get(io, :limit, false)::Bool) + screenheight = screenwidth = typemax(Int) + else + sz = displaysize(io)::Tuple{Int,Int} + screenheight, screenwidth = sz[1] - 4, sz[2] + end + + nirr = nirreps(chtbl) + nccl = nconjugacy_classes(chtbl) + pre_width = 1 + length(string(nirr)) + sep = " │ " + + if nirr > screenheight - 1 + rows_to_print = [1:screenheight - 2; nirr] + else + rows_to_print = [1:nirr;] + end + + maxpossiblecols = div(screenwidth - pre_width - 3, 3) + if nccl > maxpossiblecols + cols_to_print = [1:maxpossiblecols - 1; nccl] + else + cols_to_print = [1:nccl;] + end + + A = Base.alignment( + io, + chtbl.values, + rows_to_print, + cols_to_print, + screenwidth, + screenwidth, + 2, + ) + + hellipsis = nccl > length(A) ? " …" : "" + + print(io, " "^(pre_width), " ") + Base.print_matrix_row( + io, + reshape(cols_to_print, (1, length(cols_to_print))), + A, + 1, + cols_to_print, + " ", + ) + + println(io, hellipsis) + println(io, "─"^pre_width, "─┬─", "─"^(sum(sum, A) + length(A) - 1), hellipsis) + + if nirr > screenheight - 1 + for i in 1:screenheight - 2 + print(io, rpad("χ$(FiniteFields.subscriptify(i))", pre_width), sep) + Base.print_matrix_row(io, chtbl.values, A, i, cols_to_print, " ") + println(io, hellipsis) + end + print(io, " ⋮", " "^(pre_width - 2), sep) + Base.print_matrix_vdots(io, "⋮", A, " ", 2, 1, false) + println(io) + + print(io, rpad("χ$(FiniteFields.subscriptify(nirr))", pre_width), sep) + Base.print_matrix_row(io, chtbl.values, A, nirr, cols_to_print, " ") + print(io, hellipsis) + else + for i in 1:nirr + print(io, rpad("χ$(FiniteFields.subscriptify(i))", pre_width), sep) + Base.print_matrix_row(io, chtbl.values, A, i, cols_to_print, " ") + i != nirr && println(io, hellipsis) + end + end +end diff --git a/src/Characters/class_functions.jl b/src/Characters/class_functions.jl new file mode 100644 index 0000000..dab00ee --- /dev/null +++ b/src/Characters/class_functions.jl @@ -0,0 +1,247 @@ +""" + AbstractClassFunction{T} +Abstract type representing functions constant on conjugacy classes of a group. + +The following functionality is required for an `AbstractClassFunction`: + * `parent(χ)`: the underlying group + * `conjugacy_classes(χ)`: iterator over conjugacy classes of `χ`. + * `values(χ)`: a _typed_ iterator over values + * `getindex(χ, i::Integer)`: the value on `i`-th conjugacy class. + Note: Indexing with negative integers should return values on the class which + contains inverses of the `i`-th class. + * `χ(g::GroupElement)`: value of `χ` on a group element (only for julia < 1.3) +""" +abstract type AbstractClassFunction{T} end # <:AbstractVector{T}? + +Base.eltype(::AbstractClassFunction{T}) where T = T + +function LinearAlgebra.dot(χ::AbstractClassFunction, ψ::AbstractClassFunction) + @assert conjugacy_classes(χ) === conjugacy_classes(ψ) + val = sum( + length(cc) * χ[i] * ψ[-i] for (i, cc) in enumerate(conjugacy_classes(χ)) + ) + orderG = sum(length, conjugacy_classes(χ)) + val = _div(val, orderG) + return val +end + +_div(val, orderG) = div(val, orderG) +_div(val::AbstractFloat, orderG) = val / orderG +_div(val::Complex{<:AbstractFloat}, orderG) = val / orderG + +## Arbitrary ClassFunctions without decomposition into irreps +struct ClassFunction{T,CCl} <: AbstractClassFunction{T} + vals::Vector{T} + conjugacy_classes::CCl + inv_of::Vector{Int} +end + +ClassFunction(vals, cclasses) = ClassFunction(vals, cclasses, _inv_of(cclasses)) + +## AbstractClassFunction api +Base.parent(χ::ClassFunction) = parent(first(first(conjugacy_classes(χ)))) +conjugacy_classes(χ::ClassFunction) = χ.conjugacy_classes +Base.values(χ::ClassFunction) = χ.vals + +Base.@propagate_inbounds function Base.getindex(χ::ClassFunction, i::Integer) + i = i < 0 ? χ.inv_of[-i] : i + @boundscheck 1 ≤ i ≤ length(conjugacy_classes(χ)) + return values(χ)[i] +end + +# TODO: move to AbstractClassFunction when we drop julia-1.0 +# adding call methods to abstract types was not supported before julia-1.3 +function (χ::ClassFunction)(g::GroupElement) + for (i, cc) in enumerate(conjugacy_classes(χ)) + g ∈ cc && return χ[i] + end + throw( + DomainError( + g, + "element does not belong to conjugacy classes of $χ", + ), + ) +end + +""" + Character <: AbstractClassFunction +Struct representing (possibly virtual) character of a group. + +Characters are backed by `table(χ)::CharacterTable` which actually stores the +character values. The constituents (decomposition into the irreducible summands) +of a given character can be obtained by calling `constituents(χ)` which returns a +vector of coefficients of `χ` in the basis of `irreducible_characters(table(χ))`. + +It is assumed that equal class functions on the same group will have **identical** +(ie. `===`) character tables. +""" +struct Character{T, S, ChT<:CharacterTable} <: AbstractClassFunction{T} + table::ChT + constituents::Vector{S} +end + +function Character( + chtbl::CharacterTable{Gr, T}, + constituents::AbstractVector{S} +) where {Gr, T, S} + R = Base._return_type(*, Tuple{S,T}) + return Character{R, S, typeof(chtbl)}(chtbl, constituents) +end + +Character(chtbl::CharacterTable, i::Integer) = Character{eltype(chtbl)}(chtbl, i) + +function Character{T}(chtbl::CharacterTable, i::Integer) where T + v = zeros(Int, nconjugacy_classes(chtbl)) + v[i] = 1 + return Character{T, Int, typeof(chtbl)}(chtbl, v) +end + +function Character{T}(χ::Character) where T + S = eltype(constituents(χ)) + ChT = typeof(table(χ)) + return Character{T, S, ChT}(table(χ), constituents(χ)) +end + +## Accessors +table(χ::Character) = χ.table +constituents(χ::Character) = χ.constituents + +## AbstractClassFunction api +Base.parent(χ::Character) = parent(table(χ)) +conjugacy_classes(χ::Character) = conjugacy_classes(table(χ)) +Base.values(χ::Character) = (χ[i] for i in 1:nconjugacy_classes(table(χ))) + +Base.@propagate_inbounds function Base.getindex(χ::Character{T}, i::Integer) where T + i = i < 0 ? table(χ).inv_of[abs(i)] : i + @boundscheck 1 ≤ i ≤ nconjugacy_classes(table(χ)) + return T( + sum(c * table(χ)[idx, i] for (idx, c) in enumerate(constituents(χ)) + if !iszero(c)) + ) +end + +# TODO: move to AbstractClassFunction when we drop julia-1.0 +function (χ::Character)(g::GroupElement) + for (i, cc) in enumerate(conjugacy_classes(χ)) + g ∈ cc && return χ[i] + end + throw( + DomainError( + g, + "element does not belong to conjugacy classes of $χ", + ), + ) +end + +## Basic functionality + +Base.:(==)(χ::Character, ψ::Character) = + table(χ) === table(ψ) && constituents(χ) == constituents(ψ) +Base.hash(χ::Character, h::UInt) = hash(table(χ), hash(constituents(χ), h)) + +Base.deepcopy_internal(χ::Character, ::IdDict) = + Character(table(χ), copy(constituents(χ))) + +## Character arithmetic + +for f in (:+, :-) + @eval begin + function Base.$f(χ::Character, ψ::Character) + @assert table(χ) === table(ψ) + return Character(table(χ), $f(constituents(χ), constituents(ψ))) + end + end +end + +Base.:*(χ::Character, c::Number) = Character(table(χ), c .* constituents(χ)) +Base.:*(c::Number, χ::Character) = χ * c +Base.:/(χ::Character, c::Number) = Character(table(χ), constituents(χ) ./ c) + +## Group-theoretic functions: + +PermutationGroups.degree(χ::Character) = Int(χ(one(parent(χ)))) +PermutationGroups.degree( + χ::Character{T,CCl}, +) where {T,CCl <: AbstractOrbit{<:AbstractMatrix}} = Int(χ[1]) + +function Base.conj(χ::Character{T, S}) where {T, S} + vals = collect(values(χ)) + all(isreal, vals) && return Character{T}(χ) + tbl = table(χ) + ψ = ClassFunction(vals[tbl.inv_of], conjugacy_classes(tbl), tbl.inv_of) + constituents = S[dot(ψ, χ) for χ in irreducible_characters(tbl)] + return Character{T, eltype(constituents), typeof(tbl)}(tbl, constituents) +end + +isvirtual(χ::Character) = + any(<(0), constituents(χ)) || any(!isinteger, constituents(χ)) + +function isirreducible(χ::Character) + count(isone, constituents(χ)) == 1 || return false + count(iszero, constituents(χ)) == length(constituents(χ)) - 1 || return false + return true +end + +""" + affordable_real!(χ::Character) +Return either `χ` or `2re(χ)` depending whether `χ` is afforded by a real +representation, modifying `χ` in place. +""" +function affordable_real!(χ::Character) + ι = frobenius_schur(χ) + if ι <= 0 # i.e. χ is complex or quaternionic + χ.constituents .+= constituents(conj(χ)) + end + return χ +end + +""" + frobenius_schur(χ::AbstractClassFunction[, pmap::PowerMap]) +Return Frobenius-Schur indicator of `χ`, i.e. `Σχ(g²)` where sum is taken over +the whole group. + +If χ is an irreducible `Character`, Frobenius-Schur indicator takes values in +`{1, 0, -1}` which correspond to the following situations: + 1. `χ` is real-valued and is afforded by an irreducible real representation, + 2. `χ` is a complex character which is not afforded by a real representation, and + 3. `χ` is quaternionic character, i.e. it is real valued, but is not afforded by a + real representation. + +In cases 2. and 3. `2re(χ) = χ + conj(χ)` corresponds to an irreducible character +afforded by a real representation. +""" +function frobenius_schur(χ::Character) + @assert isirreducible(χ) + + pmap = powermap(table(χ)) + ι = sum( + length(c) * χ[pmap[i, 2]] for (i, c) in enumerate(conjugacy_classes(χ)) + ) + + ι_int = Int(ι) + ordG = sum(length, conjugacy_classes(χ)) + d, r = divrem(ι_int, ordG) + @assert r == 0 "Non integral Frobenius Schur Indicator: $(ι_int) = $d * $ordG + $r" + return d +end + +Base.isreal(χ::Character) = frobenius_schur(χ) > 0 + +function Base.show(io::IO, ::MIME"text/plain", χ::Character) + println(io, "Character over ", eltype(χ)) + _print_char(io, χ) +end + +Base.show(io::IO, χ::Character) = _print_char(io, χ) + +function _print_char(io::IO, χ::Character) + first = true + for (i, c) in enumerate(constituents(χ)) + iszero(c) && continue + first || print(io, " ") + print(io, ((c < 0 || first) ? "" : '+')) + !isone(c) && print(io, c, '·') + print(io, 'χ', FiniteFields.subscriptify(i)) + first = false + end +end diff --git a/src/cmmatrix.jl b/src/Characters/cmmatrix.jl similarity index 92% rename from src/cmmatrix.jl rename to src/Characters/cmmatrix.jl index 679ecfa..3ad476c 100644 --- a/src/cmmatrix.jl +++ b/src/Characters/cmmatrix.jl @@ -45,12 +45,12 @@ function conjugacy_classes_orbit(G::GroupsCore.Group) S = gens(G) ordG = order(Int, G) - cclasses = [Orbit([id], Dict(id => nothing))] + cclasses = [PermutationGroups.Orbit([id], Dict(id => nothing))] elts_counted = 1 for g in G any(ccl -> g ∈ ccl, cclasses) && continue - ccl_g = Orbit(S, g, ^) + ccl_g = PermutationGroups.Orbit(S, g, ^) elts_counted += length(ccl_g) push!(cclasses, ccl_g) elts_counted == ordG && break diff --git a/src/Characters/dixon.jl b/src/Characters/dixon.jl new file mode 100644 index 0000000..a788e4d --- /dev/null +++ b/src/Characters/dixon.jl @@ -0,0 +1,69 @@ +Base.exponent(G::GroupsCore.Group) = exponent(conjugacy_classes(G)) +Base.exponent(cclasses::AbstractArray) = lcm(order.(first.(cclasses))) +dixon_prime(G::GroupsCore.Group) = dixon_prime(order(G), exponent(G)) + +function dixon_prime(cclasses::AbstractArray) + ordG = sum(length, cclasses) + m = exponent(cclasses) + return dixon_prime(ordG, m) +end + +function dixon_prime(ordG::Integer, exponent::Integer) + # `FiniteFields.GF{p}` needs `p` to be `Int` so no need to do the + # computations of this function over `BigInt` if `ordG` is so we convert + # to `Int`. + p = 2 * convert(Int, isqrt(ordG)) + while true + p = nextprime(p + 1) + isone(p % exponent) && break # we need -1 to be in the field + end + return p +end + +function common_esd(Ns, F::Type{FiniteFields.GF{q}}) where q + @assert isprime(q) + itr = iterate(Ns) + @assert itr !== nothing + esd = EigenSpaceDecomposition(F.(first(itr))) + for N in Iterators.rest(Ns, last(itr)) + esd = refine(esd, F.(N, false)) + @debug N esd.eigspace_ptrs + isdiag(esd) && return esd + end + return esd +end + +function _multiplicities( + chars::AbstractVector{<:Character{F}}, + cclasses = conjugacy_classes(first(chars)), +) where {F<:FiniteFields.GF} + + e = Int(exponent(cclasses)) + ie = inv(F(e)) + + ω⁻¹ = inv(FiniteFields.rootofunity(F, e)) + ωs = Matrix{typeof(ω⁻¹)}(undef, e, e) + Threads.@threads for k in 0:e-1 + ω⁻ᵏ = ω⁻¹^k + for l in 0:e-1 + ωs[l+1, k+1] = ω⁻ᵏ^l + end + end + + multiplicities = zeros(Int, length(chars), length(cclasses), e) + pmap = powermap(table(first(chars))) + # pmap is a lazy object, computed on demand; to ensure its thread-safety + # we precompute all its values here so in the loops below it is accessed read-only + collect(pmap) + + for (i, χ) in enumerate(chars) + Threads.@threads for j = 1:length(cclasses) + for k = 0:e-1 + val = Int(ie * sum(χ[pmap[j, l]] * ωs[l+1, k+1] for l = 0:e-1)) + multiplicities[i, j, k+1] = val + end + end + end + + return multiplicities +end diff --git a/src/Characters/echelon_form.jl b/src/Characters/echelon_form.jl new file mode 100644 index 0000000..581a3e9 --- /dev/null +++ b/src/Characters/echelon_form.jl @@ -0,0 +1,97 @@ +# Version over exact field: +function row_echelon_form!(A::AbstractMatrix{T}) where {T<:Union{FiniteFields.GF, Cyclotomic{<:Rational}, Rational}} + pivots = Int[] + pos = 0 + for i = 1:size(A, 2) + j = findnext(x -> !iszero(x), @view(A[:, i]), pos+1) + j === nothing && continue + pos += 1 + push!(pivots, i) + + # swap the rows so that + if (pos != j) + for k in 1:size(A, 2) + A[pos, k], A[j, k] = A[j, k], A[pos, k] + end + end + # A[pos, :] is the row with leading nz in i-th column + + w = inv(A[pos, i]) + # all columns to the left of i are zeroed (below pivots) so we start at i + for colidx in i:size(A, 2) + A[pos, colidx] *= w + end + # A[pos, i] is 1 now + + # zero the whole i-th column above and below pivot: + for rowidx = 1:size(A, 1) + rowidx == pos && continue # don't zero the active row + v = A[rowidx, i] + iszero(v) && continue + for k in i:size(A, 2) # to the left of pivot everything is zero + A[rowidx, k] -= v * A[pos, k] + end + end + end + return A, pivots +end + +function _findmax(f, A::AbstractArray) + mval, midx = f(first(A)), firstindex(A) + for idx in eachindex(A) + if (v = f(A[idx])) > mval + mval, midx = v, idx + end + end + return mval, midx +end + +function row_echelon_form!(A::AbstractMatrix{T}, atol=eps(T)*max(size(A)...)) where {T<:AbstractFloat} + pivots = Int[] + pos = 0 + sgn = 1.0 + for i = 1:size(A, 2) + # find the largest entry in the column below the last pivot + mval, midx = _findmax(abs, @view(A[pos+1:end, i])) + # j = findnext(x -> !iszero(x), , pos+1) + if mval < atol # the largest entry is below threshold so we zero everything in the column + A[pos+1:end, i] .= zero(T) + continue + end + j = midx+pos + pos += 1 + push!(pivots, i) + + # swap the rows so that + if (pos != j) + for k in 1:size(A, 2) + A[pos, k], A[j, k] = A[j, k], A[pos, k] + end + end + # A[pos, :] is the row with leading nz in i-th column + + w = inv(A[pos, i]) + sgn *= sign(w) + # all columns to the left of i are zeroed (below pivots) so we start at i + for colidx in i:size(A, 2) + A[pos, colidx] *= w + end + # A[pos, i] is 1 now + + # zero the whole i-th column above and below pivot: + for rowidx = 1:size(A, 1) + rowidx == pos && continue # don't zero the active row + v = A[rowidx, i] + iszero(v) && continue + for k in i:size(A, 2) # to the left of pivot everything is zero + A[rowidx, k] -= v * A[pos, k] + end + end + end + if sgn < 0 + A .*= sgn + end + return A, pivots +end + +row_echelon_form(A::AbstractMatrix) = row_echelon_form!(deepcopy(A)) diff --git a/src/eigenspacedecomposition.jl b/src/Characters/eigenspacedecomposition.jl similarity index 53% rename from src/eigenspacedecomposition.jl rename to src/Characters/eigenspacedecomposition.jl index 773ad5f..33aca60 100644 --- a/src/eigenspacedecomposition.jl +++ b/src/Characters/eigenspacedecomposition.jl @@ -1,109 +1,3 @@ -# Version over exact field: -function row_echelon_form!(A::AbstractMatrix{T}) where {T<:Union{FiniteFields.GF, Cyclotomic{<:Rational}, Rational}} - pivots = Int[] - pos = 0 - for i = 1:size(A, 2) - j = findnext(x -> !iszero(x), @view(A[:, i]), pos+1) - j === nothing && continue - pos += 1 - push!(pivots, i) - - # swap the rows so that - if (pos != j) - for k in 1:size(A, 2) - A[pos, k], A[j, k] = A[j, k], A[pos, k] - end - end - # A[pos, :] is the row with leading nz in i-th column - - w = inv(A[pos, i]) - # all columns to the left of i are zeroed (below pivots) so we start at i - for colidx in i:size(A, 2) - A[pos, colidx] *= w - end - # A[pos, i] is 1 now - - # zero the whole i-th column above and below pivot: - for rowidx = 1:size(A, 1) - rowidx == pos && continue # don't zero the active row - v = A[rowidx, i] - iszero(v) && continue - for k in i:size(A, 2) # to the left of pivot everything is zero - A[rowidx, k] -= v * A[pos, k] - end - end - end - return A, pivots -end - -function row_echelon_form!(A::AbstractMatrix{C}) where {T<:AbstractFloat, C<:Cyclotomic{T}} - pivots = Int[] - pos = 0 - for i = 1:size(A, 2) - - j = findnext(x -> abs(x) > sqrt(eps(T)), @view(A[:, i]), pos+1) - j === nothing && continue - pos += 1 - push!(pivots, i) - - if (pos != j) - for k in 1:size(A, 2) - A[pos, k], A[j, k] = A[j, k], A[pos, k] - end - end - - @assert abs(A[pos, i]) >= sqrt(eps(T)) - w = inv(A[pos, i]) - - for colidx in i:size(A, 2) - if abs(A[pos, colidx]) <= sqrt(eps(T)) - Cyclotomics.zero!(A[pos, colidx]) - elseif colidx == i - Cyclotomics.one!(A[pos, colidx]) - else - A[pos, colidx] *= w - end - end - - for rowidx = 1:size(A, 1) - rowidx == pos && continue - v = A[rowidx, i] - if abs(v) < sqrt(eps(T)) - Cyclotomics.zero!(v) - else - for k in i:size(A, 2) - A[rowidx, k] -= v * A[pos, k] - end - end - end - end - - for i in eachindex(A) - ndigs = floor(Int, -log10(eps(T)) - log10(max(size(A)...)) - 3) - A[i] = Cyclotomics.roundcoeffs!(A[i], digits=ndigs) - end - - return A, pivots -end - -row_echelon_form(A::AbstractMatrix) = row_echelon_form!(deepcopy(A)) - -image_basis(A::AbstractMatrix) = image_basis!(deepcopy(A)) - -image_basis!(A::AbstractMatrix) = row_echelon_form!(A) - -function image_basis!(A::AbstractMatrix{T}) where {T<:AbstractFloat} - fact = svd!(A) - A_rank = sum(fact.S .> maximum(size(A)) * eps(T)) - return fact.Vt, 1:A_rank -end - -function image_basis!(A::AbstractMatrix{T}) where {T<:Complex} - fact = svd!(A) - A_rank = sum(fact.S .> maximum(size(A)) * 2eps(real(T))) - return fact.U, 1:A_rank -end - function right_nullspace(M::AbstractMatrix{T}) where {T} A, pivots = row_echelon_form(M) ncolsA = size(A, 2) @@ -124,23 +18,6 @@ function left_nullspace(M::AbstractMatrix) return transpose(right_nullspace(transpose(M))) end -# function left_nullspace(M::AbstractMatrix{T}) where {T} -# At, pivots = row_echelon_form(transpose(M)) -# A = transpose(A) -# nrows = size(A, 1) -# length(pivots) == ncolsA && return zeros(T, 0,0) -# W = zeros(T, nrowsA - length(pivots), nrowsA) -# for (i, el) in enumerate(setdiff(1:rowsA, pivots)) -# W[i, el] += 1 -# for (j, k) in enumerate(pivots) -# if j < el -# W[i,k] -= A[el, j] -# end -# end -# end -# return W -# end - function left_eigen(M::AbstractMatrix{T}) where {T<:FiniteFields.GF} @assert ==(size(M)...) eigen = Dict{T,typeof(M)}() diff --git a/src/gf.jl b/src/Characters/gf.jl similarity index 98% rename from src/gf.jl rename to src/Characters/gf.jl index 9184e9f..89bcfb1 100644 --- a/src/gf.jl +++ b/src/Characters/gf.jl @@ -32,6 +32,7 @@ Base.hash(n::GF{q}, h::UInt) where {q} = Base.:+(n::GF{q}, m::GF{q}) where {q} = GF{q}(Int(n) + Int(m), false) Base.:-(n::GF{q}, m::GF{q}) where {q} = GF{q}(Int(n) - Int(m), false) Base.:*(n::GF{q}, m::GF{q}) where {q} = GF{q}(Int(n) * Int(m), false) +Base.:*(n::Integer, m::GF{q}) where {q} = GF{q}(n * Int(m), false) Base.:/(n::GF{q}, m::GF{q}) where {q} = n * inv(m) Base.:-(n::GF{q}) where {q} = GF{q}(q - Int(n), false) diff --git a/src/powermap.jl b/src/Characters/powermap.jl similarity index 100% rename from src/powermap.jl rename to src/Characters/powermap.jl diff --git a/src/SymbolicWedderburn.jl b/src/SymbolicWedderburn.jl index fb4b4f7..e3ca167 100644 --- a/src/SymbolicWedderburn.jl +++ b/src/SymbolicWedderburn.jl @@ -1,28 +1,41 @@ module SymbolicWedderburn using LinearAlgebra +using SparseArrays using Primes using Cyclotomics using GroupsCore - -import PermutationGroups: - AbstractOrbit, AbstractPerm, AbstractPermutationGroup, Orbit, Perm, degree -import PermutationGroups - -export conjugacy_classes, symmetry_adapted_basis - -include("gf.jl") -include("eigenspacedecomposition.jl") -include("cmmatrix.jl") - -include("powermap.jl") -include("characters.jl") -include("characters_arith.jl") -include("dixon.jl") +using PermutationGroups +using StarAlgebras + +export symmetry_adapted_basis, WedderburnDecomposition +export basis, + # degree, # too common name to export + direct_summands, + invariant_vectors, + issimple, + multiplicity + +macro spawn_compat(expr) + @static if VERSION < v"1.3.0" + return :(@async $(esc(expr))) + else + return :(Threads.@spawn $(esc(expr))) + end +end + +include("Characters/Characters.jl") +using .Characters +import .Characters: row_echelon_form! +import .Characters.FiniteFields include("actions.jl") -include("projections.jl") +include("action_characters.jl") +include("matrix_projections.jl") +include("minimal_projections.jl") +include("direct_summands.jl") include("sa_basis.jl") +include("wedderburn_decomposition.jl") end # module diff --git a/src/action_characters.jl b/src/action_characters.jl new file mode 100644 index 0000000..a9522e0 --- /dev/null +++ b/src/action_characters.jl @@ -0,0 +1,63 @@ +## characters defined by actions/homomorphisms + +function _action_class_fun( + conjugacy_cls::AbstractVector{CCl}, +) where {CCl <: AbstractOrbit{<:PermutationGroups.AbstractPerm}} + vals = Int[PermutationGroups.nfixedpoints(first(cc)) for cc in conjugacy_cls] + # in general: + # vals = [tr(matrix_representative(first(cc))) for cc in conjugacy_cls] + return Characters.ClassFunction(vals, conjugacy_cls) +end + +function _action_class_fun( + conjugacy_cls::AbstractVector{CCl}, +) where {CCl <: AbstractOrbit{<:AbstractMatrix}} + vals = [tr(first(cc)) for cc in conjugacy_cls] + return Characters.ClassFunction(vals, conjugacy_cls) +end + +function _action_class_fun( + hom::InducedActionHomomorphism{<:ByPermutations}, + conjugacy_cls +) + vals = Int[PermutationGroups.nfixedpoints(induce(hom, first(cc))) for cc in conjugacy_cls] + # in general: + # vals = [tr(matrix_representative(first(cc))) for cc in conjugacy_cls] + return Characters.ClassFunction(vals, conjugacy_cls) +end + +function _action_class_fun( + hom::InducedActionHomomorphism{<:ByLinearTransformation}, + conjugacy_cls, +) + vals = [tr(induce(hom, first(cc))) for cc in conjugacy_cls] + return Characters.ClassFunction(vals, conjugacy_cls) +end + +""" + action_character(conjugacy_cls, tbl::CharacterTable) +Return the character of the representation given by the elements in the +conjugacy classes `conjugacy_cls`. + +This corresponds to the classical definition of characters as a traces of the +representation matrices. +""" +function action_character(conjugacy_cls, tbl::CharacterTable) + ac_char = _action_class_fun(conjugacy_cls) + constituents = Int[dot(ac_char, χ) for χ in irreducible_characters(tbl)] + return Character(tbl, constituents) +end + +""" + action_character(hom::InducedActionHomomorphism, tbl::CharacterTable) +Return the character of the representation given by the images under `hom` of +elements in `conjugacy_classes(tbl)`. + +This corresponds to the classical definition of characters as a traces of the +representation matrices. +""" +function action_character(hom::InducedActionHomomorphism, tbl::CharacterTable) + ac_char = _action_class_fun(hom, conjugacy_classes(tbl)) + constituents = Int[dot(ac_char, χ) for χ in irreducible_characters(tbl)] + return Character(tbl, constituents) +end diff --git a/src/actions.jl b/src/actions.jl index a20003c..750e55a 100644 --- a/src/actions.jl +++ b/src/actions.jl @@ -2,7 +2,13 @@ abstract type Action end abstract type ByPermutations <: Action end abstract type ByLinearTransformation <: Action end -abstract type InducedActionHomomorphism{A, T} end +abstract type InducedActionHomomorphism{A,T} end +#= +implements: +* features(action_hom) +* action(hom::InducedActionHomomorphism) → Action +* action(hom::InducedActionHomomorphism, g::GroupElement, feature) → the action of `g` on `feature` +=# function induce(ac::Action, hom::InducedActionHomomorphism, g::GroupElement) throw( @@ -10,86 +16,20 @@ function induce(ac::Action, hom::InducedActionHomomorphism, g::GroupElement) ) end +induce(hom::InducedActionHomomorphism, g::GroupElement) = induce(action(hom), hom, g) + decompose(x::Any, hom::InducedActionHomomorphism) = throw( "No fallback is provided for $(typeof(x)). You need to implement `decompose(::$(typeof(x)), ::$(typeof(hom)))`.", ) coeff_type(ac::ByLinearTransformation) = throw( - "No fallback is provided for $(typeof(ac)). You need to implement `coeff_type(::$(typeof(ac)))`.", + "No fallback is provided for $(typeof(ac)). You need to implement `_coeff_type(::$(typeof(ac)))`.", ) -#= -implements: -* features(action_hom) -* action(hom::InducedActionHomomorphism) → Action -* action(hom::InducedActionHomomorphism, g::GroupElement, feature) → return action of `g` on `feature` -=# - -struct ExtensionHomomorphism{A,T,I,V} <: InducedActionHomomorphism{A, T} - ac::A - features::V # supports linear indexing - reversef::Dict{T,I} -end - -_inttype(::ExtensionHomomorphism{A,T,I}) where {A,T,I} = I - -@static if VERSION < v"1.3.0" - # cannot add methods to an abstract type in pre Julia v1.3 - (hom::ExtensionHomomorphism)(g::GroupElement) = induce(action(hom), hom, g) - - function (hom::ExtensionHomomorphism)(orb::O) where {T,O<:AbstractOrbit{T,Nothing}} - elts = hom.(orb) - dict = Dict(e => nothing for e in elts) - return Orbit(elts, dict) - end - - function (hom::ExtensionHomomorphism)(chars::AbstractArray{<:Character}) - χ = first(chars) - iccG = hom.(conjugacy_classes(χ)) - return [Character(values(χ), χ.inv_of, iccG) for χ in chars] - end -else - (hom::InducedActionHomomorphism)(g::GroupElement) = induce(action(hom), hom, g) - - function (hom::InducedActionHomomorphism)(orb::O) where {T,O<:AbstractOrbit{T,Nothing}} - S = typeof(hom(one(first(orb)))) - elts = Vector{S}(undef, length(orb)) - for i in 1:length(orb) - elts[i] = hom(orb.elts[i]) - end - return Orbit(elts) - end - - function (hom::InducedActionHomomorphism)(chars::AbstractArray{<:Character}) - χ = first(chars) - iccG = hom.(conjugacy_classes(χ)) - return [Character(values(χ), χ.inv_of, iccG) for χ in chars] - end -end - -function ExtensionHomomorphism{I}(ac::Action, features) where {I<:Integer} - @assert typemax(I) >= length(features) "Use wider Integer type!" - reversef = Dict{eltype(features),I}(f => idx for (idx, f) in pairs(features)) - return ExtensionHomomorphism(ac, features, reversef) -end - -# Exceeding typemax(UInt16) here would mean e.g. that you're trying to block-diagonalize SDP constraint of size 65535×65535, which is highly unlikely ;) -ExtensionHomomorphism(ac::Action, features) = ExtensionHomomorphism{UInt16}(ac, features) - -# interface: -features(hom::ExtensionHomomorphism) = hom.features -action(hom::ExtensionHomomorphism) = hom.ac +# Exceeding typemax(UInt16) here would mean e.g. that you're trying to block-diagonalize +# an SDP constraint of size 65535×65535, which is highly unlikely ;) +_int_type(::InducedActionHomomorphism) = UInt16 -# user convenience functions: -Base.getindex(ehom::ExtensionHomomorphism, i::Integer) = ehom.features[i] -Base.getindex(ehom::ExtensionHomomorphism{A,T}, f::T) where {A,T} = ehom.reversef[f] - -function induce(::ByPermutations, hom::ExtensionHomomorphism, g::GroupElement) - I = _inttype(hom) - return Perm(vec(I[hom[action(action(hom), g, f)] for f in features(hom)])) -end - -using SparseArrays function induce(ac::ByLinearTransformation, hom::InducedActionHomomorphism, g::GroupElement) n = length(features(hom)) @@ -109,21 +49,67 @@ function induce(ac::ByLinearTransformation, hom::InducedActionHomomorphism, g::G end # ala permutation action -function decompose(m::T, hom::InducedActionHomomorphism{A, T}) where {A, T} +function decompose(m::T, hom::InducedActionHomomorphism{A,T}) where {A,T} return [hom[m]], [one(coeff_type(action(hom)))] end -struct CachedExtensionHomomorphism{A, T, G, H, E<:InducedActionHomomorphism{A, T}} <: InducedActionHomomorphism{A,T} + +struct ExtensionHomomorphism{A,T,I,V} <: InducedActionHomomorphism{A,T} + ac::A + features::V # supports linear indexing + reversef::Dict{T,I} +end + +# see the comment about UInt16 above +ExtensionHomomorphism(ac::Action, features) = ExtensionHomomorphism{UInt16}(ac, features) + +function ExtensionHomomorphism{I}(ac::Action, features) where {I <: Integer} + @assert typemax(I) >= length(features) "Use wider Integer type!" + reversef = Dict{eltype(features),I}(f => idx for (idx, f) in pairs(features)) + return ExtensionHomomorphism(ac, features, reversef) +end + +_int_type(::ExtensionHomomorphism{A,T,I}) where {A,T,I} = I + +Base.getindex(ehom::ExtensionHomomorphism, i::Integer) = ehom.features[i] +Base.getindex(ehom::ExtensionHomomorphism{A,T}, f::T) where {A,T} = ehom.reversef[f] + +# interface: +features(hom::ExtensionHomomorphism) = hom.features +action(hom::ExtensionHomomorphism) = hom.ac + +function induce(::ByPermutations, hom::ExtensionHomomorphism, g::GroupElement) + I = _int_type(hom) + return Perm(vec(I[hom[action(action(hom), g, f)] for f in features(hom)])) +end + +struct CachedExtensionHomomorphism{A,T,G,H,E <: InducedActionHomomorphism{A,T}} <: InducedActionHomomorphism{A,T} ehom::E - cache::Dict{G, H} + cache::Dict{G,H} end -for f in (:_inttype, :features, :action) +for f in (:action, :features, :_int_type) @eval $f(h::CachedExtensionHomomorphism) = $f(h.ehom) end -CachedExtensionHomomorphism{G, H}(hom::InducedActionHomomorphism) where {G,H} = - CachedExtensionHomomorphism(hom, Dict{G, H}()) +Base.getindex(h::CachedExtensionHomomorphism, x::Any) = h.ehom[x] + +function CachedExtensionHomomorphism(G::Group, action::Action, basis; precompute=false) + hom = ExtensionHomomorphism(action, basis) + S = typeof(induce(hom, one(G))) + chom = CachedExtensionHomomorphism{eltype(G), S}(hom) + if precompute + # to make sure that the future access to chom is read-only, i.e. thread-safe + # one may choose to precompute the images + for g in G + induce(chom, g) + end + end + return chom +end + +CachedExtensionHomomorphism{G,H}(hom::InducedActionHomomorphism) where {G,H} = + CachedExtensionHomomorphism(hom, Dict{G,H}()) function induce(ac::Action, hom::CachedExtensionHomomorphism, g::GroupElement) if !haskey(hom.cache, g) @@ -131,3 +117,10 @@ function induce(ac::Action, hom::CachedExtensionHomomorphism, g::GroupElement) end return hom.cache[g] end +#disambiguation: +function induce(ac::ByLinearTransformation, hom::CachedExtensionHomomorphism, g::GroupElement) + if !haskey(hom.cache, g) + hom.cache[g] = induce(ac, hom.ehom, g) + end + return hom.cache[g] +end diff --git a/src/characters.jl b/src/characters.jl deleted file mode 100644 index bbd2986..0000000 --- a/src/characters.jl +++ /dev/null @@ -1,271 +0,0 @@ -""" - AbstractClassFunction -Abstract type representing functions constant on conjugacy classes of a group (e.g. characters). - -Subtypes need to implement: - * `getindex(χ, i::Integer)`: return the value on `i`-th conjugacy class. - * Indexing with negative integers should return values on the class which contains inverses of the `i`-th class - * `χ(g::GroupElem)`: return value of `χ` on a group element (only for julia < 1.3) - * `values(χ)`: return an iterator over values - * `conjugacy_classes(χ)`: return iterator over conjugacy classes of `χ`. - * `conj(χ)`: return the conjugate function - -It is assumed that two class functions on the same group will return **identical** (ie. `===`) conjugacy_classes. -""" -abstract type AbstractClassFunction{T} end # <: AbstractVector{T} ?? - -Base.eltype(::Type{<:AbstractClassFunction{T}}) where {T} = T - -function LinearAlgebra.dot(χ::AbstractClassFunction, ψ::AbstractClassFunction) - val = sum( - length(cc) * χ[i] * ψ[-i] for - (i, cc) in enumerate(conjugacy_classes(χ)) - ) - orderG = sum(length, conjugacy_classes(χ)) - val = _div(val, orderG) - return val -end - -_div(val, orderG) = div(val, orderG) -_div(val::AbstractFloat, orderG) = val/orderG -_div(val::ComplexF64, orderG) = val/orderG - -Base.:(==)(χ::AbstractClassFunction, ψ::AbstractClassFunction) = - conjugacy_classes(χ) === conjugacy_classes(ψ) && values(χ) == values(ψ) - -Base.hash(χ::AbstractClassFunction, h::UInt) = - hash(conjugacy_classes(χ), hash(values(χ), hash(AbstractClassFunction, h))) - -Base.parent(χ::AbstractClassFunction) = parent(first(first(conjugacy_classes(χ)))) - -#################################### -# Characters - -struct VirtualCharacter{T,CCl} <: AbstractClassFunction{T} - vals::Vector{T} - inv_of::Vector{Int} - cc::Vector{CCl} -end - -struct Character{T,CCl} <: AbstractClassFunction{T} - vals::Vector{T} - inv_of::Vector{Int} - cc::Vector{CCl} -end - -const ClassFunction = Union{Character,VirtualCharacter} - -""" - affordable_real!(χ::ClassFunction[, pmap::PowerMap]) -Return either `χ` or `2re(χ)` depending whether `χ` is afforded by a real representation, modifying `χ` in place. -""" -function affordable_real!( - χ::ClassFunction, - pmap = PowerMap(conjugacy_classes(χ)), -) - ι = frobenius_schur_indicator(χ, pmap) - if ι <= 0 # i.e. χ is complex or quaternionic - χ.vals .+= values(conj(χ)) - return χ - else - return χ - end -end - -""" - frobenius_schur_indicator(χ::AbstractClassFunction[, pmap::PowerMap]) -Return Frobenius-Schur indicator of `χ`, i.e. `Σχ(g²)` where sum is taken over -the whole group. - -If χ is an irreducible `Character`, Frobenius-Schur indicator takes values in -`{1, 0, -1}` which correspond to the following situations: - 1. `χ` is real-valued and is afforded by an irreducible real representation, - 2. `χ` is a complex character which is not afforded by a real representation, and - 3. `χ` is quaternionic character, i.e. it is real valued, but is not afforded a real representation. - -In cases 2. and 3. `2re(χ) = χ + conj(χ)` corresponds to an irreducible character -afforded by a real representation. -""" -function frobenius_schur_indicator( - χ::AbstractClassFunction, - pmap::PowerMap = PowerMap(conjugacy_classes(χ)), -) - ι = sum( - length(c) * χ[pmap[i, 2]] for (i, c) in enumerate(conjugacy_classes(χ)) - ) - - ι_int = Int(ι) - ordG = sum(length, conjugacy_classes(χ)) - d, r = divrem(ι_int, ordG) - @assert r == 0 "Non integral Frobenius Schur Indicator: $(ι_int) = $d * $ordG + $r" - return d -end - -Base.isreal(χ::AbstractClassFunction) = frobenius_schur_indicator(χ) > 0 - -if VERSION >= v"1.3.0" - function (χ::AbstractClassFunction)(g::GroupsCore.GroupElement) - for (i, cc) in enumerate(conjugacy_classes(χ)) - g ∈ cc && return χ[i] - end - throw( - DomainError( - g, - "element does not belong to conjugacy classes of $χ", - ), - ) - end -end - -################################### -# Specific definitions for ClassFunction -const ClassFunction = Union{Character,VirtualCharacter} - -function Character( - vals::AbstractVector{T}, - ccls::AbstractVector{CCl}, -) where {T,CCl} - χ = Character{T,CCl}(vals, _inv_of(ccls), ccls) - return χ -end - -function VirtualCharacter( - vals::AbstractVector{T}, - ccls::AbstractVector{CCl}, -) where {T,CCl} - χ = VirtualCharacter{T,CCl}(vals, _inv_of(ccls), ccls) - return χ -end - -function Base.deepcopy_internal(χ::CF, ::IdDict) where {CF<:ClassFunction} - return CF(copy(χ.vals), χ.inv_of, conjugacy_classes(χ)) -end - -Character{S}(χ::Character{T,Cl}) where {S,T,Cl} = - Character{S, Cl}(values(χ), χ.inv_of, conjugacy_classes(χ)) - -VirtualCharacter(χ::Character{T,Cl}) where {T,Cl} = VirtualCharacter{T,Cl}(χ) -VirtualCharacter{T}(χ::Character{S,Cl}) where {T,S,Cl} = - VirtualCharacter{T,Cl}(χ) -VirtualCharacter(χ::VirtualCharacter) = deepcopy(χ) - -VirtualCharacter{T,Cl}(χ::ClassFunction) where {T,S,Cl} = - VirtualCharacter{T,Cl}(values(χ), χ.inv_of, conjugacy_classes(χ)) - -Base.@propagate_inbounds function Base.getindex(χ::ClassFunction, i::Integer) - @boundscheck checkbounds(values(χ), abs(i)) - if i < 0 - return @inbounds values(χ)[χ.inv_of[abs(i)]] - else - return @inbounds values(χ)[i] - end -end - -if VERSION < v"1.3.0" - function (χ::Character)(g::GroupsCore.GroupElement) - for (i, cc) in enumerate(conjugacy_classes(χ)) - g ∈ cc && return χ[i] - end - throw( - DomainError( - g, - "element does not belong to conjugacy classes of $χ", - ), - ) - end - function (χ::VirtualCharacter)(g::GroupsCore.GroupElement) - for (i, cc) in enumerate(conjugacy_classes(χ)) - g ∈ cc && return χ[i] - end - throw( - DomainError( - g, - "element does not belong to conjugacy classes of $χ", - ), - ) - end -end - -Base.values(χ::ClassFunction) = χ.vals -conjugacy_classes(χ::ClassFunction) = χ.cc - -Base.conj(χ::Cf) where {Cf<:ClassFunction} = - Cf(values(χ)[χ.inv_of], χ.inv_of, conjugacy_classes(χ)) - -PermutationGroups.degree(χ::Character) = Int(χ(one(parent(χ)))) -PermutationGroups.degree( - χ::Character{T,CCl}, -) where {T,CCl<:AbstractOrbit{<:AbstractMatrix}} = Int(χ[1]) - -function _inv_of(cc::AbstractVector{<:AbstractOrbit}) - inv_of = zeros(Int, size(cc)) - for (i, c) in enumerate(cc) - g = inv(first(c)) - inv_of[i] = something(findfirst(k -> g in k, cc), 0) - end - any(iszero, inv_of) && throw( - ArgumentError( - "Could not find the conjugacy class for inverse of $(first(cc[findfirst(iszero, inv_of)])).", - ), - ) - return inv_of -end - -function normalize!(χ::Character) - - id = one(parent(χ)) - - k = χ(id) - if !isone(k) - χ.vals .*= inv(k) - end - - # ⟨χ, χ⟩ = 1/d² - - deg = sqrt(inv(dot(χ, χ))) - @debug "normalizing with" n dot(χ, χ) χ(id) χ - - # normalizing χ - χ.vals .*= deg - - return χ -end - -for C in (:Character, :VirtualCharacter) - @eval begin - function Base.show(io::IO, ::MIME"text/plain", χ::$C) - println(io, string($C) * " over $(eltype(χ))") - cc_reps = string.(first.(conjugacy_classes(χ))) - k = maximum(length.(cc_reps)) - - vals = values(χ) - for i = 1:length(vals) - (c, v) = cc_reps[i], vals[i] - if i == length(values(χ)) - print(io, rpad("$c^G", k + 3), "→ \t", v) - else - print(io, rpad("$c^G", k + 3), "→ \t", v, "\n") - end - end - end - - function Base.show(io::IO, χ::$C) - v = values(χ) - io = IOContext(io, :typeinfo => eltype(v)) - limited = get(io, :limit, false) - opn, cls = '[', ']' - - if limited && length(v) > 20 - Base.show_delim_array(io, v, opn, ",", "", false, 1, 9) - print(io, " … ") - l = length(v) - Base.show_delim_array(io, v, "", ",", cls, false, l - 9, l) - else - Base.show_delim_array(io, v, opn, ",", cls, false) - end - end - - Base.show(io::IO, ::Type{<:$C{T}}) where T = - print(io, $C, "{", T, ", …}") - end -end diff --git a/src/characters_arith.jl b/src/characters_arith.jl deleted file mode 100644 index 13d16ef..0000000 --- a/src/characters_arith.jl +++ /dev/null @@ -1,31 +0,0 @@ -Base.promote_rule( - ::Type{Character{T,Cl}}, - ::Type{VirtualCharacter{S,Cl}}, -) where {T,S,Cl} = VirtualCharacter{promote_type(T, S),Cl} - -for f in (:+, :-) - @eval begin - function Base.$f(χ::ClassFunction, ψ::ClassFunction) - @assert χ.inv_of == ψ.inv_of - @assert conjugacy_classes(χ) === conjugacy_classes(ψ) - return VirtualCharacter( - $f.(values(χ), values(ψ)), - χ.inv_of, - conjugacy_classes(χ), - ) - end - end -end - -for C in (:Character, :VirtualCharacter) - @eval begin - Base.:*(χ::$C, c::Number) = - VirtualCharacter(c .* values(χ), χ.inv_of, conjugacy_classes(χ)) - - Base.:*(c::Number, χ::$C) = - VirtualCharacter(c .* values(χ), χ.inv_of, conjugacy_classes(χ)) - - Base.:/(χ::$C, c::Number) = - VirtualCharacter(values(χ) ./ c, χ.inv_of, conjugacy_classes(χ)) - end -end diff --git a/src/direct_summands.jl b/src/direct_summands.jl new file mode 100644 index 0000000..87be901 --- /dev/null +++ b/src/direct_summands.jl @@ -0,0 +1,28 @@ +struct DirectSummand{T, M<:AbstractMatrix{T}} <: AbstractMatrix{T} + basis::M + multiplicity::Int + degree::Int + simple::Bool +end + +image_basis(ds::DirectSummand) = ds.basis +issimple(ds::DirectSummand) = ds.simple +PermutationGroups.degree(ds::DirectSummand) = ds.degree +multiplicity(ds::DirectSummand) = ds.multiplicity + +Base.size(ds::DirectSummand) = size(ds.basis) +Base.@propagate_inbounds Base.getindex(ds::DirectSummand, i...) = + image_basis(ds)[i...] + +function SparseArrays.sparse(ds::DirectSummand) + sp = sparse(image_basis(ds)) + return DirectSummand(sp, multiplicity(ds), degree(ds), issimple(ds)) +end + +function SparseArrays.droptol!( + ds::DirectSummand{T, <:AbstractSparseArray}, + tol +) where T + droptol!(image_basis(ds), tol) + return ds +end diff --git a/src/dixon.jl b/src/dixon.jl deleted file mode 100644 index 58ffd81..0000000 --- a/src/dixon.jl +++ /dev/null @@ -1,120 +0,0 @@ -Base.exponent(G::GroupsCore.Group) = exponent(conjugacy_classes(G)) -Base.exponent(cclasses::AbstractArray) = lcm(order.(first.(cclasses))) -dixon_prime(G::GroupsCore.Group) = dixon_prime(order(G), exponent(G)) - -function dixon_prime(cclasses::AbstractArray) - ordG = sum(length, cclasses) - m = exponent(cclasses) - return dixon_prime(ordG, m) -end - -function dixon_prime(ordG::Integer, exponent::Integer) - # `FiniteFields.GF{p}` needs `p` to be `Int` so no need to do the - # computations of this function over `BigInt` if `ordG` is so we convert - # to `Int`. - p = 2 * convert(Int, isqrt(ordG)) - while true - p = nextprime(p + 1) - isone(p % exponent) && break # we need -1 to be in the field - end - return p -end - -function common_esd(Ns, F::Type{<:FiniteFields.GF}) - itr = iterate(Ns) - @assert itr !== nothing - esd = EigenSpaceDecomposition(F.(first(itr))) - for N in Iterators.rest(Ns, last(itr)) - esd = refine(esd, F.(N)) - @debug N esd.eigspace_ptrs - isdiag(esd) && return esd - end - return esd -end - -characters_dixon(G::GroupsCore.Group) = - characters_dixon(Rational{Int}, G) -characters_dixon(::Type{R}, G::GroupsCore.Group) where {R<:Real} = - characters_dixon(R, conjugacy_classes(G)) - -function characters_dixon(::Type{R}, cclasses::AbstractVector) where {R} - p = dixon_prime(cclasses) - chars_𝔽p = characters_dixon(FiniteFields.GF{p}, cclasses) - return complex_characters(R, chars_𝔽p) -end - -function characters_dixon(F::Type{<:FiniteFields.GF}, cclasses::AbstractVector) - Ns = [CMMatrix(cclasses, i) for i = 1:length(cclasses)] - esd = common_esd(Ns, F) - @assert isdiag(esd) "Class Matrices failed to diagonalize! $esd" - inv_ccls = _inv_of(cclasses) - return [ - normalize!(Character(vec(eigensubspace), inv_ccls, cclasses)) for - eigensubspace in esd - ] -end - -function _multiplicities( - chars::AbstractVector{<:Character{F}}, - cclasses = conjugacy_classes(first(chars)), -) where {F<:FiniteFields.GF} - - e = Int(exponent(cclasses)) - ie = inv(F(e)) - - ω⁻¹ = inv(FiniteFields.rootofunity(F, e)) - ωs = Matrix{typeof(ω⁻¹)}(undef, e, e) - Threads.@threads for k in 0:e-1 - ω⁻ᵏ = ω⁻¹^k - for l in 0:e-1 - ωs[l+1, k+1] = ω⁻ᵏ^l - end - end - - multiplicities = zeros(Int, length(chars), length(cclasses), e) - powermap = PowerMap(cclasses) - - for (i, χ) in enumerate(chars) - Threads.@threads for j = 1:length(cclasses) - for k = 0:e-1 - val = Int(ie * sum(χ[powermap[j, l]] * ωs[l+1, k+1] for l = 0:e-1)) - multiplicities[i, j, k+1] = val - end - end - end - - return multiplicities -end - -function complex_characters( - ::Type{R}, - chars::AbstractVector{<:Character{F,CCl}}, -) where {R,F<:FiniteFields.GF,CCl} - - cclasses = conjugacy_classes(first(chars)) - lccl = length(cclasses) - mult_c = _multiplicities(chars, cclasses) - e = size(mult_c, 3) # the exponent - - inv_of_cls = first(chars).inv_of - - C = Cyclotomics.Cyclotomic{R,Cyclotomics.SparseVector{R,Int}} - # C = typeof(Cyclotomics.E(5)) - - complex_chars = Vector{Character{C,CCl}}(undef, length(chars)) - - Es = [E(e, k) for k in 0:e-1] - Threads.@threads for i = 1:length(complex_chars) - complex_chars[i] = Character{C,CCl}( - [ - Cyclotomics.reduced_embedding( - sum(mult_c[i, j, k+1] * Es[k+1] for k = 0:e-1), - ) for j = 1:lccl - ], - inv_of_cls, - cclasses, - ) - end - - return complex_chars -end diff --git a/src/matrix_projections.jl b/src/matrix_projections.jl new file mode 100644 index 0000000..df18892 --- /dev/null +++ b/src/matrix_projections.jl @@ -0,0 +1,183 @@ +""" + matrix_projection_irr([hom::InducedActionHomomorphism, ]χ::Character) +Compute matrix projection associated to the irreducible character `χ`. + +If the homomorphism is not passed, the dimension `d` of the projection is +derived from `conjugacy_classes(χ)`. E.g. if `conjugacy_classes(χ)` consist +of permutations of degree `d` (i.e. acting naturally on the set `1:d`) the +result will be a matrix of size `(d,d)`. + +If the homomorphism is passed, the dimension will be derived in similar +manner from the elements of the image of the homomorphism. +""" +function matrix_projection_irr(χ::Character) + @assert isirreducible(χ) + mproj = matrix_projection_irr(collect(values(χ)), conjugacy_classes(χ)) + mproj .*= degree(χ) // sum(length, conjugacy_classes(χ)) # χ(1)/order(G) + return mproj +end + +function matrix_projection_irr( + vals, + ccls::AbstractVector{<:AbstractOrbit{<:PermutationGroups.AbstractPerm}}, +) + dim = degree(first(first(ccls))) + result = zeros(eltype(vals), dim, dim) + + for (val, cc) in zip(vals, ccls) + for g in cc + for i in 1:dim + result[i, i^g] += val + end + end + end + + return result +end + +function matrix_projection_irr( + vals, + ccls::AbstractVector{<:AbstractOrbit{<:AbstractMatrix}}, +) + # TODO: call to inv(Matrix(g)) is a dirty hack, since if `g` + # is given by a sparse matrix `inv(g)` will fail. + return sum(val .* sum(g -> inv(Matrix(g)), cc) for (val, cc) in zip(vals, ccls)) +end + +## versions with InducedActionHomomorphisms + +function matrix_projection_irr(hom::InducedActionHomomorphism, χ::Character{T}) where T + @assert isirreducible(χ) + mproj = matrix_projection_irr(hom, collect(values(χ)), conjugacy_classes(χ)) + mproj .*= degree(χ) // sum(length, conjugacy_classes(χ)) # χ(1)/order(G) + return eltype(mproj) == T ? mproj : T.(mproj) +end + +function matrix_projection_irr( + hom::InducedActionHomomorphism{<:ByPermutations}, + class_values, + conjugacy_cls, + ) + dim = degree(induce(hom, first(first(conjugacy_cls)))) + result = zeros(eltype(class_values), dim, dim) + for (val, ccl) in zip(class_values, conjugacy_cls) + for g in ccl + h = induce(hom, g) + for i in 1:dim + result[i, i^h] += val + end + end + end + return result +end + +function matrix_projection_irr( + hom::InducedActionHomomorphism{<:ByLinearTransformation}, + class_values, + conjugacy_cls, +) + return sum( + val .* sum(g -> induce(hom, inv(g)), cc) + for (val, cc) in zip(class_values, conjugacy_cls) + ) +end + +function matrix_projection(χ::Character{T}) where T + tbl = table(χ) + res = sum( + c .* matrix_projection_irr(ψ) + for (c, ψ) in zip(constituents(χ), irreducible_characters(tbl)) if !iszero(c) + ) + return eltype(res) == T ? res : T.(res) +end + +function matrix_projection( + hom::InducedActionHomomorphism{<:ByPermutations}, + α::AlgebraElement, + dim::Integer = length(features(hom)), +) + result = zeros(eltype(α), dim, dim) + b = basis(parent(α)) + + @inbounds for (j, a) in StarAlgebras._nzpairs(StarAlgebras.coeffs(α)) + g = induce(hom, b[j]) + for i in 1:dim + result[i, i^g] += a + end + end + + return result +end + +function matrix_projection( + hom::InducedActionHomomorphism{<:ByLinearTransformation}, + α::AlgebraElement, + dim::Integer = length(features(hom)), +) + + result = zeros(Base._return_type(*, Tuple{eltype(α), coeff_type(action(hom))}), dim, dim) + b = basis(parent(α)) + + for (j, a) in StarAlgebras._nzpairs(StarAlgebras.coeffs(α)) + result += a*induce(hom, inv(b[j])) + end + + return result +end + +function matrix_projection( + α::AlgebraElement{<:StarAlgebra{<:PermutationGroups.AbstractPermutationGroup}}, + dim::Integer = degree(parent(parent(α))) +) + result = zeros(eltype(α), dim, dim) + b = basis(parent(α)) + + @inbounds for (j, a) in StarAlgebras._nzpairs(StarAlgebras.coeffs(α)) + g = b[j] + for i in 1:dim + result[i, i^g] += a + end + end + + return result +end + +image_basis!(A::AbstractMatrix) = row_echelon_form!(A) + +function image_basis!(A::AbstractMatrix{T}) where {T<:AbstractFloat} + A, p = row_echelon_form!(A) + fact = svd!(Matrix(@view A[1:length(p), :])) + A_rank = sum(fact.S .> maximum(size(A)) * eps(T)) + return fact.Vt, 1:A_rank +end + +function image_basis!(A::AbstractMatrix{T}) where {T<:Complex} + fact = svd!(A) + A_rank = sum(fact.S .> maximum(size(A)) * 2eps(real(T))) + return fact.Vt, 1:A_rank +end + +image_basis(A::AbstractMatrix) = + ((m, p) = image_basis!(deepcopy(A)); m[1:length(p), :]) + +function image_basis(χ::Character) + mpr = isirreducible(χ) ? matrix_projection_irr(χ) : matrix_projection(χ) + image, pivots = image_basis!(mpr) + return image[1:length(pivots), :] +end + +function image_basis(hom::InducedActionHomomorphism, χ::Character) + mpr = isirreducible(χ) ? matrix_projection_irr(hom, χ) : matrix_projection(hom, χ) + image, pivots = image_basis!(mpr) + return image[1:length(pivots), :] +end + +function image_basis(α::AlgebraElement) + image, pivots = image_basis!(matrix_projection(α)) + return image[1:length(pivots), :] +end + +function image_basis(hom::InducedActionHomomorphism, α::AlgebraElement) + image, pivots = image_basis!(matrix_projection(hom, α)) + return image[1:length(pivots), :] +end diff --git a/src/minimal_projections.jl b/src/minimal_projections.jl new file mode 100644 index 0000000..c2c6c84 --- /dev/null +++ b/src/minimal_projections.jl @@ -0,0 +1,138 @@ +Base.parent(A::StarAlgebra{<:Group}) = A.object + +function StarAlgebras.AlgebraElement( + χ::AbstractClassFunction, + RG::StarAlgebra{<:Group}, +) + G = parent(RG) + @assert G === parent(χ) + b = basis(RG) + + dim = degree(χ) + ord = order(Int, G) + c = dim // ord + + T = Base._return_type(*, Tuple{typeof(c), eltype(χ)}) + + x = AlgebraElement(zeros(T, length(b)), RG) + + for (v, cc) in zip(values(χ), conjugacy_classes(χ)) + for g in cc + x[inv(g)] = c * v + end + end + return x +end + +function algebra_elt_from_support(support, RG::StarAlgebra) + b = basis(RG) + I = [b[s] for s in support] + V = fill(1, length(support)) + return AlgebraElement(sparsevec(I, V, length(b)), RG) +end +struct CyclicSubgroups{Gr, GEl} + group::Gr + seen::Dict{Int,Set{GEl}} + min_order::Int + max_order::Int + + function CyclicSubgroups(G::Group; min_order=1, max_order=order(Int, G)) + seen = Dict{Int, Set{eltype(G)}}() + return new{typeof(G), eltype(G)}(G, seen, min_order, max_order) + end +end + +Base.eltype(citr::CyclicSubgroups) = valtype(citr.seen) +Base.IteratorSize(::CyclicSubgroups) = Base.SizeUnknown() + +function Base.iterate(citr::CyclicSubgroups) + g, state = iterate(citr.group) # g is identity here + @assert isone(g) + if citr.min_order ≤ 1 ≤ citr.max_order + citr.seen[ord] = Set([g]) + return Set([g]), state + end + return iterate(citr, state) +end + +function Base.iterate(citr::CyclicSubgroups, state) + k = iterate(citr.group, state) + k === nothing && return nothing + g, state = k + ord = order(Int, g) + if citr.min_order ≤ ord ≤ citr.max_order + if !haskey(citr.seen, ord) + citr.seen[ord] = Set([g]) + return Set(g^i for i in 1:ord), state + else + if any(g^i in citr.seen[ord] for i in 1:ord-1) + return iterate(citr, state) + else + push!(citr.seen[ord], g) + return Set(g^i for i in 1:ord), state + end + end + end + return iterate(citr, state) +end + +function small_idempotents( + RG::StarAlgebra{<:Group}, + subgroups = CyclicSubgroups(parent(RG), min_order = 2), +) + return (algebra_elt_from_support(H, RG) // length(H) for H in subgroups) +end + +@static if VERSION < v"1.3.0" + function (χ::Character)(α::AlgebraElement{<:StarAlgebra{<:Group}}) + @assert parent(χ) === parent(parent(α)) + return sum(α(g) * χ(g) for g in supp(α)) + end +else + function (χ::AbstractClassFunction)(α::AlgebraElement{<:StarAlgebra{<:Group}}) + @assert parent(χ) === parent(parent(α)) + return sum(α(g) * χ(g) for g in supp(α)) + end +end + +function minimal_rank_projection(χ::Character, RG::StarAlgebra{<:Group}; + idems = small_idempotents(RG), # idems are AlgebraElements over Rational{Int} + iters=3 +) + + k = if isirreducible(χ) + 1 + else + sum(constituents(χ).>0) + # the minimal rank projection for the composite character + end + + degree(χ) == k && return one(RG, Rational{Int}), true + + for µ in idems + χ(µ) == k && µ^2 == µ && return µ, true + end + + for n in 2:iters + for elts in Iterators.product(ntuple(i -> idems, n)...) + µ = *(elts...) + χ(µ) == k && µ^2 == µ && return µ, true + end + end + @debug "Could not find minimal projection for $χ" + return one(RG, Rational{Int}), false +end + +function minimal_projection_system( + chars::AbstractVector{<:AbstractClassFunction}, + RG::StarAlgebra{<:Group} +) + res = fetch.([@spawn_compat minimal_rank_projection(χ, RG) for χ in chars]) + + r1p, simple = first.(res), last.(res) # rp1 are sparse storage + + mps = [isone(µ) ? AlgebraElement(χ, RG) : µ*AlgebraElement(χ, RG) + for (µ,χ) in zip(r1p, chars)] # dense storage + + return mps, simple +end diff --git a/src/projections.jl b/src/projections.jl deleted file mode 100644 index 4667865..0000000 --- a/src/projections.jl +++ /dev/null @@ -1,94 +0,0 @@ -""" - matrix_projection(χ::AbstractClassFunction) -Compute matrix projection associated to the abstract class function `χ`. - -The dimension `d` of the projection is derived from `conjugacy_classes(χ)`. -E.g. if `conjugacy_classes(χ)` consist of permutations of degree `d` (i.e. -acting naturally on the set `1:d`) the result will be a matrix of size `(d,d)`. - -Returned tuple consist of the matrix realization of the projection and a coefficient (its weight). -""" -function matrix_projection(χ::AbstractClassFunction{T}) where T - deg = degree(χ) - if iszero(deg) # short circuting the trivial case - return zeros(eltype(χ), 0, 0), zero(deg) - end - ordG = sum(length, conjugacy_classes(χ)) - U = matrix_projection(values(χ), conjugacy_classes(χ)) - - return U, T(deg) / ordG -end - -""" - matrix_projection(vals, ccls) -Return matrix projection defined by an abstract class function with values `vals` -attained on (the corresponding) conjugacy classes `ccls`. - -The dimension of the projection is equal to the degree of the permutations in `ccls`. -""" -function matrix_projection( - vals::AbstractVector{T}, - ccls::AbstractVector{<:AbstractOrbit{<:AbstractPerm}}, - dim = degree(first(first(ccls))), -) where {T} - - result = zeros(T, dim, dim) - - for (val, cc) in zip(vals, ccls) - for g in cc - for i in 1:dim - result[i, i^g] += val - end - end - end - - return result -end - -function matrix_projection( - vals::AbstractVector{T}, - ccls::AbstractVector{<:AbstractOrbit{<:AbstractMatrix}}, - dim = size(first(first(ccls)), 1), -) where {T} - - return sum(val .* sum(g -> inv(Matrix(g)), cc) for (val, cc) in zip(vals, ccls)) -end - -""" - action_character(conjugacy_cls) -Return the character of the representation given by the elements in the -conjugacy classes `conjugacy_cls`. - -This corresponds to the classical definition of characters as a traces of the -representation matrices. -""" -function action_character( - conjugacy_cls::AbstractVector{CCl}, - inv_of = _inv_of(conjugacy_cls), -) where {CCl<:AbstractOrbit{<:AbstractPerm}} - vals = Int[PermutationGroups.nfixedpoints(first(cc)) for cc in conjugacy_cls] - # in general: - # vals = [tr(matrix_representative(first(cc))) for cc in conjugacy_cls] - return Character(vals, inv_of, conjugacy_cls) -end - -function action_character( - conjugacy_cls::AbstractVector{CCl}, - inv_of = _inv_of(conjugacy_cls), -) where {CCl<:AbstractOrbit{<:AbstractMatrix}} - vals = [tr(first(cc)) for cc in conjugacy_cls] - return Character(vals, inv_of, conjugacy_cls) -end - -""" - affordable_real(chars::AbstractVector{<:AbstractClassFunction}) -Return _real_ characters formed from `chars` by replacing `χ` with `2re(χ)` when necessary. -""" -function affordable_real(chars::AbstractVector{<:AbstractClassFunction}) - pmap = PowerMap(conjugacy_classes(first(chars))) - return unique!(map(chars) do χ - χ = deepcopy(χ) - affordable_real!(χ, pmap) - end - ) -end diff --git a/src/sa_basis.jl b/src/sa_basis.jl index cf225d5..0ddb75a 100644 --- a/src/sa_basis.jl +++ b/src/sa_basis.jl @@ -1,160 +1,232 @@ -""" - isotypical_basis(χ::AbstractClassFunction) -Compute a basis of the image of the projection corresponding to a class function `χ`. - -Return the coefficients of basis vectors in an invariant subspace corresponding to `χ` -(so called _isotypical subspace_) in the action on `ℝⁿ` encoded by the conjugacy classes of `χ`. -""" -function isotypical_basis(χ::AbstractClassFunction) - u, weight = matrix_projection(χ) - image, pivots = if iszero(weight) - u_ = similar(u, 0, size(u, 2)) - image_basis!(u_) - else - u .*= weight - image_basis!(u) +function affordable_real( + irreducible_characters, + multiplicities=fill(1, length(irreducible_characters)), +) + irr_real = similar(irreducible_characters, 0) + mls_real = similar(multiplicities, 0) + for (i, χ) in pairs(irreducible_characters) + ι = Characters.frobenius_schur(χ) + if abs(ι) == 1 # real or quaternionic + @debug "real/quaternionic:" χ + push!(irr_real, χ) + push!(mls_real, multiplicities[i]) + else # complex one... + cχ = conj(χ) + k = findfirst(==(cχ), irreducible_characters) + @assert k !== nothing + @debug "complex" χ conj(χ)=irreducible_characters[k] + if k > i # ... we haven't already observed a conjugate of + @assert multiplicities[i] == multiplicities[k] + push!(irr_real, χ + cχ) + push!(mls_real, multiplicities[i]) + end + end end - @debug "isotypical subspace corresponding to χ has dimension $(length(pivots))" χ - return image[1:length(pivots), :] -end - -struct SemisimpleSummand{T} - basis::T - multiplicity::Int -end - -basis(b::SemisimpleSummand) = b.basis -Base.convert(::Type{M}, b::SemisimpleSummand) where {M<:AbstractMatrix} = - convert(M, basis(b)) - -multiplicity(b::SemisimpleSummand) = b.multiplicity - -function degree(b::SemisimpleSummand) - d, r = divrem(size(basis(b), 1), multiplicity(b)) - @assert iszero(r) - return d -end - -function extended_characters(::Type{T}, G::Group, basis, action) where {T} - chars = characters_dixon(T, G) - ehom = ExtensionHomomorphism(action, basis) - return ehom(chars) # result: Cyclotomic{T} valued Characters + return irr_real, mls_real end """ - symmetry_adapted_basis([T::Type,] G::AbstractPermutationGroup[, S=Rational{Int}]) -Compute a basis for the linear space `ℝⁿ` which is invariant under the symmetry of `G`. + symmetry_adapted_basis([T::Type,] G::AbstractPermutationGroup[, S=Rational{Int}; + semisimple=false]) +Compute a basis for the linear space `ℝⁿ` (where `n = degree(G)`) which is +invariant under the symmetry of `G`. -The permutation group is acting naturally on `1:degree(G)`. The coefficients of -the invariant basis are returned in (orthogonal) blocks corresponding to irreducible -characters of `G`. +The basis is returned as a vector of `DirectSummand{T}`s (blocks) corresponding +to the irreducible characters of `G`. The blocks are orthogonal to **each other**, +however vectors within a single block may *not* be orthogonal. +If `T<:LinearAlgebra.BlasFloat` BLAS routines will be used to orthogonalize +vectors within each `DirectSummand`. Arguments: * `S` controls the types of `Cyclotomic`s used in the computation of character table. Exact type are preferred. For larger groups `G` `Rational{BigInt}` might be necessary. * `T` controls the type of coefficients of the returned basis. +* `semisimple`: if set to `false` (the default) an effort to find minimal +projection system is made, so that each block defines a projection to a single +simple summand within each (isotypical) block. Otherwise an isomorphism to the +semisimple decomposition is computed. !!! Note: -Each block is invariant under the action of `G`, i.e. the action may permute -vectors from symmetry adapted basis within each block. The blocks are guaranteed -to be orthogonal. If `T<:LinearAlgebra.BlasFloat` BLAS routines will be used to -orthogonalize vectors within each block. +Each returned block (a `DirectSummand`) is invariant under the action of `G`, +which means that the action may still e.g. permute (row) vectors , but only +*within* each block. + +When `semisimple=true` the blocks constitute an isomorphism. Otherwise blocks +may represent only a projection onto the commutant of the appropriate matrix +algebra (which has in general lower dimension). This happens precisely when +`issimple` on a block returns `true` and `SymbolicWedderburn.degree(ds) == 1`. """ -symmetry_adapted_basis(G::AbstractPermutationGroup, S::Type = Rational{Int}) = - _symmetry_adapted_basis(characters_dixon(S, G)) +function symmetry_adapted_basis( + G::PermutationGroups.AbstractPermutationGroup, + S::Type = Rational{Int}; + semisimple::Bool = false, +) + tbl = CharacterTable(S, G) + return symmetry_adapted_basis(eltype(tbl), tbl, semisimple = semisimple) +end + +symmetry_adapted_basis( + T::Type, + G::PermutationGroups.AbstractPermutationGroup, + S::Type = Rational{Int}; + semisimple::Bool = false, +) = symmetry_adapted_basis(T, CharacterTable(S, G), semisimple = semisimple) + +function symmetry_adapted_basis( + T::Type, + tbl::CharacterTable; + semisimple::Bool = false, +) + irr, multips = _constituents_decomposition( + action_character(conjugacy_classes(tbl), tbl), + tbl, + ) + + if T <: Real + irr, multips = affordable_real(irr, multips) + end -symmetry_adapted_basis(T::Type, G::AbstractPermutationGroup, S::Type = Rational{Int}) = - symmetry_adapted_basis(T, characters_dixon(S, G)) + if semisimple || all(isone ∘ degree, irr) + return _symmetry_adapted_basis(T, irr, multips) + else + RG = _group_algebra(parent(tbl)) + return _symmetry_adapted_basis(T, irr, multips, RG) + end +end """ - symmetry_adapted_basis([T::Type,] G::Group, basis, action[, S=Rational{Int}]) -Compute a basis for the linear space spanned by `basis` which is invariant under -the symmetry of `G`. - -* The action used in these computations is -> `(b,g) → action(b,g)` for `b ∈ basis`, `g ∈ G` -and needs to be defined by the user. -* It is assumed that `G` acts on a subset of basis and the action needs to be + symmetry_adapted_basis([T::Type,] G::Group, action, basis[, S=Rational{Int}]; + semisimple=false) +Compute a decomposition of basis into (semi)simple subspaces which are invariant under +the action of `G`. + +It is assumed that `G` acts on a subset of basis and the action needs to be extended to the whole `basis`. If `G` is a permutation group already acting on the whole `basis`, a call to `symmetry_adapted_basis(G)` is preferred. * For inducing the action `basis` needs to be indexable and iterable (e.g. in the form of an `AbstractVector`). - -Arguments: -* `S` controlls the types of `Cyclotomic`s used in the computation of -character table. Exact type are preferred. For larger groups `G` `Rational{BigInt}` -might be necessary. -* `T` controls the type of coefficients of the returned basis. - -!!! Note: -Each block is invariant under the action of `G`, i.e. the action may permute -vectors from symmetry adapted basis within each block. The blocks are guaranteed -to be orthogonal. If `T<:LinearAlgebra.BlasFloat` BLAS routines will be used to -orthogonalize vectors within each block. """ -function symmetry_adapted_basis(G::Group, basis, action, S::Type = Rational{Int}) - chars_ext = extended_characters(S, G, basis, action) - return symmetry_adapted_basis(chars_ext) +function symmetry_adapted_basis( + G::Group, + action::Action, + basis, + S::Type = Rational{Int}; + semisimple=false, +) + tbl = CharacterTable(S, G) + ehom = CachedExtensionHomomorphism(parent(tbl), action, basis, precompute=true) + return symmetry_adapted_basis(eltype(tbl), tbl, ehom, semisimple=semisimple) end function symmetry_adapted_basis( - ::Type{T}, + T::Type, G::Group, + action::Action, basis, - action, - ::Type{S} = Rational{Int}, -) where {T,S} - chars_ext = extended_characters(S, G, basis, action) - return symmetry_adapted_basis(T, chars_ext) + S::Type = Rational{Int}; + semisimple=false, +) + tbl = CharacterTable(S, G) + ehom = CachedExtensionHomomorphism(parent(tbl), action, basis, precompute=true) + return symmetry_adapted_basis(T, tbl, ehom, semisimple=semisimple) end -symmetry_adapted_basis(chars::AbstractVector{<:AbstractClassFunction{T}}) where T = - symmetry_adapted_basis(T, chars) +function symmetry_adapted_basis( + T::Type, + tbl::CharacterTable, + ehom::InducedActionHomomorphism; + semisimple=false, +) + ψ = action_character(ehom, tbl) + + irr, multips = _constituents_decomposition(ψ, tbl) + if T <: Real + irr, multips = affordable_real(irr, multips) + end -symmetry_adapted_basis(::Type{T}, chars::AbstractVector{<:AbstractClassFunction}) where T = - _symmetry_adapted_basis(Character{T}.(chars)) + if semisimple || all(isone ∘ degree, irr) + return _symmetry_adapted_basis(T, irr, multips, ehom) + else + RG = _group_algebra(parent(tbl)) + return _symmetry_adapted_basis(T, irr, multips, RG, ehom) + end +end -symmetry_adapted_basis(T::Type{<:Real}, chars::AbstractVector{<:AbstractClassFunction}) = - _symmetry_adapted_basis(Character{T}.(affordable_real(chars))) +function _constituents_decomposition(ψ::Character, tbl::CharacterTable) + irr = irreducible_characters(tbl) + degrees = degree.(irr) + multiplicities = constituents(ψ) -symmetry_adapted_basis(::Type{T}, chars::AbstractVector{<:AbstractClassFunction{T}}) where T = - _symmetry_adapted_basis(chars) + @debug "Decomposition into character spaces: + degrees: $(join([lpad(d, 6) for d in degrees], "")) + multiplicities: $(join([lpad(m, 6) for m in multiplicities], ""))" + @assert dot(multiplicities, degrees) == degree(ψ) + "Something went wrong: characters do not constitute a complete basis for action: + $(dot(multiplicities, degrees)) ≠ $(degree(ψ))" -_multiplicities(ψ, chars) = Int[div(Int(dot(ψ, χ)), Int(dot(χ, χ))) for χ in chars] -_multiplicities(ψ, chars::AbstractVector{<:AbstractClassFunction{T}}) where T<:AbstractFloat = - [round(Int, dot(ψ, χ) / dot(χ, χ)) for χ in chars] -_multiplicities(ψ, chars::AbstractVector{<:AbstractClassFunction{T}}) where T<:Complex = - [round(Int, real(dot(ψ, χ) / dot(χ, χ))) for χ in chars] + present_irreps = [i for (i, m) in enumerate(multiplicities) if m ≠ 0] + return irr[present_irreps], multiplicities[present_irreps] +end -macro spawn_compat(expr) - @static if VERSION < v"1.3.0" - return :(@async $(esc(expr))) +function _group_algebra(G::Group) + @assert isfinite(G) + b = StarAlgebras.Basis{UInt16}(vec(collect(G))) + RG = if order(Int, G) <= (typemax(UInt16)>>2) + StarAlgebra(G, b, (length(b), length(b)), precompute=true) + # cache is about ~ 1Gb else - return :(Threads.@spawn $(esc(expr))) + StarAlgebra(G, b) end + return RG end -function _symmetry_adapted_basis(chars::AbstractVector{<:AbstractClassFunction}) - ψ = let χ = first(chars) - action_character(conjugacy_classes(χ), χ.inv_of) +function _symmetry_adapted_basis( + T::Type, + irr::AbstractVector{<:Character}, + multiplicities::AbstractVector{<:Integer}, + hom=nothing +) + res = map(zip(irr, multiplicities)) do (µ, m) + @spawn_compat begin + µT = eltype(µ) == T ? µ : Character{T}(µ) + image = hom === nothing ? image_basis(µT) : image_basis(hom, µT) + simple = size(image, 1) == m + deg = degree(µ) + if deg == 1 + @assert simple "Central projection associated to character is not simple unless its degree == 1" + end + DirectSummand(image, m, deg, simple) + end end + return fetch.(res) +end - multiplicities = _multiplicities(ψ, chars) - degrees = degree.(chars) - - @debug "Decomposition into character spaces: - degrees: $(join([lpad(d, 6) for d in degrees], "")) - multiplicities: $(join([lpad(m, 6) for m in multiplicities], ""))" +function _symmetry_adapted_basis( + T::Type, + irr::AbstractVector{<:Character}, + multiplicities::AbstractVector{<:Integer}, + RG::StarAlgebra{<:Group}, + hom=nothing, +) + mps, simples = minimal_projection_system(irr, RG) + degrees = degree.(irr) + res = map(zip(mps, multiplicities, degrees, simples)) do (µ, m, deg, simple) + @spawn_compat begin + µT = eltype(µ) == T ? µ : AlgebraElement{T}(µ) + image = hom === nothing ? image_basis(µT) : image_basis(hom, µT) + DirectSummand(image, m, deg, simple) + end + end + direct_summands = fetch.(res) - dot(multiplicities, degrees) == degree(ψ) || - @error "characters do not constitute a complete basis for action: - $(dot(multiplicities, degrees)) ≠ $(degree(ψ))" + for (χ, ds) in zip(irr, direct_summands) + if issimple(ds) && (d = size(ds, 1)) != (e = multiplicity(ds)*sum(constituents(χ).>0)) + throw("The dimension of the projection doesn't match with simple summand multiplicity: $d ≠ $e") + end + end - res = [@spawn_compat(SemisimpleSummand(isotypical_basis(χ), m)) for (χ, m) - in zip(chars, multiplicities) if m ≠ 0] - - return fetch.(res) + return direct_summands end diff --git a/src/wedderburn_decomposition.jl b/src/wedderburn_decomposition.jl new file mode 100644 index 0000000..63d997c --- /dev/null +++ b/src/wedderburn_decomposition.jl @@ -0,0 +1,167 @@ +struct WedderburnDecomposition{B, iV, DS<:DirectSummand, Hom} + basis::B + invariants::iV + Uπs::Vector{DS} + hom::Hom +end + +function WedderburnDecomposition( + T::Type, + G::Group, + action::Action, + basis_full, + basis_half, + S = Rational{Int}; +) + basis = StarAlgebras.Basis{UInt32}(basis_full) + invariants = invariant_vectors(G, action, basis) + + tbl = CharacterTable(S, G) + ehom = CachedExtensionHomomorphism(G, action, basis_half, precompute=true) + + Uπs = let + sa_basis = symmetry_adapted_basis(T, tbl, ehom; semisimple=false) + @debug "fill factor of sa_basis:" round.(_fillfactor.(sa_basis), digits=3) + sp_basis = sparse.(sa_basis) + if T <: AbstractFloat + droptol!.(sp_basis, eps(T)*length(basis_half)) + @debug "fill factor after sparsification" round.(_fillfactor.(sa_basis), digits=3) + end + sp_basis + end + + return WedderburnDecomposition(basis, invariants, Uπs, ehom) +end + +invariant_vectors(wbdec::WedderburnDecomposition) = wbdec.invariants +StarAlgebras.basis(wbdec::WedderburnDecomposition) = wbdec.basis +direct_summands(wbdec::WedderburnDecomposition) = wbdec.Uπs + +_tmps(wbdec::WedderburnDecomposition) = + [zeros(eltype(U), reverse(size(U))) for U in wbdec.Uπs] + +_fillfactor(M::AbstractMatrix) = count(!iszero, M)/length(M) +_fillfactor(M::AbstractSparseMatrix) = nnz(M)/length(M) + +function diagonalize( + M::AbstractMatrix, + wbdec::WedderburnDecomposition, + tmps=_tmps(wbdec) +) + # return [degree(Uπ).*(Uπ*(M*Uπ')) for Uπ in summands(wbdec)] + + T = eltype(eltype(direct_summands(wbdec))) + Mπs = [(d = size(Uπ, 1); zeros(T, d,d)) for Uπ in direct_summands(wbdec)] + return diagonalize!(Mπs, M, direct_summands(wbdec), tmps) +end + +function diagonalize!( + Mπs, + M::AbstractMatrix, + wbdec::WedderburnDecomposition, + tmps=_tmps(wbdec)) + return diagonalize!(Mπs, M, direct_summands(wbdec), tmps) +end + +function diagonalize!( + Mπs::AbstractVector{<:AbstractMatrix}, + M::AbstractMatrix, + Uπs::AbstractVector{<:DirectSummand}, + tmps +) + for (π, Uπ) in enumerate(Uπs) + imUπ = image_basis(Uπ) # Base.Matrix to allow BLAS paths below + LinearAlgebra.mul!(tmps[π], M, imUπ') + LinearAlgebra.mul!(Mπs[π], imUπ, tmps[π]) + deg = degree(Uπ) + Mπs[π] .*= deg + end + return Mπs +end + + +function invariant_vectors( + G::Group, + act::ByPermutations, + basis::StarAlgebras.Basis{T, I}, +) where {T, I} + + tovisit = trues(length(basis)); + invariant_vs = Vector{SparseVector{Rational{Int}}}() + + ordG = order(Int, G) + elts = collect(G) + orbit = zeros(I, ordG) + + for i in eachindex(basis) + if tovisit[i] + bi = basis[i] + Threads.@threads for j in eachindex(elts) + orbit[j] = basis[action(act, elts[j], bi)] + end + tovisit[orbit] .= false + push!(invariant_vs, sparsevec(orbit, fill(1//ordG, ordG), length(basis))) + end + end + return invariant_vs +end + +_orth_proj(A::AbstractMatrix, v, QR=qr(A)) = A * (QR\v) + +function invariant_vectors( + G::Group, + act::ByLinearTransformation, + basis, + atol=1e-12 + ) + + hom = CachedExtensionHomomorphism(G, act, basis, precompute = true) + + dim = 0 + T = coeff_type(action(hom)) + @assert T==Float64 + subspaces = SparseMatrixCSC{T, Int64}[] + qrs = LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}}[] + + inv_vectors = SparseVector{T, Int64}[] + inv_v = Vector{T}(undef, length(basis)) + + for b in basis + b_orth = sparsevec(decompose(b, hom)..., length(basis)) + if dim > 0 + b_orth -= sum(_orth_proj(V, b_orth, QR) for (V, QR) in zip(subspaces, qrs)) + @debug "orthogonalizing $b:" norm(b_orth) + norm(b_orth, 2) < atol && continue + end + + # we found a new subspace! + b_orth ./= norm(b_orth, 2) + orbit_subspace = Matrix{T}(undef, length(basis), 1) + orbit_subspace[:, 1] .= b_orth + inv_v .= zero(T) + + # for the orbit of b_orth find + # * `inv_v`, the invariant (average) vector, and + # * `tmp_subspace`, the subspace spanned by the orbit + + for g in G + k = induce(hom, g) * b_orth + @assert isapprox(norm(k, 2), 1.0, atol=atol) + inv_v .+= k + + residual = k - _orth_proj(orbit_subspace, k) + (rnorm = norm(residual, 2)) < atol && continue + orbit_subspace = hcat(orbit_subspace, Vector(residual./rnorm)) + end + + push!(inv_vectors, inv_v/order(Int, G)) # a copy if inv_v + push!(subspaces, orbit_subspace) # orbit_subspace is freshly allocated every nonzero b_orth + dim += size(orbit_subspace, 2) + dim == length(basis) && break + dim > length(basis) && throw("This should never occur") + # moving qr here to save one qr decomposition + push!(qrs, qr(orbit_subspace)) + end + @assert dim == length(basis) "Encountered incomplete decomposition into orbits!" + return inv_vectors +end diff --git a/test/action_dihedral.jl b/test/action_dihedral.jl new file mode 100644 index 0000000..8bf7244 --- /dev/null +++ b/test/action_dihedral.jl @@ -0,0 +1,92 @@ +using LinearAlgebra +using SparseArrays + +using SymbolicWedderburn +using DynamicPolynomials + +include(joinpath(dirname(@__DIR__), "examples", "action_polynomials.jl")) +include(joinpath(dirname(@__DIR__), "examples", "dihedral.jl")) +include(joinpath(dirname(@__DIR__), "examples", "util.jl")) + +using JuMP +using SCS + +const OPTIMIZER = optimizer_with_attributes( + SCS.Optimizer, + "acceleration_lookback" => 10, + "max_iters" => 3_000, + "alpha" => 1.5, + "eps" => 1e-6, + "linear_solver" => SCS.DirectSolver, +) + +@polyvar x y +const robinson_form = x^6 + y^6 - x^4 * y^2 - y^4 * x^2 - x^4 - y^4 - x^2 - y^2 + 3x^2 * y^2 + 1 + +struct DihedralAction <: OnMonomials end + +SymbolicWedderburn.coeff_type(::DihedralAction) = Float64 +function SymbolicWedderburn.action(::DihedralAction, el::DihedralElement, mono::AbstractMonomial) + if iseven(el.reflection + el.id) + var_x, var_y = x, y + else + var_x, var_y = y, x + end + sign_x = 1 <= el.id <= 2 ? -1 : 1 + sign_y = 2 <= el.id ? -1 : 1 + return mono([x, y] => [sign_x * var_x, sign_y * var_y]) +end + +@testset "Dihedral Action" begin + G = DihedralGroup(4) + @test all( + SymbolicWedderburn.action(DihedralAction(), g, robinson_form) == robinson_form + for g in DihedralGroup(4) + ) + + T = Float64 + + wedderburn = WedderburnDecomposition( + T, + G, + DihedralAction(), + monomials([x,y], 0:DynamicPolynomials.maxdegree(robinson_form)), + monomials([x,y], 0:(DynamicPolynomials.maxdegree(robinson_form) ÷ 2)) + ) + + M = let bh = monomials([x,y], 0:(DynamicPolynomials.maxdegree(robinson_form) ÷ 2)) + [SymbolicWedderburn.basis(wedderburn)[x * y] for x in bh, y in bh] + end + + m = let m = JuMP.Model(OPTIMIZER) + JuMP.@variable m t + JuMP.@objective m Max t + psds = [ + JuMP.@variable(m, [1:d, 1:d] in PSDCone()) + for d in size.(direct_summands(wedderburn), 1) + ] + + # preallocating + Mπs = zeros.(T, size.(psds)) + M_orb = similar(M, T) + tmps = zeros.(T, reverse.(size.(direct_summands(wedderburn)))) + + C = DynamicPolynomials.coefficients(robinson_form - t, 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, tmps) + + JuMP.@constraint m sum( + dot(Mπ, Pπ) for (Mπ, Pπ) in zip(Mπs, psds) if !iszero(Mπ) + ) == c + end + m + end + + optimize!(m) + + @test isapprox(value(m[:t]), -3825 / 4096, rtol = 1e-4) + @test termination_status(m) == MOI.OPTIMAL +end diff --git a/test/action_linear.jl b/test/action_linear.jl index 70f58f2..7cce826 100644 --- a/test/action_linear.jl +++ b/test/action_linear.jl @@ -4,15 +4,14 @@ using GroupsCore include(joinpath(pathof(GroupsCore), "..", "..", "test", "cyclic.jl")) x,y = @polyvar x y -struct MyAction <: SymbolicWedderburn.ByLinearTransformation end -# (x, y) → (x+y, x-y) +struct By90Rotation <: SymbolicWedderburn.ByLinearTransformation end +# (x, y) → (x+y, x-y)/sqrt(2) -SymbolicWedderburn.coeff_type(m::MyAction) = Float64 -function SymbolicWedderburn.action(a::MyAction, g::CyclicGroupElement, m::Monomial) +SymbolicWedderburn.coeff_type(::By90Rotation) = Float64 +function SymbolicWedderburn.action(::By90Rotation, g::CyclicGroupElement, m::Monomial) isone(g) && return m x, y = variables(m) - # we need the action to be orthogonal - return m(x => (x+y)/sqrt(2), y => (x-y)/sqrt(2)) + return m(x => (x+y)/sqrt(2), y => (x-y)/sqrt(2)) # we need the action to be orthogonal end function SymbolicWedderburn.decompose(k::AbstractPolynomial, hom::SymbolicWedderburn.InducedActionHomomorphism) @@ -27,34 +26,34 @@ end @testset "Linear Actions" begin G = CyclicGroup(2) - basis = monomials([x,y], 0:4) - ehom = SymbolicWedderburn.ExtensionHomomorphism(MyAction(), basis) + monomial_basis = monomials([x,y], 0:4) + ehom = SymbolicWedderburn.ExtensionHomomorphism(By90Rotation(), monomial_basis) @testset "decomposition in basis" begin g = gens(G, 1) - k = SymbolicWedderburn.action(MyAction(), g, basis[2]) + k = SymbolicWedderburn.action(By90Rotation(), g, monomial_basis[2]) idcs, vals = SymbolicWedderburn.decompose(k, ehom) - @test sum(c*basis[i] for (c,i) in zip(vals, idcs)) == k + @test sum(c*monomial_basis[i] for (c,i) in zip(vals, idcs)) == k - for b in basis - k = SymbolicWedderburn.action(MyAction(), g, b) + for b in monomial_basis + k = SymbolicWedderburn.action(By90Rotation(), g, b) idcs, vals = SymbolicWedderburn.decompose(k, ehom) - @test sum(c*basis[i] for (c,i) in zip(vals, idcs)) == k + @test sum(c*monomial_basis[i] for (c,i) in zip(vals, idcs)) == k end end @testset "induced Matrix Representation" begin g = gens(G, 1) - M = SymbolicWedderburn.induce(MyAction(), ehom, g) + M = SymbolicWedderburn.induce(By90Rotation(), ehom, g) @test M isa AbstractMatrix{Float64} - @test size(M) == (length(basis), length(basis)) + @test size(M) == (length(monomial_basis), length(monomial_basis)) @test !(M ≈ one(M)) @test M^2 ≈ one(M) - B = SymbolicWedderburn.symmetry_adapted_basis(G, basis, MyAction()); - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, B) == length(basis) + sa_basis = symmetry_adapted_basis(Float64, G, By90Rotation(), monomial_basis); + @test sum(first ∘ size, sa_basis) == length(monomial_basis) end end diff --git a/test/action_permutation.jl b/test/action_permutation.jl index b162a65..cc1495f 100644 --- a/test/action_permutation.jl +++ b/test/action_permutation.jl @@ -13,18 +13,16 @@ Base.:(==)(w::Word, v::Word) = w.alphabet == v.alphabet && w.letters == v.letter Base.hash(w::Word, h::UInt=UInt(0)) = hash(w.alphabet, hash(w.letters, hash(Word, h))) struct OnLetters <: SymbolicWedderburn.ByPermutations end -function SymbolicWedderburn.action(a::OnLetters, p::PermutationGroups.AbstractPerm, w::Word) +function SymbolicWedderburn.action(::OnLetters, p::PermutationGroups.AbstractPerm, w::Word) return Word(w.alphabet, [l^p for l in w.letters]) end -my_action(w::Word, g) = SymbolicWedderburn.action(OnLetters(), g, w) - @testset "Extending homomorphism" begin words = let A = [:a, :b, :c] w = Word(A, [1,2,3,2,1]) # (a·b·c·b·a)^(2,3) == a·c·b·c·a - @test my_action(w, perm"(2,3)") == Word(A, [1,3,2,3,1]) + @test SymbolicWedderburn.action(OnLetters(), perm"(2,3)", w) == Word(A, [1,3,2,3,1]) words = [Word(A, [1]), Word(A, [2]), Word(A, [3])] for r in 2:4 @@ -37,10 +35,119 @@ my_action(w::Word, g) = SymbolicWedderburn.action(OnLetters(), g, w) end G = PermGroup(perm"(1,2,3)",perm"(1,2)") # G acts on words permuting letters - sa_basis = symmetry_adapted_basis(G, words, OnLetters()) + basis = words + action = OnLetters() + tbl = SymbolicWedderburn.CharacterTable(Rational{Int}, G) + ehom = SymbolicWedderburn.CachedExtensionHomomorphism(G, action, basis, precompute=true) + @test all(g ∈ keys(ehom.cache) for g in G) # we actually cached + + ψ = SymbolicWedderburn.action_character(ehom, tbl) + @test SymbolicWedderburn.constituents(ψ) == [40, 22, 18] + irr = SymbolicWedderburn.irreducible_characters(tbl) + multips = SymbolicWedderburn.constituents(ψ) + @test dot(SymbolicWedderburn.degree.(irr), multips) == length(basis) + simple = isone.(SymbolicWedderburn.degree.(irr)) + @test simple == [false, true, true] + + @testset "semisimple decomposition" begin + let i = 1 + χ, m, s = irr[i], multips[i], simple[i] + b = SymbolicWedderburn.image_basis(ehom, χ) + @test size(b, 1) == SymbolicWedderburn.degree(χ)*m == 80 + end + + let i = 2 + χ, m, s = irr[i], multips[i], simple[i] + b = SymbolicWedderburn.image_basis(ehom, χ) + @test size(b, 1) == SymbolicWedderburn.degree(χ)*m == 22 + end + + let i = 3 + χ, m, s = irr[i], multips[i], simple[i] + b = SymbolicWedderburn.image_basis(ehom, χ) + @test size(b, 1) == SymbolicWedderburn.degree(χ)*m == 18 + end + + @test symmetry_adapted_basis(G, action, basis, semisimple=true) isa + Vector{<:SymbolicWedderburn.DirectSummand{<:SymbolicWedderburn.Cyclotomic}} + @test symmetry_adapted_basis(Rational{Int}, G, action, basis, semisimple=true) isa + Vector{<:SymbolicWedderburn.DirectSummand{Rational{Int}}} + @test symmetry_adapted_basis(Float64, G, action, basis, semisimple=true) isa + Vector{<:SymbolicWedderburn.DirectSummand{Float64}} + + sa_basis = symmetry_adapted_basis(G, action, basis, semisimple=true) + + @test [convert(Matrix{Rational{Int}}, b) for b in sa_basis] isa Vector{Matrix{Rational{Int}}} + @test [convert(Matrix{Float64}, b) for b in sa_basis] isa Vector{Matrix{Float64}} + + @test length(sa_basis) == 3 + @test multiplicity.(sa_basis) == [40, 22, 18] + @test SymbolicWedderburn.degree.(sa_basis) == [2,1,1] + @test size.(sa_basis, 1) == + multips .* SymbolicWedderburn.degree.(irr) == + [80, 22, 18] + @test sum(first ∘ size, sa_basis) == length(words) + end + + @testset "simple decomposition" begin + RG = let G = G + b = StarAlgebras.Basis{UInt16}(collect(G)) + StarAlgebra(G, b, (length(b), length(b))) + end - @test sa_basis isa Vector{<:SymbolicWedderburn.SemisimpleSummand{<:Matrix{<:SymbolicWedderburn.Cyclotomic}}} - @test [convert(Matrix{Int}, b) for b in sa_basis] isa Vector{Matrix{Int}} - @test [convert(Matrix{Float64}, b) for b in sa_basis] isa Vector{Matrix{Float64}} - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, sa_basis) == length(words) + let χ = irr[1], m = multips[1] + (a, fl) = SymbolicWedderburn.minimal_rank_projection(χ, RG) + @test fl == isone(χ(a)) + + µ = AlgebraElement(χ, RG)*a + mpr = SymbolicWedderburn.image_basis(ehom, µ) + @test mpr isa Matrix{eltype(µ)} + @test size(mpr, 1) == m + + µR = AlgebraElement{Rational{Int}}(µ) + mpr = SymbolicWedderburn.image_basis(ehom, µR) + @test mpr isa Matrix{eltype(µR)} + @test size(mpr, 1) == m + + µFl = AlgebraElement{Float64}(µ) + mpr = SymbolicWedderburn.image_basis(ehom, µFl) + @test mpr isa Matrix{eltype(µFl)} + @test size(mpr, 1) == m + end + + let χ = irr[2], m = multips[2] + (a, fl) = SymbolicWedderburn.minimal_rank_projection(χ, RG) + @test fl == isone(χ(a)) + + µ = AlgebraElement(χ, RG)*a + mpr = SymbolicWedderburn.image_basis(ehom, µ) + @test mpr isa Matrix{eltype(µ)} + @test size(mpr, 1) == m + end + + let χ = irr[3], m = multips[3] + (a, fl) = SymbolicWedderburn.minimal_rank_projection(χ, RG) + @test fl == isone(χ(a)) + + µ = AlgebraElement(χ, RG)*a + mpr = SymbolicWedderburn.image_basis(ehom, µ) + @test mpr isa Matrix{eltype(µ)} + @test size(mpr, 1) == m + end + + @test symmetry_adapted_basis(G, action, basis, semisimple=false) isa + Vector{<:SymbolicWedderburn.DirectSummand{<:SymbolicWedderburn.Cyclotomic}} + @test symmetry_adapted_basis(Rational{Int}, G, action, basis, semisimple=false) isa + Vector{<:SymbolicWedderburn.DirectSummand{Rational{Int}}} + @test symmetry_adapted_basis(Float64, G, action, basis, semisimple=false) isa + Vector{<:SymbolicWedderburn.DirectSummand{Float64}} + + sa_basis = symmetry_adapted_basis(G, action, basis) + + @test length(sa_basis) == 3 + @test multiplicity.(sa_basis) == [40, 22, 18] + @test SymbolicWedderburn.degree.(sa_basis) == [2,1,1] + @test all(issimple, sa_basis) + @test size.(sa_basis, 1) == multips == [40, 22, 18] + end end diff --git a/test/ccmatrix.jl b/test/ccmatrix.jl index d14b0f5..69c0592 100644 --- a/test/ccmatrix.jl +++ b/test/ccmatrix.jl @@ -1,11 +1,11 @@ function generic_tests_ccmatrix(C) r = length(C) for i = 1:r - ccm = SymbolicWedderburn.CMMatrix(C, i) + ccm = Characters.CMMatrix(C, i) l = length(C[i]) @test all([sum(ccm[:,j]) == l for j = 1:r]) end - ccm1 = SymbolicWedderburn.CMMatrix(C, 1) + ccm1 = Characters.CMMatrix(C, 1) @test all([isone(ccm1[i,i]) for i = 1:r]) end @@ -13,30 +13,30 @@ end @testset "example: Symmetric group (4)" begin G = PermutationGroups.SymmetricGroup(4) - @test conjugacy_classes(G) isa AbstractVector{<:AbstractOrbit} + @test Characters.conjugacy_classes(G) isa AbstractVector{<:AbstractOrbit} - for cc in conjugacy_classes(G) + for cc in Characters.conjugacy_classes(G) @test all(permtype(g) == permtype(first(cc)) for g in cc) end - C = conjugacy_classes(G) + C = Characters.conjugacy_classes(G) generic_tests_ccmatrix(C) - @test SymbolicWedderburn.CMMatrix(C, 2) == [0 1 0 0 0; + @test Characters.CMMatrix(C, 2) == [0 1 0 0 0; 6 0 3 0 2; 0 4 0 4 0; 0 0 3 0 4; 0 1 0 2 0] - @test SymbolicWedderburn.CMMatrix(C, 3) == [0 0 1 0 0; + @test Characters.CMMatrix(C, 3) == [0 0 1 0 0; 0 4 0 4 0; 8 0 4 0 8; 0 4 0 4 0; 0 0 3 0 0] - @test SymbolicWedderburn.CMMatrix(C, 4) == [0 0 0 1 0; + @test Characters.CMMatrix(C, 4) == [0 0 0 1 0; 0 0 3 0 4; 0 4 0 4 0; 6 0 3 0 2; 0 2 0 1 0] - @test SymbolicWedderburn.CMMatrix(C, 5) == [0 0 0 0 1; + @test Characters.CMMatrix(C, 5) == [0 0 0 0 1; 0 1 0 2 0; 0 0 3 0 0; 0 2 0 1 0; @@ -54,19 +54,19 @@ end ] generic_tests_ccmatrix(C) - @test SymbolicWedderburn.CMMatrix(C, 1) == [1 0 0 0 + @test Characters.CMMatrix(C, 1) == [1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1] - @test SymbolicWedderburn.CMMatrix(C, 2) == [0 1 0 0 + @test Characters.CMMatrix(C, 2) == [0 1 0 0 3 2 0 0 0 0 3 0 0 0 0 3] - @test SymbolicWedderburn.CMMatrix(C, 3) == [0 0 1 0 + @test Characters.CMMatrix(C, 3) == [0 0 1 0 0 0 3 0 0 0 0 4 4 4 0 0] - @test SymbolicWedderburn.CMMatrix(C, 4) == [0 0 0 1 + @test Characters.CMMatrix(C, 4) == [0 0 0 1 0 0 0 3 4 4 0 0 0 0 4 0] @@ -87,31 +87,31 @@ end @assert sum(length, ccG) == order(G) - @test SymbolicWedderburn.CMMatrix(ccG, 1) == + @test Characters.CMMatrix(ccG, 1) == [ 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 ] - @test SymbolicWedderburn.CMMatrix(ccG, 2) == + @test Characters.CMMatrix(ccG, 2) == [ 0 1 0 0 0 5 0 0 0 5 0 0 0 5 0 0 0 5 0 0 0 4 0 0 0 ] - @test SymbolicWedderburn.CMMatrix(ccG, 3) == + @test Characters.CMMatrix(ccG, 3) == [ 0 0 1 0 0 0 0 0 5 0 0 5 0 0 0 5 0 0 0 5 0 0 4 0 0 ] - @test SymbolicWedderburn.CMMatrix(ccG, 4) == + @test Characters.CMMatrix(ccG, 4) == [ 0 0 0 1 0 0 0 5 0 0 5 0 0 0 5 0 5 0 0 0 0 0 0 4 0 ] - @test SymbolicWedderburn.CMMatrix(ccG, 5) == + @test Characters.CMMatrix(ccG, 5) == [ 0 0 0 0 1 0 4 0 0 0 0 0 4 0 0 @@ -119,7 +119,7 @@ end 4 0 0 0 3 ] generic_tests_ccmatrix(ccG) - generic_tests_ccmatrix(conjugacy_classes(G)) # might be in different order + generic_tests_ccmatrix(Characters.conjugacy_classes(G)) # might be in different order end @testset "random subgroups of SymetricGroup(N)" begin @@ -127,7 +127,7 @@ end G = PermutationGroups.SymmetricGroup(i) for _ in 1:5 PG = PermGroup(rand(G, 2)) - generic_tests_ccmatrix(SymbolicWedderburn.conjugacy_classes(PG)) + generic_tests_ccmatrix(Characters.conjugacy_classes(PG)) end end end diff --git a/test/characters.jl b/test/characters.jl index d2803db..d6b8eee 100644 --- a/test/characters.jl +++ b/test/characters.jl @@ -1,41 +1,58 @@ @testset "Frobenius-Schur & projections" begin G = SmallPermGroups[8][4] - ccG = conjugacy_classes(G) + irr = Characters.irreducible_characters(G) # quaternionic character - χ = SymbolicWedderburn.Character(Float64[2, 0, 0, -2, 0], ccG) + χ = irr[1] + @test collect(values(χ)) == [2, 0, 0, -2, 0] - ι = SymbolicWedderburn.frobenius_schur_indicator + ι = Characters.frobenius_schur @test ι(χ) == -1 @test deepcopy(χ) == χ @test deepcopy(χ) !== χ - ψ = SymbolicWedderburn.Character(copy(values(χ)), conjugacy_classes(χ)) + ψ = deepcopy(χ) @test ψ == χ @test ψ !== χ @test ψ == deepcopy(χ) - @test ψ.cc === χ.cc - @test conjugacy_classes(ψ) === conjugacy_classes(χ) + @test Characters.conjugacy_classes(ψ) === Characters.conjugacy_classes(χ) + @test Characters.table(ψ) === Characters.table(χ) - @test values(ψ) == values(χ) - @test values(ψ) !== values(χ) + @test collect(values(ψ)) == collect(values(χ)) + @test Characters.constituents(ψ) == Characters.constituents(χ) + @test Characters.constituents(ψ) !== Characters.constituents(χ) @test hash(ψ) == hash(χ) @test 2ψ == 2χ @test 2ψ/2 == χ @test 2ψ != χ - @test hash((2ψ/2).vals) == hash(χ.vals) @test PermutationGroups.degree(χ) == 2 - @test size(SymbolicWedderburn.isotypical_basis(χ), 1) == 4 + # @test size(SymbolicWedderburn.image_basis(χ), 1) == 4 - @test SymbolicWedderburn.affordable_real!(deepcopy(χ)) isa SymbolicWedderburn.Character + @test_throws AssertionError ι(2χ) + a = 2χ + @test sum(length(cc)*a(first(cc)^2) for cc in Characters.conjugacy_classes(a))//order(G) == 2ι(χ) +end + +@testset "Characters io" begin + G = PermGroup([perm"(2,3)(4,5)"]) + chars = Characters.irreducible_characters(G) + + @test sprint(show, chars[1]) == "χ₁" + @test sprint(show, chars[2]) == "χ₂" + @test sprint(show, chars[1]+2chars[2]) == "χ₁ +2·χ₂" + @test sprint(show, 2chars[1]-chars[2]) == "2·χ₁ -1·χ₂" + @test sprint(show, -2chars[1]+3chars[2]) == "-2·χ₁ +3·χ₂" - @test ι(SymbolicWedderburn.affordable_real!(deepcopy(χ))) == 2*ι(χ) + if VERSION > v"1.3.0" + @test sprint(show, MIME"text/plain"(), chars[1]) == + "Character over Cyclotomic{Rational{Int64}, SparseArrays.SparseVector{Rational{Int64}, Int64}}\nχ₁" + end end diff --git a/test/dixon.jl b/test/dixon.jl index cdf44a5..0f8d5b8 100644 --- a/test/dixon.jl +++ b/test/dixon.jl @@ -1,9 +1,9 @@ -function generictest_dixon_Fp(G, p = SymbolicWedderburn.dixon_prime(G)) - F = SymbolicWedderburn.FiniteFields.GF{p} - ccG = conjugacy_classes(G) - Ns = [SymbolicWedderburn.CMMatrix(ccG, i) for i = 1:length(ccG)] - @test isdiag(SymbolicWedderburn.common_esd(Ns, F)) - esd = SymbolicWedderburn.common_esd(Ns, F) +function generictest_dixon_Fp(G, p = Characters.dixon_prime(G)) + F = Characters.FiniteFields.GF{p} + ccG = Characters.conjugacy_classes(G) + Ns = [Characters.CMMatrix(ccG, i) for i = 1:length(ccG)] + @test isdiag(Characters.common_esd(Ns, F)) + esd = Characters.common_esd(Ns, F) let W = esd.basis @test inv(W) isa Matrix{F} @@ -16,53 +16,56 @@ function generictest_dixon_Fp(G, p = SymbolicWedderburn.dixon_prime(G)) end end - @test SymbolicWedderburn.characters_dixon(F, ccG) isa - Vector{<:SymbolicWedderburn.Character} + tbl = Characters.CharacterTable(F, G, ccG) + @test Characters.irreducible_characters(tbl) isa + Vector{<:Characters.Character} - chars = SymbolicWedderburn.characters_dixon(F, ccG) + chars = Characters.irreducible_characters(tbl) # checking the degrees - @test sum(Int.(degree.(chars)) .^ 2) == order(G) + @test sum(Int.(SymbolicWedderburn.degree.(chars)) .^ 2) == order(G) # orthogonality of characters over F @test [dot(χ, ψ) for χ in chars, ψ in chars] == I end -function generictest_dixon_C(G, p = SymbolicWedderburn.dixon_prime(G)) - F = SymbolicWedderburn.FiniteFields.GF{p} - ccG = conjugacy_classes(G) - chars_Fp = SymbolicWedderburn.characters_dixon(F, ccG) +function generictest_dixon_C(G, p = Characters.dixon_prime(G)) + F = Characters.FiniteFields.GF{p} + ccG = Characters.conjugacy_classes(G) + tbl = Characters.CharacterTable(F, G, ccG) + chars_Fp = Characters.irreducible_characters(tbl) - degrees = degree.(chars_Fp) + degrees = SymbolicWedderburn.degree.(chars_Fp) - m = SymbolicWedderburn._multiplicities(chars_Fp) + m = Characters._multiplicities(chars_Fp) for i = 1:size(m, 1) @test all(Int.(m[i, :, :]) .<= degrees[i]) end - @test SymbolicWedderburn.complex_characters(Int, chars_Fp) isa - Vector{<:SymbolicWedderburn.Character} - @test length(SymbolicWedderburn.complex_characters(Int, chars_Fp[1:1])) == 1 + @test Characters.complex_character_table(Rational{Int}, tbl) isa + Characters.CharacterTable{<:Group, <:Cyclotomic{Rational{Int}}} - chars_CC = SymbolicWedderburn.complex_characters(Int, chars_Fp) - id = one(first(first(ccG))) + tblCC = Characters.complex_character_table(Rational{Int}, tbl) + chars_CC = Characters.irreducible_characters(tblCC) + + id = one(G) @test [ψ(id) for ψ in chars_CC] == degrees # orthogonality of characters over ℂ @test [dot(χ, ψ) for χ in chars_CC, ψ in chars_CC] == I - @test sum(χ -> degree(χ)^2, chars_CC) == order(G) + @test sum(χ -> SymbolicWedderburn.degree(χ)^2, chars_CC) == order(G) end @testset "Dixon Algorithm" begin @testset "DixonPrimes" begin - @test SymbolicWedderburn.dixon_prime(20, 20) == 41 + @test Characters.dixon_prime(20, 20) == 41 @testset "random" begin for n in rand(2:1000, 10) - F = SymbolicWedderburn.Primes.factor(n) + F = Characters.Primes.factor(n) e = lcm(collect(keys(F))) - p = SymbolicWedderburn.dixon_prime(n, e) + p = Characters.dixon_prime(n, e) @test mod(p, e) == 1 @test p > 2 * floor(sqrt(n)) end @@ -70,11 +73,11 @@ end @testset "DixonPrimeGroups" begin G = PermutationGroups.SymmetricGroup(4) - ccG = conjugacy_classes(G) + ccG = Characters.conjugacy_classes(G) @test exponent(G) == 12 @test exponent(ccG) == 12 - @test SymbolicWedderburn.dixon_prime(G) == - SymbolicWedderburn.dixon_prime(ccG) + @test Characters.dixon_prime(G) == + Characters.dixon_prime(ccG) end end @@ -91,34 +94,35 @@ end end let ccG = ccG, p = 29 - F = SymbolicWedderburn.FiniteFields.GF{p} - Ns = [SymbolicWedderburn.CMMatrix(ccG, i) for i = 1:length(ccG)] - @test_throws AssertionError SymbolicWedderburn.common_esd(Ns, F) + F = Characters.FiniteFields.GF{p} + Ns = [Characters.CMMatrix(ccG, i) for i = 1:length(ccG)] + @test_throws AssertionError Characters.common_esd(Ns, F) end let ccG = ccG, p = 31 - F = SymbolicWedderburn.FiniteFields.GF{p} - Ns = [SymbolicWedderburn.CMMatrix(ccG, i) for i = 1:length(ccG)] - esd = SymbolicWedderburn.EigenSpaceDecomposition(F.(Ns[1])) - esd = SymbolicWedderburn.refine(esd, F.(Ns[2])) - esd = SymbolicWedderburn.refine(esd, F.(Ns[3])) + F = Characters.FiniteFields.GF{p} + Ns = [Characters.CMMatrix(ccG, i) for i = 1:length(ccG)] + esd = Characters.EigenSpaceDecomposition(F.(Ns[1])) + esd = Characters.refine(esd, F.(Ns[2])) + esd = Characters.refine(esd, F.(Ns[3])) @test isdiag(esd) - @test isdiag(SymbolicWedderburn.common_esd(Ns, F)) + @test isdiag(Characters.common_esd(Ns, F)) end # generic tests generictest_dixon_Fp(G) # tests specific to G - p = SymbolicWedderburn.dixon_prime(ccG) - F = SymbolicWedderburn.FiniteFields.GF{p} + p = Characters.dixon_prime(ccG) + F = Characters.FiniteFields.GF{p} - chars = SymbolicWedderburn.characters_dixon(F, ccG) + tbl = Characters.CharacterTable(F, G) + chars = Characters.irreducible_characters(tbl) - @test sort(degree.(chars)) == [1, 1, 1, 3] + @test sort(SymbolicWedderburn.degree.(chars)) == [1, 1, 1, 3] - @test Set(Int.(χ.vals) for χ in chars) == + @test Set(Int.(values(χ)) for χ in chars) == Set([[3, 0, 0, 6], [1, 4, 2, 1], [1, 2, 4, 1], [1, 1, 1, 1]]) end @@ -137,7 +141,7 @@ end end @testset "PowerMap" begin - pmG = SymbolicWedderburn.PowerMap(ccG) + pmG = Characters.PowerMap(ccG) @test size(pmG) == (4, 6) @test pmG[:, 0] == ones(Int, 4) @test pmG[:, 1] == 1:4 @@ -151,10 +155,10 @@ end # generic tests generictest_dixon_C(G) - chars_C = SymbolicWedderburn.characters_dixon(Int, G) - E = SymbolicWedderburn.Cyclotomics.E + chars_C = Characters.irreducible_characters(G, ccG) + E = Characters.Cyclotomics.E - @test [χ.vals for χ in chars_C] == [ + @test [collect(values(χ)) for χ in chars_C] == [ E(3, 0) .* [3, 0, 0, -1], E(3, 0) .* [1, E(3, 2), E(3, 1), 1], E(3, 0) .* [1, E(3, 1), E(3, 2), 1], @@ -176,10 +180,10 @@ end generictest_dixon_Fp(G) generictest_dixon_C(G) - chars = SymbolicWedderburn.characters_dixon(Int, ccG) + chars = Characters.irreducible_characters(G, ccG) - @test sort(degree.(chars)) == [1, 1, 2, 3, 3] - @test [χ.vals for χ in chars] == [ + @test sort(SymbolicWedderburn.degree.(chars)) == [1, 1, 2, 3, 3] + @test [collect(values(χ)) for χ in chars] == [ [2, 0, 2, -1, 0], [3, 1, -1, 0, -1], [1, 1, 1, 1, 1], @@ -203,10 +207,10 @@ end generictest_dixon_Fp(G) generictest_dixon_C(G) - chars = SymbolicWedderburn.characters_dixon(Int, ccG) + chars = Characters.irreducible_characters(G, ccG) - @test sort(degree.(chars)) == [1, 1, 1, 1, 4] - @test [χ.vals for χ in chars] == [ + @test sort(SymbolicWedderburn.degree.(chars)) == [1, 1, 1, 1, 4] + @test [collect(values(χ)) for χ in chars] == [ E(4, 0) .* [4, 0, 0, 0, -1], E(4, 0) .* [1, 1, 1, 1, 1], E(4, 0) .* [1, 1, -1, -1, 1], @@ -218,40 +222,20 @@ end @testset "Different base rings for characters_dixon" begin G = PermGroup(perm"(1,2,3,4)") - @test eltype(first(SymbolicWedderburn.characters_dixon(Float64, G))) <: - Cyclotomic{Float64} - @test eltype(first(SymbolicWedderburn.characters_dixon(Int, G))) <: - Cyclotomic{Int} - @test eltype( - first(SymbolicWedderburn.characters_dixon(Rational{Int}, G)), - ) <: Cyclotomic{Rational{Int}} + @test eltype(Characters.irreducible_characters(Rational{Int}, G)) <: + Characters.Character{<:Cyclotomic{Rational{Int}}} + @test eltype(Characters.irreducible_characters(Rational{BigInt}, G)) <: + Characters.Character{<:Cyclotomic{Rational{BigInt}}} end @time @testset "SmallPermGroups" begin for (ord, groups) in SmallPermGroups @testset "SmallGroup($ord, $n)" for (n, G) in enumerate(groups) - @test SymbolicWedderburn.characters_dixon(G) isa - Vector{<:SymbolicWedderburn.Character} + @test Characters.irreducible_characters(G) isa + Vector{<:Characters.Character} generictest_dixon_Fp(G) generictest_dixon_C(G) end end end end - -@testset "Characters io" begin - G = PermGroup([perm"(2,3)(4,5)"]) - chars = SymbolicWedderburn.characters_dixon(Int, G) - - @test sprint(show, chars[1]) == "[1, 1]" - @test sprint(show, chars[2]) == "[1, -1]" - - @test occursin( - "()^G → \t 1*E(1)^0\n(2,3)(4,5)^G → \t 1*E(1)^0", - sprint(show, MIME"text/plain"(), chars[1]), - ) - @test occursin( - "()^G → \t 1*E(1)^0\n(2,3)(4,5)^G → \t-1*E(1)^0", - sprint(show, MIME"text/plain"(), chars[2]) - ) -end diff --git a/test/eigenspacedecomposition.jl b/test/eigenspacedecomposition.jl index ffb56c0..dbb4887 100644 --- a/test/eigenspacedecomposition.jl +++ b/test/eigenspacedecomposition.jl @@ -1,5 +1,5 @@ function _row_echelon(A; L = []) - A, l = SymbolicWedderburn.row_echelon_form(A) + A, l = Characters.row_echelon_form(A) @test all([x ∈ l for x in L]) && length(L) <= length(l) for (i, j) in enumerate(l) @test all([iszero(A[k, j]) for k = 1:size(A, 1) if k != i]) @@ -8,14 +8,14 @@ function _row_echelon(A; L = []) end function _nullspace(A) - V = SymbolicWedderburn.left_nullspace(A) - W = SymbolicWedderburn.right_nullspace(A) + V = Characters.left_nullspace(A) + W = Characters.right_nullspace(A) @test all(iszero.(V*A)) @test all(iszero.(A*W)) end function _left_eigen(A) - E = SymbolicWedderburn.left_eigen(A) + E = Characters.left_eigen(A) for (val, space) in E @test space * A == space .* val end @@ -47,7 +47,7 @@ end 0 0 0 0 0 1 0 2; 0 0 0 0 0 0 0 0]) - @test SymbolicWedderburn.row_echelon_form(A)[1] == N + @test Characters.row_echelon_form(A)[1] == N end @testset "random" begin @@ -90,8 +90,8 @@ end 0 0 0 1 1 0 2 1; 0 0 0 0 0 1 0 2]) - K = SymbolicWedderburn.left_nullspace(A) - @test SymbolicWedderburn.row_echelon_form(K)[1] == N + K = Characters.left_nullspace(A) + @test Characters.row_echelon_form(K)[1] == N end @testset "random" begin @@ -120,17 +120,17 @@ end 0 0 0 3; 4 4 0 0; 0 0 4 0]) - - E2 = SymbolicWedderburn.left_eigen(M2) + + E2 = Characters.left_eigen(M2) @test size(E2[GF{7}(3)], 1) == 3 - E3 = SymbolicWedderburn.left_eigen(M3) - E4 = SymbolicWedderburn.left_eigen(M4) + E3 = Characters.left_eigen(M3) + E4 = Characters.left_eigen(M4) for val in GF{7}.([0, -3, -5, -6]) @test haskey(E3, val) @test haskey(E4, val) end - + @test E3[GF{7}(0)] == GF{7}.([1 2 0 0]) @test E3[GF{7}(-3)] == GF{7}.([1 1 1 1]) @test E3[GF{7}(-5)] == GF{7}.([1 1 2 4]) @@ -155,8 +155,8 @@ end 0 0 3 0; 0 0 0 4; 4 4 0 0]) - esd = SymbolicWedderburn.EigenSpaceDecomposition(M3) - @test SymbolicWedderburn.isdiag(esd) + esd = Characters.EigenSpaceDecomposition(M3) + @test Characters.isdiag(esd) #Handbook of Computational Group Theory p.262 M6 = GF{13}.([0 0 0 0 0 1; @@ -165,8 +165,8 @@ end 0 0 0 0 2 0; 0 0 0 2 0 0; 2 0 1 0 0 0]) - esd = SymbolicWedderburn.EigenSpaceDecomposition(M6) - + esd = Characters.EigenSpaceDecomposition(M6) + @test length(esd) == 4 @test esd[4] == GF{13}.([1 1 6 0 0 6]) @test esd[3] == GF{13}.([1 12 1 0 0 12; 0 0 0 1 12 0]) @@ -180,8 +180,8 @@ end 3 0 3 0 0 0; 0 0 0 0 2 0]) - esd = SymbolicWedderburn.refine(esd, M4) - @test SymbolicWedderburn.isdiag(esd) + esd = Characters.refine(esd, M4) + @test Characters.isdiag(esd) @test esd[1] == GF{13}.([1 12 6 0 0 7]) @test esd[2] == GF{13}.([1 1 1 1 1 1]) @test esd[3] == GF{13}.([1 1 1 12 12 1]) diff --git a/test/projections.jl b/test/projections.jl index b1369a5..59cbb92 100644 --- a/test/projections.jl +++ b/test/projections.jl @@ -1,13 +1,17 @@ function test_orthogonality(chars) - projs = [*(SymbolicWedderburn.matrix_projection(χ)...) for χ in chars] + projs = [SymbolicWedderburn.matrix_projection(χ) for χ in chars] - res = map(Iterators.product(projs, projs)) do (p, q) - if p == q - # p^2 == p - all(x -> isapprox(0.0, x; atol = 1e-12), p*p - p) - else - res = p * q - all(x -> isapprox(0.0, x; atol = 1e-12), p*q) + res = Matrix{Bool}(undef, length(projs), length(projs)) + + for j in axes(res, 2) + for i in axes(res, 1) + p, q = projs[i], projs[j] + if p == q + # p^2 == p + res[i,j] = all(x -> isapprox(0.0, x; atol = 1e-12), p*p - p) + else + res[i,j] = all(x -> isapprox(0.0, x; atol = 1e-12), p*q) + end end end return all(res) @@ -15,33 +19,38 @@ end @testset "Orthogonality of projections" begin G = PermGroup([perm"(1,2,3,4)"]) - chars = SymbolicWedderburn.characters_dixon(Float64, G) + chars = let irr = SymbolicWedderburn.irreducible_characters(G) + if !all(isreal, irr) + irr, _ = SymbolicWedderburn.affordable_real(irr) + end + SymbolicWedderburn.Character{Float64}.(irr) + end @test test_orthogonality(chars) - @test sum(first ∘ size, SymbolicWedderburn.isotypical_basis.(chars)) == - degree(G) + @test sum(first ∘ size, SymbolicWedderburn.image_basis.(chars)) == + PermutationGroups.degree(G) G = PermGroup([perm"(1,2,3,4)", perm"(1,2)"]) - chars = SymbolicWedderburn.characters_dixon(Float64, G) + chars = SymbolicWedderburn.irreducible_characters(G) @test test_orthogonality(chars) - @test sum(first ∘ size, SymbolicWedderburn.isotypical_basis.(chars)) == - degree(G) + @test sum(first ∘ size, SymbolicWedderburn.image_basis.(chars)) == + PermutationGroups.degree(G) G = PermGroup([perm"(1,2,3)", perm"(2,3,4)"]) - chars = SymbolicWedderburn.characters_dixon(Float64, G) + chars = SymbolicWedderburn.irreducible_characters(G) @test test_orthogonality(chars) - @test sum(first ∘ size, SymbolicWedderburn.isotypical_basis.(chars)) == - degree(G) + @test sum(first ∘ size, SymbolicWedderburn.image_basis.(chars)) == + PermutationGroups.degree(G) @time for ord = 2:16 # @testset "SmallGroup($ord, $n)" for (n, G) in enumerate(SmallPermGroups[ord]) - chars = SymbolicWedderburn.characters_dixon(Rational{Int}, G) + chars = SymbolicWedderburn.irreducible_characters(Rational{Int}, G) @test test_orthogonality(chars) @test sum( first ∘ size, - SymbolicWedderburn.isotypical_basis.(chars), - ) == degree(G) + SymbolicWedderburn.image_basis.(chars), + ) == PermutationGroups.degree(G) end end end diff --git a/test/runtests.jl b/test/runtests.jl index c5b972d..0449c3f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,15 +11,19 @@ using Cyclotomics include("smallgroups.jl") -include("gf.jl") -include("eigenspacedecomposition.jl") -include("ccmatrix.jl") -include("dixon.jl") -include("characters.jl") +@testset "Characters" begin + import SymbolicWedderburn.Characters + include("gf.jl") + include("eigenspacedecomposition.jl") + include("ccmatrix.jl") + include("dixon.jl") + include("characters.jl") +end include("projections.jl") include("sa_basis.jl") include("action_permutation.jl") include("action_linear.jl") +include("action_dihedral.jl") if VERSION >= v"1.6.0" && !haskey(ENV, "CI") @testset "Examples" begin @@ -27,7 +31,8 @@ if VERSION >= v"1.6.0" && !haskey(ENV, "CI") Pkg.activate(joinpath(@__DIR__, "..", "examples")) Pkg.instantiate() include("../examples/ex_C2_linear.jl") - include("../examples/ex_C4.jl") + include("../examples/ex_S4.jl") include("../examples/ex_motzkin.jl") + include("../examples/ex_robinson_form.jl") end end diff --git a/test/sa_basis.jl b/test/sa_basis.jl index 1beaffb..807450d 100644 --- a/test/sa_basis.jl +++ b/test/sa_basis.jl @@ -1,83 +1,140 @@ +using StarAlgebras + @testset "affordable real degrees/dot" begin G = SmallPermGroups[10][2] # C₂⊕C₅ + tbl = SymbolicWedderburn.CharacterTable(Rational{Int}, G) + chars = SymbolicWedderburn.irreducible_characters(tbl) - chars = SymbolicWedderburn.characters_dixon(G) - charsR = SymbolicWedderburn.affordable_real(chars) - - @test all(isone ∘ degree, chars) + @test all(isone ∘ SymbolicWedderburn.degree, chars) @test all(χ -> isone(dot(χ, χ)), chars) chars_fl = SymbolicWedderburn.Character{ComplexF64}.(chars) - @test all(isone ∘ degree, chars_fl) + @test all(isone ∘ SymbolicWedderburn.degree, chars_fl) @test all(χ -> isapprox(dot(χ, χ), 1), chars_fl) - @test degree.(charsR) == [dot(χ, χ) for χ in charsR] == [1,1,2,2,2,2] + charsR, _ = SymbolicWedderburn.affordable_real(chars) + @test SymbolicWedderburn.degree.(charsR) == + [dot(χ, χ) for χ in charsR] == + [1,1,2,2,2,2] charsR_fl = SymbolicWedderburn.Character{Float64}.(charsR) - @test all(degree.(charsR_fl) .≈ [dot(χ, χ) for χ in charsR_fl] .≈ [1,1,2,2,2,2]) + @test all([1,1,2,2,2,2] .== + SymbolicWedderburn.degree.(charsR_fl) .≈ + [dot(χ, χ) for χ in charsR_fl] + ) end @testset "Symmetry adapted basis" begin + @testset "step by step" begin + G = PermGroup(perm"(1,2)", perm"(1,2,3)") + irr = SymbolicWedderburn.irreducible_characters(G) + @test irr isa AbstractVector{<:SymbolicWedderburn.Character{<:Cyclotomic}} + @test SymbolicWedderburn.degree.(irr) == [2,1,1] + + RG = let G = G + b = StarAlgebras.Basis{UInt16}(collect(G)) + StarAlgebra(G, b, (length(b), length(b))) + end + + µ, s = SymbolicWedderburn.minimal_rank_projection(irr[1], RG) + @test µ isa AlgebraElement + @test s + + @test SymbolicWedderburn.Character{Rational{Int}}(irr[1]) isa + SymbolicWedderburn.Character{Rational{Int}} + @test collect(values(SymbolicWedderburn.Character{Float64}(irr[1]))) isa Vector{Float64} + + mps, simple = SymbolicWedderburn.minimal_projection_system(irr, RG) + @test all(simple) + + @test rank(float.(SymbolicWedderburn.matrix_projection(irr[1]))) == 2 + @test rank(float.(SymbolicWedderburn.matrix_projection(mps[1]))) == 1 + + sa_basis_ssimple = symmetry_adapted_basis(Rational{Int}, G, Rational{Int}, semisimple=true) + + @test issimple.(sa_basis_ssimple) == [false, true] + @test rank.(convert.(Matrix, sa_basis_ssimple)) == [2, 1] + @test dot( + multiplicity.(sa_basis_ssimple), + SymbolicWedderburn.degree.(sa_basis_ssimple) + ) == PermutationGroups.degree(G) + + sa_basis = symmetry_adapted_basis(Rational{Int}, G, Rational{Int}, semisimple=false) + @test all(issimple.(sa_basis)) + @test rank.(convert.(Matrix{Float64}, sa_basis)) == [1,1] + end + C₄ = PermGroup([perm"(1,2,3,4)"]) A₄ = PermGroup([perm"(1,2,3)", perm"(2,3,4)"]) S₄ = PermGroup([perm"(1,2,3,4)", perm"(1,2)"]) - basis = symmetry_adapted_basis(C₄) - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basis) == degree(C₄) - basis = symmetry_adapted_basis(A₄) - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basis) == degree(A₄) - basis = symmetry_adapted_basis(S₄) - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basis) == degree(S₄) - - basis = symmetry_adapted_basis(ComplexF64, C₄) - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basis) == degree(C₄) - basis = symmetry_adapted_basis(ComplexF64, A₄) - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basis) == degree(A₄) - basis = symmetry_adapted_basis(ComplexF64, S₄) - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basis) == degree(S₄) - - basis = symmetry_adapted_basis(Float64, C₄) - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basis) == degree(C₄) - basis = symmetry_adapted_basis(Float64, A₄) - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basis) == degree(A₄) - basis = symmetry_adapted_basis(Float64, S₄) - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basis) == degree(S₄) + sa_basis = symmetry_adapted_basis(C₄, semisimple=true) + @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(C₄) + sa_basis = symmetry_adapted_basis(A₄, semisimple=true) + @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(A₄) + sa_basis = symmetry_adapted_basis(S₄, semisimple=true) + @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(S₄) + + sa_basis = symmetry_adapted_basis(ComplexF64, C₄, semisimple=true) + @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(C₄) + sa_basis = symmetry_adapted_basis(ComplexF64, A₄, semisimple=true) + @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(A₄) + sa_basis = symmetry_adapted_basis(ComplexF64, S₄, semisimple=true) + @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(S₄) + + sa_basis = symmetry_adapted_basis(Float64, C₄, semisimple=true) + @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(C₄) + sa_basis = symmetry_adapted_basis(Float64, A₄, semisimple=true) + @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(A₄) + sa_basis = symmetry_adapted_basis(Float64, S₄, semisimple=true) + @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(S₄) end - @testset "Symmetry adapted basis: small groups" begin @time @testset "Rational" begin for ord = 2:30 @testset "SmallGroup($ord, $n)" for (n, G) in enumerate(SmallPermGroups[ord]) - @testset "complex" begin - S = if (ord, n) in ((21,1),) - Rational{BigInt} - else - Rational{Int} - end + S = Rational{Int} - basis = symmetry_adapted_basis(G, S) + sa_basis = symmetry_adapted_basis(G, S, semisimple=true) - @test dot(SymbolicWedderburn.degree.(basis), SymbolicWedderburn.multiplicity.(basis)) == degree(G) + @test dot(SymbolicWedderburn.degree.(sa_basis), multiplicity.(sa_basis)) == + PermutationGroups.degree(G) - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basis) == degree(G) - end + @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(G) - @testset "real" begin - S = if (ord, n) in ((17,1), (19,1), (21,1), (23,1), (24,3), (25,1), (29,1)) + S = if (ord, n) in ((21,1),) Rational{BigInt} else Rational{Int} end + sa_basis = symmetry_adapted_basis(G, S, semisimple=false) + for b in sa_basis + if issimple(b) + @test multiplicity(b) == size(b, 1) + else + @test multiplicity(b)*SymbolicWedderburn.degree(b) == size(b, 1) + end + end + end - chars = SymbolicWedderburn.affordable_real( - SymbolicWedderburn.characters_dixon(S, G) - ) - - basisR = symmetry_adapted_basis(chars) - @test dot(SymbolicWedderburn.degree.(basisR), SymbolicWedderburn.multiplicity.(basisR)) == degree(G) - - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basisR) == degree(G) + @testset "real" begin + S = Rational{Int} + + sa_basisR = symmetry_adapted_basis(Float64, G, S, semisimple=true) + @test dot(SymbolicWedderburn.degree.(sa_basisR), multiplicity.(sa_basisR)) == + PermutationGroups.degree(G) + @test sum(first ∘ size, sa_basisR) == PermutationGroups.degree(G) + + sa_basisR = symmetry_adapted_basis(Float64, G, S, semisimple=false) + for b in sa_basisR + if issimple(b) + @test multiplicity(b) == size(b, 1) || 2*multiplicity(b) == size(b, 1) + # the first condiditon doesn't hold for realified characters; + else + @test multiplicity(b)*SymbolicWedderburn.degree(b) == size(b, 1) + end + end end end end @@ -88,17 +145,17 @@ end @testset "SmallGroup($ord, $n)" for (n, G) in enumerate(SmallPermGroups[ord]) T = Float64 - basis = symmetry_adapted_basis(T, G) + sa_basis = symmetry_adapted_basis(T, G, semisimple=true) + @test sa_basis isa Vector{<:AbstractMatrix{T}} - @test all(basis) do b - B = SymbolicWedderburn.basis(b) - rank(B) == size(B, 1) + @test all(sa_basis) do b + rank(b) == size(b, 1) end - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basis) == degree(G) + @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(G) - basisC = symmetry_adapted_basis(Complex{T}, G) - @test sum(first ∘ size ∘ SymbolicWedderburn.basis, basisC) == degree(G) + basisC = symmetry_adapted_basis(Complex{T}, G, semisimple=true) + @test sum(first ∘ size, basisC) == PermutationGroups.degree(G) end end end