diff --git a/Project.toml b/Project.toml index e681151..eea7135 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1" @@ -18,6 +19,7 @@ Cyclotomics = "0.3" GroupsCore = "0.4" PermutationGroups = "0.4.2" PrecompileTools = "1" +PrettyTables = "2" Primes = "0.4, 0.5" StarAlgebras = "0.2" julia = "1.6" diff --git a/src/Characters/Characters.jl b/src/Characters/Characters.jl index 8caca5d..41a0cc8 100644 --- a/src/Characters/Characters.jl +++ b/src/Characters/Characters.jl @@ -26,5 +26,8 @@ include("cmmatrix.jl") include("powermap.jl") include("character_tables.jl") include("class_functions.jl") +include("chars.jl") +include("io.jl") + include("dixon.jl") end diff --git a/src/Characters/character_tables.jl b/src/Characters/character_tables.jl index fb3a437..f84e894 100644 --- a/src/Characters/character_tables.jl +++ b/src/Characters/character_tables.jl @@ -162,86 +162,3 @@ function normalize!(chtbl::CharacterTable{<:Group,<:FiniteFields.GF}) 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/chars.jl b/src/Characters/chars.jl new file mode 100644 index 0000000..90fac44 --- /dev/null +++ b/src/Characters/chars.jl @@ -0,0 +1,217 @@ +""" + Character <: AbstractClassFunction +Struct representing (possibly virtual) character of a group. + +Characters are backed by `table(χ)::CharacterTable` which actually stores the +character values. The multiplicities (decomposition into the irreducible +summands) of a given character can be obtained by calling `multiplicities(χ)` +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 + multips::Vector{S} +end + +function Character{R}( + chtbl::CharacterTable{Gr,T}, + multips::AbstractVector{S}, +) where {R,Gr,T,S} + return Character{R,S,typeof(chtbl)}(chtbl, multips) +end + +function Character(chtbl::CharacterTable, multips::AbstractVector) + R = Base._return_type(*, Tuple{eltype(chtbl),eltype(multips)}) + @assert R ≠ Any + return Character{R}(chtbl, multips) +end + +function Character(chtbl::CharacterTable, i::Integer) + return Character{eltype(chtbl)}(chtbl, i) +end + +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(multiplicities(χ)) + ChT = typeof(table(χ)) + return Character{T,S,ChT}(table(χ), multiplicities(χ)) +end + +## Accessors +table(χ::Character) = χ.table +multiplicities(χ::Character) = χ.multips + +## AbstractClassFunction api +Base.parent(χ::Character) = parent(table(χ)) +conjugacy_classes(χ::Character) = conjugacy_classes(table(χ)) +function Base.values(χ::Character{T}) where {T} + return T[χ[i] for i in 1:nconjugacy_classes(table(χ))] +end + +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 convert( + T, + sum( + c * table(χ)[idx, i] for + (idx, c) in enumerate(multiplicities(χ)) if !iszero(c); + init = zero(T), + ), + ) +end + +## Basic functionality + +function Base.:(==)(χ::Character, ψ::Character) + return table(χ) === table(ψ) && multiplicities(χ) == multiplicities(ψ) +end +Base.hash(χ::Character, h::UInt) = hash(table(χ), hash(multiplicities(χ), h)) + +function Base.deepcopy_internal(χ::Character{T}, d::IdDict) where {T} + haskey(d, χ) && return d[χ] + return Character{T}(table(χ), copy(multiplicities(χ))) +end + +## Character arithmetic + +for f in (:+, :-) + @eval begin + function Base.$f(χ::Character, ψ::Character) + @assert table(χ) === table(ψ) + return Character(table(χ), $f(multiplicities(χ), multiplicities(ψ))) + end + end +end + +Base.:*(χ::Character, c::Number) = Character(table(χ), c .* multiplicities(χ)) +Base.:*(c::Number, χ::Character) = χ * c +Base.:/(χ::Character, c::Number) = Character(table(χ), multiplicities(χ) ./ c) + +Base.zero(χ::Character) = 0 * χ + +function __decompose(T::Type, values::AbstractVector, tbl::CharacterTable) + ψ = ClassFunction(values, conjugacy_classes(tbl), tbl.inv_of) + return decompose(T, ψ, tbl) +end + +function decompose(cfun::AbstractClassFunction, tbl::CharacterTable) + return decompose(eltype(cfun), cfun, tbl) +end + +function decompose(T::Type, cfun::AbstractClassFunction, tbl::CharacterTable) + vals = Vector{T}(undef, length(conjugacy_classes(tbl))) + return decompose!(vals, cfun, tbl) +end + +function decompose!( + vals::AbstractVector, + cfun::AbstractClassFunction, + tbl::CharacterTable, +) + @assert length(vals) == length(conjugacy_classes(tbl)) + @assert conjugacy_classes(tbl) === conjugacy_classes(cfun) + + for (i, idx) in enumerate(eachindex(vals)) + χ = Character(tbl, i) + vals[idx] = dot(χ, cfun) + end + return vals +end + +function Base.:*(χ::Character, ψ::Character) + @assert table(χ) == table(ψ) + values = Characters.values(χ) .* Characters.values(ψ) + return Character(table(χ), __decompose(Int, values, table(χ))) +end + +Base.:^(χ::Character, n::Integer) = Base.power_by_squaring(χ, n) + +## Group-theoretic functions: + +PermutationGroups.degree(χ::Character) = Int(χ(one(parent(χ)))) +function PermutationGroups.degree( + χ::Character{T,CCl}, +) where {T,CCl<:AbstractOrbit{<:AbstractMatrix}} + return Int(χ[1]) +end + +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) + multips = S[dot(ψ, χ) for χ in irreducible_characters(tbl)] + return Character{T,eltype(multips),typeof(tbl)}(tbl, multips) +end + +function isvirtual(χ::Character) + return any(<(0), multiplicities(χ)) || any(!isinteger, multiplicities(χ)) +end + +function isirreducible(χ::Character) + C = multiplicities(χ) + k = findfirst(!iszero, C) + k !== nothing || return false # χ is zero + isone(C[k]) || return false # muliplicity is ≠ 1 + kn = findnext(!iszero, C, k + 1) + kn === nothing && return true # there is only one ≠ 0 entry + return false +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 + χ.multips .+= multiplicities(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 diff --git a/src/Characters/class_functions.jl b/src/Characters/class_functions.jl index 307d8db..73d71ea 100644 --- a/src/Characters/class_functions.jl +++ b/src/Characters/class_functions.jl @@ -54,239 +54,3 @@ Base.@propagate_inbounds function Base.getindex(χ::ClassFunction, i::Integer) return values(χ)[i] end -""" - Character <: AbstractClassFunction -Struct representing (possibly virtual) character of a group. - -Characters are backed by `table(χ)::CharacterTable` which actually stores the -character values. The multiplicities (decomposition into the irreducible -summands) of a given character can be obtained by calling `multiplicities(χ)` -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 - multips::Vector{S} -end - -function Character{R}( - chtbl::CharacterTable{Gr,T}, - multips::AbstractVector{S}, -) where {R,Gr,T,S} - return Character{R,S,typeof(chtbl)}(chtbl, multips) -end - -function Character(chtbl::CharacterTable, multips::AbstractVector) - R = Base._return_type(*, Tuple{eltype(chtbl),eltype(multips)}) - @assert R ≠ Any - return Character{R}(chtbl, multips) -end - -function Character(chtbl::CharacterTable, i::Integer) - return Character{eltype(chtbl)}(chtbl, i) -end - -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(multiplicities(χ)) - ChT = typeof(table(χ)) - return Character{T,S,ChT}(table(χ), multiplicities(χ)) -end - -## Accessors -table(χ::Character) = χ.table -multiplicities(χ::Character) = χ.multips - -## AbstractClassFunction api -Base.parent(χ::Character) = parent(table(χ)) -conjugacy_classes(χ::Character) = conjugacy_classes(table(χ)) -function Base.values(χ::Character{T}) where {T} - return T[χ[i] for i in 1:nconjugacy_classes(table(χ))] -end - -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 convert( - T, - sum( - c * table(χ)[idx, i] for - (idx, c) in enumerate(multiplicities(χ)) if !iszero(c); - init = zero(T), - ), - ) -end - -## Basic functionality - -function Base.:(==)(χ::Character, ψ::Character) - return table(χ) === table(ψ) && multiplicities(χ) == multiplicities(ψ) -end -Base.hash(χ::Character, h::UInt) = hash(table(χ), hash(multiplicities(χ), h)) - -function Base.deepcopy_internal(χ::Character{T}, d::IdDict) where {T} - haskey(d, χ) && return d[χ] - return Character{T}(table(χ), copy(multiplicities(χ))) -end - -## Character arithmetic - -for f in (:+, :-) - @eval begin - function Base.$f(χ::Character, ψ::Character) - @assert table(χ) === table(ψ) - return Character(table(χ), $f(multiplicities(χ), multiplicities(ψ))) - end - end -end - -Base.:*(χ::Character, c::Number) = Character(table(χ), c .* multiplicities(χ)) -Base.:*(c::Number, χ::Character) = χ * c -Base.:/(χ::Character, c::Number) = Character(table(χ), multiplicities(χ) ./ c) - -Base.zero(χ::Character) = 0 * χ - -function __decompose(T::Type, values::AbstractVector, tbl::CharacterTable) - ψ = ClassFunction(values, conjugacy_classes(tbl), tbl.inv_of) - return decompose(T, ψ, tbl) -end - -function decompose(cfun::AbstractClassFunction, tbl::CharacterTable) - return decompose(eltype(cfun), cfun, tbl) -end - -function decompose(T::Type, cfun::AbstractClassFunction, tbl::CharacterTable) - vals = Vector{T}(undef, length(conjugacy_classes(tbl))) - return decompose!(vals, cfun, tbl) -end - -function decompose!( - vals::AbstractVector, - cfun::AbstractClassFunction, - tbl::CharacterTable, -) - @assert length(vals) == length(conjugacy_classes(tbl)) - @assert conjugacy_classes(tbl) === conjugacy_classes(cfun) - - for (i, idx) in enumerate(eachindex(vals)) - χ = Character(tbl, i) - vals[idx] = dot(χ, cfun) - end - return vals -end - -function Base.:*(χ::Character, ψ::Character) - @assert table(χ) == table(ψ) - values = Characters.values(χ) .* Characters.values(ψ) - return Character(table(χ), __decompose(Int, values, table(χ))) -end - -Base.:^(χ::Character, n::Integer) = Base.power_by_squaring(χ, n) - -## Group-theoretic functions: - -PermutationGroups.degree(χ::Character) = Int(χ(one(parent(χ)))) -function PermutationGroups.degree( - χ::Character{T,CCl}, -) where {T,CCl<:AbstractOrbit{<:AbstractMatrix}} - return Int(χ[1]) -end - -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) - multips = S[dot(ψ, χ) for χ in irreducible_characters(tbl)] - return Character{T,eltype(multips),typeof(tbl)}(tbl, multips) -end - -function isvirtual(χ::Character) - return any(<(0), multiplicities(χ)) || any(!isinteger, multiplicities(χ)) -end - -function isirreducible(χ::Character) - C = multiplicities(χ) - k = findfirst(!iszero, C) - k !== nothing || return false # χ is zero - isone(C[k]) || return false # muliplicity is ≠ 1 - kn = findnext(!iszero, C, k + 1) - kn === nothing && return true # there is only one ≠ 0 entry - return false -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 - χ.multips .+= multiplicities(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(χ)) - return _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(multiplicities(χ)) - 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/Characters/io.jl b/src/Characters/io.jl new file mode 100644 index 0000000..29f9abe --- /dev/null +++ b/src/Characters/io.jl @@ -0,0 +1,75 @@ +import PrettyTables +import PrettyTables.Tables + +Tables.istable(::Type{<:CharacterTable}) = true +Tables.rowaccess(::Type{<:CharacterTable}) = true +Tables.rows(chtbl::CharacterTable) = irreducible_characters(chtbl) +Tables.getcolumn(row::Character, i::Int) = row[i] + +# ugly hack to get this going +function Tables.columnnames(char::Character) + return [Symbol(i) for i in axes(conjugacy_classes(char), 1)] +end +function Tables.getcolumn(row::Character, nm::Symbol) + i = parse(Int, string(nm)) + return Tables.getcolumn(row, i) +end + +__coefs_desc(T::Type) = T +__coefs_desc(::Type{<:Rational{T}}) where {T} = "rationals ($T)" +__coefs_desc(::Type{<:Cyclotomics.Cyclotomic{T}}) where {T} = "cyclotomics ($T)" + +function Base.show(io::IO, ::MIME"text/plain", chtbl::CharacterTable) + hl_odd = PrettyTables.Highlighter(; + f = (rule, i, j) -> i % 2 == 0, + crayon = PrettyTables.Crayon(; + foreground = :dark_gray, + negative = true, + ), + ) + + fmt = (v, args...) -> sprint(show, MIME"text/plain"(), v) + + return PrettyTables.pretty_table( + io, + chtbl; + title = "Character table of $(parent(chtbl)) over $(__coefs_desc(eltype(chtbl)))", + header = ["$(first(cc))^G" for cc in conjugacy_classes(chtbl)], + row_labels = [ + Symbol('χ', FiniteFields.subscriptify(i)) for i in axes(chtbl, 1) + ], + row_label_column_title = "", + # hlines = [:header, :end], + # vlines = [1], + formatters = fmt, + autowrap = true, + linebreaks = true, + columns_width = displaysize(io)[2] ÷ size(chtbl, 2) - 8, + reserved_display_lines = 3[], + # vcrop_mode = :middle, + # equal_columns_width = true, + # crop = :vertical, + ellipsis_line_skip = 1, + # alignment = [:r, :l], + highlighters = hl_odd, + ) +end + +function Base.show(io::IO, ::MIME"text/plain", χ::Character) + println(io, "Character over ", __coefs_desc(eltype(χ))) + return _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(multiplicities(χ)) + 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/SymbolicWedderburn.jl b/src/SymbolicWedderburn.jl index bb400c5..4c4f913 100644 --- a/src/SymbolicWedderburn.jl +++ b/src/SymbolicWedderburn.jl @@ -27,6 +27,7 @@ include("actions.jl") include("group_action_error.jl") include("action_characters.jl") include("matrix_projections.jl") +include("image_basis.jl") include("minimal_projections.jl") include("direct_summands.jl") include("sa_basis.jl") diff --git a/src/image_basis.jl b/src/image_basis.jl new file mode 100644 index 0000000..389f800 --- /dev/null +++ b/src/image_basis.jl @@ -0,0 +1,127 @@ +## Finding basis of the row-space (right image) of an AbstractMatrix + +""" + image_basis(A::AbstractMatrix) + image_basis([hom::InducedActionHomomorphism, ]χ::Character[, rank]) + image_basis([hom::InducedActionHomomorphism, ]α::AlgebraElement[, rank]) +Return basis of the row-space of a matrix. + +For characters or algebra elements return a basis of the row-space of +`matrix_projection([hom, ]χ)` or `matrix_representation([hom, ]α)`. + +Two methods are employed to achieve the goal. +* By default (symbolic, exact) row echelon form is produced, and therefore there +is no guarantee on the orthogonality of the returned basis vectors (rows). +* If `eltype(A) <: AbstractFloat` a thin `svd` (or `qr` in the sparse case) +decomposition is computed and the appropriate rows of its orthogonal factor are +returned, thus the basis is (numerically) orthonormal. + +If the dimension (`=rank` of the matrix) is known beforehand (because of e.g. +algebra) it can be passed to `image_basis` to use more efficient methods. + +# Examples: +```julia +julia> a = rand(-3:3, 3,3).//rand(1:4, 3,3); + +julia> a[3, :] .= 2a[1, :] .- 1a[2, :]; a +3×3 Matrix{Rational{Int64}}: + 3//2 0//1 1//1 + -1//2 -2//1 1//2 + 7//2 2//1 3//2 + +julia> ib = SymbolicWedderburn.image_basis(a) +2×3 Matrix{Rational{Int64}}: + 1//1 0//1 2//3 + 0//1 1//1 -5//12 + +julia> ibf = SymbolicWedderburn.image_basis(float.(a)) +2×3 Matrix{Float64}: + -0.666651 0.416657 -0.618041 + 0.529999 0.847998 -8.85356e-17 + +``` +""" +image_basis(A::AbstractMatrix, rank = nothing) = image_basis!(deepcopy(A), rank) + +function image_basis(α::Union{Character,AlgebraElement}, rank = nothing) + mpr = matrix_representation(α) + return image_basis!(mpr, rank) +end + +function image_basis( + hom::InducedActionHomomorphism, + α::Union{Character,AlgebraElement}, + rank = nothing, +) + mpr = matrix_representation(hom, α) + return image_basis!(mpr, rank) +end + +## +# the internal, (possibly) modifying versions + +_eps(T::Type{<:AbstractFloat}) = eps(T) +_eps(::Type{Complex{T}}) where {T} = 2 * _eps(real(T)) +_eps(::Type{<:Cyclotomic{T}}) where {T} = 2_eps(T) +_eps(::Any) = 0.0 + +function image_basis!(A::AbstractMatrix, rank = nothing) + img, pivots = row_echelon_form!(A) + if img isa AbstractSparseArray + ε = _eps(eltype(img)) + if iszero(ε) + dropzeros!(img) + else + droptol!(img, _eps(eltype(img)) * max(size(img)...)) + end + end + # TODO: orthogonalize the result + return img[pivots, :] +end + +function image_basis!( + A::AbstractMatrix{T}, + ::Nothing, +) where {T<:Union{AbstractFloat,Complex}} + F = svd!(A) + tol = _eps(T) * first(F.S) + rk = count(x -> x > tol, F.S) + return Matrix((@view F.U[1:rk, :])') +end + +function image_basis!( + A::AbstractSparseMatrix{T}, + ::Nothing, +) where {T<:Union{AbstractFloat,Complex}} + F = qr(A) + rk = rank(F) + img = let tmp = F.Q * Matrix(I, size(A, 1), rk) + pinv = getfield(F, :rpivinv) + sparse((@view tmp[pinv, :])') + end + return droptol!(img, _eps(T) * max(size(img)...)) +end + +# to disambiguate +function image_basis!( + M::AbstractMatrix{T}, + rank::Integer, +) where {T<:Union{AbstractFloat,Complex}} + img = image_basis!(M, nothing) + if size(img, 1) ≠ rank + @warn "Possibly wrong numerical rank: (numerical) $(size(img, 1)) ≠ $rank (expected)" + end + return img +end + +function image_basis!( + M::AbstractSparseMatrix{T}, + rank::Integer, +) where {T<:Union{AbstractFloat,Complex}} + # return _image_basis!(M, rank) + img = image_basis!(M, nothing) + if size(img, 1) ≠ rank + @warn "Possibly wrong rank estimate? (numerical) $(size(img, 1)) ≠ $rank (expected)" + end + return img +end diff --git a/src/matrix_projections.jl b/src/matrix_projections.jl index f6133d3..42afb33 100644 --- a/src/matrix_projections.jl +++ b/src/matrix_projections.jl @@ -360,131 +360,3 @@ function matrix_representation(hom::InducedActionHomomorphism, χ::Character) return isirreducible(χ) ? matrix_projection_irr(hom, χ) : matrix_projection(hom, χ) end - -## Finding basis of the row-space (right image) of an AbstractMatrix - -""" - image_basis(A::AbstractMatrix) - image_basis([hom::InducedActionHomomorphism, ]χ::Character[, rank]) - image_basis([hom::InducedActionHomomorphism, ]α::AlgebraElement[, rank]) -Return basis of the row-space of a matrix. - -For characters or algebra elements return a basis of the row-space of -`matrix_projection([hom, ]χ)` or `matrix_representation([hom, ]α)`. - -Two methods are employed to achieve the goal. -* By default (symbolic, exact) row echelon form is produced, and therefore there -is no guarantee on the orthogonality of the returned basis vectors (rows). -* If `eltype(A) <: AbstractFloat` a thin `svd` (or `qr` in the sparse case) -decomposition is computed and the appropriate rows of its orthogonal factor are -returned, thus the basis is (numerically) orthonormal. - -If the dimension (`=rank` of the matrix) is known beforehand (because of e.g. -algebra) it can be passed to `image_basis` to use more efficient methods. - -# Examples: -```julia -julia> a = rand(-3:3, 3,3).//rand(1:4, 3,3); - -julia> a[3, :] .= 2a[1, :] .- 1a[2, :]; a -3×3 Matrix{Rational{Int64}}: - 3//2 0//1 1//1 - -1//2 -2//1 1//2 - 7//2 2//1 3//2 - -julia> ib = SymbolicWedderburn.image_basis(a) -2×3 Matrix{Rational{Int64}}: - 1//1 0//1 2//3 - 0//1 1//1 -5//12 - -julia> ibf = SymbolicWedderburn.image_basis(float.(a)) -2×3 Matrix{Float64}: - -0.666651 0.416657 -0.618041 - 0.529999 0.847998 -8.85356e-17 - -``` -""" -image_basis(A::AbstractMatrix, rank = nothing) = image_basis!(deepcopy(A), rank) - -function image_basis(α::Union{Character,AlgebraElement}, rank = nothing) - mpr = matrix_representation(α) - return image_basis!(mpr, rank) -end - -function image_basis( - hom::InducedActionHomomorphism, - α::Union{Character,AlgebraElement}, - rank = nothing, -) - mpr = matrix_representation(hom, α) - return image_basis!(mpr, rank) -end - -## -# the internal, (possibly) modifying versions - -_eps(T::Type{<:AbstractFloat}) = eps(T) -_eps(::Type{Complex{T}}) where {T} = 2 * _eps(real(T)) -_eps(::Type{<:Cyclotomic{T}}) where {T} = 2_eps(T) -_eps(::Any) = 0.0 - -function image_basis!(A::AbstractMatrix, rank = nothing) - img, pivots = row_echelon_form!(A) - if img isa AbstractSparseArray - ε = _eps(eltype(img)) - if iszero(ε) - dropzeros!(img) - else - droptol!(img, _eps(eltype(img)) * max(size(img)...)) - end - end - # TODO: orthogonalize the result - return img[pivots, :] -end - -function image_basis!( - A::AbstractMatrix{T}, - ::Nothing, -) where {T<:Union{AbstractFloat,Complex}} - F = svd!(A) - tol = _eps(T) * first(F.S) - rk = count(x -> x > tol, F.S) - return Matrix((@view F.U[1:rk, :])') -end - -function image_basis!( - A::AbstractSparseMatrix{T}, - ::Nothing, -) where {T<:Union{AbstractFloat,Complex}} - F = qr(A) - rk = rank(F) - img = let tmp = F.Q * Matrix(I, size(A, 1), rk) - pinv = getfield(F, :rpivinv) - sparse((@view tmp[pinv, :])') - end - return droptol!(img, _eps(T) * max(size(img)...)) -end - -# to disambiguate -function image_basis!( - M::AbstractMatrix{T}, - rank::Integer, -) where {T<:Union{AbstractFloat,Complex}} - img = image_basis!(M, nothing) - if size(img, 1) ≠ rank - @warn "Possibly wrong numerical rank: (numerical) $(size(img, 1)) ≠ $rank (expected)" - end - return img -end - -function image_basis!( - M::AbstractSparseMatrix{T}, - rank::Integer, -) where {T<:Union{AbstractFloat,Complex}} - # return _image_basis!(M, rank) - img = image_basis!(M, nothing) - if size(img, 1) ≠ rank - @warn "Possibly wrong rank estimate? (numerical) $(size(img, 1)) ≠ $rank (expected)" - end - return img -end diff --git a/test/characters.jl b/test/characters.jl index 74f0cff..36f5977 100644 --- a/test/characters.jl +++ b/test/characters.jl @@ -102,5 +102,12 @@ end @test sprint(show, -2chars[1] + 3chars[2]) == "-2·χ₁ +3·χ₂" @test sprint(show, MIME"text/plain"(), chars[1]) == - "Character over Cyclotomic{Rational{$Int}, SparseVector{Rational{$Int}, $Int}}\nχ₁" + "Character over cyclotomics (Rational{$Int})\nχ₁" + + χ = Characters.Character{Rational{BigInt}}(chars[1]) + @test sprint(show, MIME"text/plain"(), χ) == + "Character over rationals ($BigInt)\nχ₁" + + @test sprint(show, MIME"text/plain"(), Characters.table(χ)) isa + AbstractString end