diff --git a/Project.toml b/Project.toml index 1002b26..dfd4def 100644 --- a/Project.toml +++ b/Project.toml @@ -1,26 +1,14 @@ name = "StarAlgebras" uuid = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1" authors = ["Marek Kaluba "] -version = "0.2.1" +version = "0.3.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] -Groups = "0.8" -GroupsCore = "0.5" -PermutationGroups = "0.6" +MutableArithmetics = "1.4.3" SparseArrays = "1" julia = "1.6" - -[extras] -Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557" -GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120" -PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420" -Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[targets] -test = ["Pkg", "Groups", "GroupsCore", "PermutationGroups", "Random", "Test"] diff --git a/README.md b/README.md index 74ed129..6b8773c 100644 --- a/README.md +++ b/README.md @@ -7,10 +7,12 @@ ---- -The package implements `*`-algebras with basis. The prime example use is group/monoid algebras (or rings) (or their finite dimensional subspaces). An example usage can be as follows. +The package implements `*`-algebras with basis. The prime example use are group +or monoid algebras (or their finite dimensional subspaces). +An example usage can be as follows. ```julia -julia> using StarAlgebras +julia> using StarAlgebras; import StarAlgebras as SA julia> using PermutationGroups @@ -19,21 +21,15 @@ Permutation group on 2 generators generated by (1,2) (1,2,3) -julia> b = StarAlgebras.Basis{UInt8}(collect(G)) -6-element StarAlgebras.Basis{Permutation{Int64, …}, UInt8, Vector{Permutation{Int64, …}}}: - () - (2,3) - (1,2) - (1,3,2) - (1,3) - (1,2,3) +julia> b = SA.DiracBasis{UInt32}(G); julia> RG = StarAlgebra(G, b) -*-algebra of Permutation group on 2 generators of order 6 +*-algebra of PermGroup( (1,2), (1,2,3) ) ``` -This creates the group algebra of the symmetric group. How do we compute inside the group algebra? There are a few ways to comstruct elements: +This creates the group algebra of the symmetric group. There are a few ways to +comstruct its elements: ```julia julia> zero(RG) 0·() @@ -48,77 +44,86 @@ julia> RG(-5.0) # coerce a scalar to the ring -5.0·() julia> RG(rand(G)) # the indicator function on a random element of G -1·(1,3,2) +1·() -julia> f = AlgebraElement(rand(-3:3, length(b)), RG) # an element given by vectors of coefficients in the basis -1·() -1·(2,3) +3·(1,2) -1·(1,3,2) -3·(1,3) +3·(1,2,3) +julia> f = sum(rand(-3:3)*RG(rand(G)) for _ in 1:12) # an element given by vectors of coefficients in the basis +-1·() +4·(1,2) +5·(1,2,3) -2·(1,3,2) -2·(1,3) + +julia> SA.star(g::PermutationGroups.AbstractPermutation) = inv(g); star(f) # the star involution +-1·() +4·(1,2) +5·(1,3,2) -2·(1,2,3) -2·(1,3) + +julia> f' # the same +-1·() +4·(1,2) +5·(1,3,2) -2·(1,2,3) -2·(1,3) + +julia> f' == f +false + +julia> f - f' # the alternating part +7·(1,2,3) -7·(1,3,2) + +julia> StarAlgebras.aug(p::PermutationGroups.Permutation) = 1 + +julia> StarAlgebras.aug(f) # sum of coefficients +4 + +julia> StarAlgebras.aug(f - f') # sum of coefficients +0 + +julia> using LinearAlgebra; norm(f, 2) # 2-norm +25.0 ``` -One may work with such element using the following functions: +By default `f` here is stored as (sparse) coefficients against `basis(RG)`: ```julia -julia> StarAlgebras.coeffs(f) -6-element Vector{Int64}: - 1 - -1 - 3 - -1 - -3 - 3 - -julia> StarAlgebras.star(p::PermutationGroups.AbstractPerm) = inv(p); star(f) # the star involution -1·() -1·(2,3) +3·(1,2) +3·(1,3,2) -3·(1,3) -1·(1,2,3) +julia> coeffs(f) # same as SA.coeffs(f, basis(RG)) +StarAlgebras.SparseCoefficients{...}(Permutation{...}[(), (1,2), (1,2,3), (1,3,2), (1,3)], [-1, 4, 5, -2, -2]) +``` -julia> f' # the same -1·() -1·(2,3) +3·(1,2) +3·(1,3,2) -3·(1,3) -1·(1,2,3) +Working directly with elements of `G` one can also write -julia> g = rand(G); g +```julia +julia> g = Permutation(perm"(1,2,3)", G) (1,2,3) -julia> StarAlgebras.coeffs(RG(g)) # note the type of coefficients -6-element SparseArrays.SparseVector{Int64, UInt8} with 1 stored entry: - [6] = 1 - julia> x = one(RG) - 3RG(g); supp(x) # support of the funtion -2-element Vector{Permutation{...}: +2-element Vector{Permutation{…}}: () (1,2,3) julia> x(g) # value of x at g -3 -julia> x[g] += 3; x # modification of x in-place -1·() - -julia> aug(f) # sum of coefficients -2 - -julia> using LinearAlgebra; norm(f, 2) # 2-norm -5.477225575051661 +julia> x(inv(g)) +0 ``` -Using this we can define e.g. a few projections in `RG` and check their orthogonality: +## Example: Projections in group algebra + +Using these functions let's define a few projections in `RG` and check their +orthogonality: ```julia julia> using Test -julia> Base.sign(p::PermutationGroups.Permutation) = sign(p.perm); - -julia> l = length(b) +julia> l = length(basis(RG)) 6 -julia> P = sum(RG(g) for g in b) // l # projection to the subspace fixed by G -1//6·() +1//6·(2,3) +1//6·(1,2) +1//6·(1,3,2) +1//6·(1,3) +1//6·(1,2,3) +julia> P = sum(RG(g) for g in b) // l # projection to the subspace fixed by all elements of G +1//6·() +1//6·(2,3) +1//6·(1,2) +1//6·(1,2,3) +1//6·(1,3,2) +1//6·(1,3) julia> @test P * P == P Test Passed julia> P3 = 2 * sum(RG(g) for g in b if sign(g) > 0) // l # projection to the subspace fixed by Alt(3) = C₃ -1//3·() +1//3·(1,3,2) +1//3·(1,2,3) +1//3·() +1//3·(1,2,3) +1//3·(1,3,2) julia> @test P3 * P3 == P3 Test Passed -julia> P2 = (RG(1) + RG(b[2])) // 2 # projection to the C₂-fixed subspace -1//2·() +1//2·(2,3) +julia> h = PermutationGroups.gens(G,1) +(1,2) + +julia> P2 = (RG(one(G)) + RG(h)) // 2 # projection to the C₂-fixed subspace +1//2·() +1//2·(1,2) julia> @test P2 * P2 == P2 Test Passed @@ -126,20 +131,25 @@ Test Passed julia> @test P2 * P3 == P3 * P2 == P # their intersection is precisely the same as the one for G Test Passed -julia> P2m = (RG(1) - RG(b[2])) // 2 # orthogonal C₂-fixed subspace -1//2·() -1//2·(2,3) +julia> P2m = (RG(1) - RG(h)) // 2 # projection onto the subspace orthogonal to C₂-fixed subspace +1//2·() -1//2·(1,2) julia> @test P2m * P2m == P2m Test Passed julia> @test iszero(P2m * P2) # indeed P2 and P2m are orthogonal Test Passed - ``` -This package originated as a tool to compute sum of hermitian squares in `*`-algebras. These consist not of standard `f*f` summands, but rather `star(f)*f`. You may think of semi-definite matrices: their Cholesky decomposition determines `P = Q'·Q`, where `Q'` denotes transpose. Algebra of matrices with transpose is an (the?) example of a `*`-algebra. To compute such sums of squares one may either sprinkle the code with `star`s, or `'` (aka `Base.adjoint` postfix symbol): +This package originated as a tool to compute sum of hermitian squares in +`*`-algebras. These consist not of standard `f*f` summands, but rather +`star(f)*f`. You may think of semi-definite matrices: their Cholesky +decomposition determines `P = Q'·Q`, where `Q'` denotes (conjugate) transpose. +Algebra of matrices with transpose is an (the?) example of a `*`-algebra. +To compute such sums of squares one may either sprinkle the code with `star`s, +or `'` (aka `Base.adjoint` postfix symbol): ```julia -julia> x = RG(G(perm"(1,2,3)")) +julia> x = RG(Permutation(perm"(1,2,3)", G)) 1·(1,2,3) julia> X = one(RG) - x @@ -149,99 +159,77 @@ julia> X' 1·() -1·(1,3,2) julia> X'*X -2·() -1·(1,3,2) -1·(1,2,3) +2·() -1·(1,2,3) -1·(1,3,2) julia> @test X'*X == star(X)*X == 2one(X) - x - star(x) Test Passed ``` -### More advanced use +## Fixed basis and translation between bases -`RG = StarAlgebra(G, b)` creates the algebra with `TrivialMStructure`, i.e. a multiplicative structure which computes product of basis elements every time it needs it. This of course may be wastefull, e.g. the computed products could be stored in a cache for future use. There are two options here: -```julia -julia> mt = StarAlgebras.MTable(b, table_size=(length(b), length(b))) -6×6 StarAlgebras.MTable{UInt8, false, Matrix{UInt8}}: - 0x01 0x02 0x03 0x04 0x05 0x06 - 0x02 0x01 0x06 0x05 0x04 0x03 - 0x03 0x04 0x01 0x02 0x06 0x05 - 0x04 0x03 0x05 0x06 0x02 0x01 - 0x05 0x06 0x04 0x03 0x01 0x02 - 0x06 0x05 0x02 0x01 0x03 0x04 +`b = SA.DiracBasis{UInt32}(G)` takes an iterator (in this case a finite +permutation group `G`) and creates a basis with Dirac multiplicative structure. +This means a basis of a linear space with the multiplicative structure where +multiplication of two basis elements results in a third one. This computation +does not depend on the finiteness of `G`, i.e. is fully lazy. -``` -creates an eagerly computed multiplication table on elements of `b`. Keyword `table_size` is used to specify the table size (above: it's the whole multiplication table). Since `MTable<:AbstractMatrix`, one can use the indexing syntax `mt[i,j]` to compute the **index** of the product of `i`-th and `j`-th elements of the basis. For example -```julia -julia> g = G(perm"(1,2,3)"); h = G(perm"(2,3)"); +When the basis is known a priori one can create an efficient `SA.FixedBasis` +which caches (memoizes) the results of multiplications. For example: -julia> i, j = b[g], b[h] # indices of g and h in basis b -(0x06, 0x02) - -julia> k = mt[i,j] # the index of the product -0x05 +```julia +julia> fb = let l = length(basis(RG)) + StarAlgebras.FixedBasis(b; n=l, mt=l) # mt can be also smaller than l +end; -julia> @test b[k] == g*h -Test Passed +julia> fRG = StarAlgebra(G, fb) +*-algebra of PermGroup( (1,2), (1,2,3) ) ``` -The second option is -```julia -julia> cmt = StarAlgebras.CachedMTable(b, table_size=(length(b), length(b))); -``` -This multiplication table is lazy, i.e. products will be computed and stored only when actually needed. Additionally, one may call +Since the length of the basis is known this time, algebra elements can be stored simply as (sparse) vectors: + ```julia -julia> using SparseArrays +julia> coeffs(one(fRG)) +6-element SparseArrays.SparseVector{Int64, Int64} with 1 stored entry: + [1] = 1 -julia> StarAlgebras.CachedMTable(b, spzeros(UInt8, length(b), length(b))); +julia> AlgebraElement(ones(Int, length(basis(fRG)))//6, fRG) +1//6·() +1//6·(2,3) +1//6·(1,2) +1//6·(1,2,3) +1//6·(1,3,2) +1//6·(1,3) ``` -to specify storage type of the matrix (by default it's a simple dense `Matrix`). -This may be advisable when a few products are computed repeatedly on a quite large basis. -```julia -julia> RGc = StarAlgebra(G, b, cmt) -*-algebra of Permutation group on 2 generators of order 6 +To translate coefficients between bases one may call -``` -should be functinally equivalent to `RG` above, however it will cache computation of products lazily. A word of caution is needed here though. Even though `RGc` and `RG` are functionally equivalent, they are not **comparable** in the sense that e.g. ```julia -julia> @test one(RGc) != one(RG) -Test Passed +julia> coeffs(P2, fb) +6-element SparseArrays.SparseVector{Rational{Int64}, Int64} with 2 stored entries: + [1] = 1//2 + [3] = 1//2 -``` - -This is a conscious decision on our part, as comparing algebraic structures is easier said than done ;) To avoid solving this conundrum (are bases equal? are multiplicative structures equal? are these permuted by a compatible permutation? or maybe a linear transformation was applied to the basis, resulting in a different, but equivalent multiplicative structure?), elements could be mixed together **only if their parents are identically** (i.e. `===`) **equal**. - -Finally, if the group is infinite (or just too large), but we need specific products, we may reduce the table_size to the required size (it doesn't have to be `length(b) × length(b)`). Note that in such case asking for a product outside of multiplication table will rise `ProductNotDefined` exception. +julia> coeffs(coeffs(P2), b, fb) # same as above, from source basis to target basis +6-element SparseArrays.SparseVector{Rational{Int64}, Int64} with 2 stored entries: + [1] = 1//2 + [3] = 1//2 -### Even more advanced use (for experts only) +``` -For low-level usage `MultiplicativeStructures` follow the sign convention: -```julia -julia> mt = StarAlgebras.CachedMTable(b, table_size=(length(b), length(b))); +Translation in the opposite direction is also possible ```julia -julia> k = mt[-i,j] -0x06 +julia> fP2 = AlgebraElement(StarAlgebras.coeffs(P2,fb), fRG) +1//2·() +1//2·(1,2) -julia> @test star(b[i])*b[j] == b[k] -Test Passed +julia> StarAlgebras.coeffs(fP2, b) +StarAlgebras.SparseCoefficients{…}(Permutation{…}[(), (1,2)], Rational{Int64}[1//2, 1//2]) -``` -Note that this (minus-twisted) "product" is no longer associative! Observe: -```julia -julia> @test mt[mt[3, 5], 4] == mt[3, mt[5, 4]] # (b[3]*b[4])*b[5] == b[3]*(b[4]*b[5]) -Test Passed +julia> P2_ = AlgebraElement(StarAlgebras.coeffs(fP2, b), RG) +1//2·() +1//2·(1,2) -julia> @test mt[-signed(mt[-3, 5]), 4] == 0x06 # star(star(b[3])*b[5])*b[4] = star(b[5])*b[3]*b[4] +julia> @test P2 == P2_ Test Passed -julia> @test mt[-3, mt[-5, 4]] == 0x01 # star(b[3])*star(b[5])*b[4] -Test Passed ``` - - ----- If you happen to use this package please cite either [1712.07167](https://arxiv.org/abs/1712.07167) or [1812.03456](https://arxiv.org/abs/1812.03456). This package superseeds [GroupRings.jl](https://github.com/kalmarek/GroupRings.jl) which was developed and used there. It served its purpose well. Let it rest peacefully. diff --git a/src/StarAlgebras.jl b/src/StarAlgebras.jl index aaeadad..7540279 100644 --- a/src/StarAlgebras.jl +++ b/src/StarAlgebras.jl @@ -3,17 +3,37 @@ module StarAlgebras using SparseArrays import LinearAlgebra +import MutableArithmetics as MA + export StarAlgebra, AlgebraElement -export aug, basis, coeffs, star, supp +export basis, coeffs, star -include("bases.jl") +# AbstractCoefficients +## abstract definitions +include("coefficients.jl") +## concrete implementation +include("sparse_coeffs.jl") +# MultiplicativeStructures include("mstructures.jl") include("mtables.jl") +# AbstractBases +## abstract definitions +include("bases.jl") +# concrete implementations +include("bases_dirac.jl") +include("bases_fixed.jl") + +# Algebras and elts include("types.jl") include("algebra_elts.jl") +include("star.jl") + include("arithmetic.jl") include("show.jl") +# augmented basis implementation +include("diracs_augmented.jl") + end diff --git a/src/algebra_elts.jl b/src/algebra_elts.jl index 467b470..e2e160b 100644 --- a/src/algebra_elts.jl +++ b/src/algebra_elts.jl @@ -5,66 +5,37 @@ function Base.:(==)(X::AlgebraElement, Y::AlgebraElement) return coeffs(X) == coeffs(Y) end -Base.getindex(a::AlgebraElement, i::Integer) = coeffs(a)[i] -Base.getindex(a::AlgebraElement{<:StarAlgebra{O,T}}, x::T) where {O,T} = - (b = basis(parent(a)); a[b[x]]) +Base.copy(a::AlgebraElement) = AlgebraElement(copy(coeffs(a)), parent(a)) + +function Base.deepcopy_internal(a::AlgebraElement, id::IdDict) + if !haskey(id, a) + id[a] = AlgebraElement(Base.deepcopy_internal(coeffs(a), id), parent(a)) + end + return id[a] +end # call overload: -(a::AlgebraElement)(x) = a[x] +(a::AlgebraElement)(x) = coeffs(a)[basis(parent(a))[x]] +Base.setindex!(a::AlgebraElement, v, idx) = a.coeffs[basis(parent(a))[idx]] = v -Base.setindex!(a::AlgebraElement, v, i::Integer) = a.coeffs[i] = v +# AlgebraElement specific functions -function Base.setindex!( - a::AlgebraElement{<:StarAlgebra{O,T}}, - v, - t::T, -) where {O,T} +function supp(a::AlgebraElement) b = basis(parent(a)) - return a[b[t]] = v + return [b[i] for (i, _) in nonzero_pairs(coeffs(a))] end -# AlgebraElement specific functions - -supp_ind(a::AlgebraElement) = findall(!iszero, coeffs(a)) -supp_ind(a::AlgebraElement{A,T,<:SparseVector}) where {A,T} = - (dropzeros!(coeffs(a)); SparseArrays.nonzeroinds(coeffs(a))) -supp(a::AlgebraElement) = (b = basis(parent(a)); [b[i] for i in supp_ind(a)]) - -function star(A::StarAlgebra, i::Integer) - @assert i > 0 - if i < max(size(A.mstructure)...) - return _get(A.mstructure, -signed(i)) - else - b = basis(A) - return b[star(b[i])] - end +function LinearAlgebra.norm(a::AlgebraElement, p::Real) + return LinearAlgebra.norm(coeffs(a), p) end -function star(X::AlgebraElement) - A = parent(X) - b = basis(A) - supp_X = supp_ind(X) - idcs = similar(supp_X) - vals = similar(idcs, eltype(X)) - for (i, idx) in enumerate(supp_X) - idcs[i] = star(parent(X), idx) - vals[i] = X[idx] - end - return AlgebraElement(sparsevec(idcs, vals, length(b)), A) +function LinearAlgebra.dot(a::AlgebraElement, b::AlgebraElement) + return LinearAlgebra.dot(coeffs(a), coeffs(b)) end -Base.adjoint(a::AlgebraElement) = star(a) - -LinearAlgebra.norm(a::AlgebraElement, p::Real) = - LinearAlgebra.norm(coeffs(a), p) -aug(a::AlgebraElement) = sum(coeffs(a)) - -LinearAlgebra.dot(a::AlgebraElement, v::AbstractVector) = - LinearAlgebra.dot(StarAlgebras.coeffs(a), v) -LinearAlgebra.dot(v::AbstractVector, a::AlgebraElement) = LinearAlgebra.dot(a, v) - -Base.copy(a::AlgebraElement) = AlgebraElement(copy(coeffs(a)), parent(a)) -function Base.deepcopy_internal(a::AlgebraElement, stackdict::IdDict) - haskey(stackdict, a) && return stackdict[a] - return AlgebraElement(Base.deepcopy_internal(coeffs(a), stackdict), parent(a)) +function LinearAlgebra.dot(w::AbstractVector, b::AlgebraElement) + return LinearAlgebra.dot(w, coeffs(b)) +end +function LinearAlgebra.dot(a::AlgebraElement, w::AbstractVector) + return LinearAlgebra.dot(coeffs(a), w) end diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 40c1450..0ae8e1b 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -1,149 +1,109 @@ -# module structure: - -function Base.:*(a::Number, X::AlgebraElement) - T = Base._return_type(*, Tuple{eltype(X),typeof(a)}) - return mul!(similar(X, T), X, a) +function _preallocate_output(X::AlgebraElement, a::Number, op) + T = MA.promote_operation(op, eltype(X), typeof(a)) + return similar(X, T) end -Base.:*(X::AlgebraElement, a::Number) = a * X -Base.:(/)(X::AlgebraElement, a::Number) = inv(a) * X - -# TODO: handle this through mul!? -Base.:(//)(X::AlgebraElement, a::Number) = AlgebraElement(coeffs(X) .// a, parent(X)) -Base.:div(X::AlgebraElement, a::Number) = AlgebraElement(div.(coeffs(X), a), parent(X)) - -# ring structure: -Base.:-(X::AlgebraElement) = neg!(similar(X), X) - function _preallocate_output(X::AlgebraElement, Y::AlgebraElement, op) - T = Base._return_type(op, Tuple{eltype(X),eltype(Y)}) - if coeffs(Y) isa DenseArray + T = MA.promote_operation(op, eltype(X), eltype(Y)) + if coeffs(Y) isa DenseArray # what a hack :) return similar(Y, T) end return similar(X, T) end -Base.:+(X::AlgebraElement, Y::AlgebraElement) = add!(_preallocate_output(X, Y, +), X, Y) -Base.:-(X::AlgebraElement, Y::AlgebraElement) = sub!(_preallocate_output(X, Y, -), X, Y) -Base.:*(X::AlgebraElement, Y::AlgebraElement) = mul!(_preallocate_output(X, Y, *), X, Y) +# module structure: + +Base.:*(X::AlgebraElement, a::Number) = a * X +Base.:(/)(X::AlgebraElement, a::Number) = inv(a) * X +Base.:(//)(X::AlgebraElement, a::Number) = 1 // a * X +function Base.:-(X::AlgebraElement) + return MA.operate_to!(_preallocate_output(X, -1, *), -, X) +end +function Base.:*(a::Number, X::AlgebraElement) + return MA.operate_to!(_preallocate_output(X, a, *), *, X, a) +end +function Base.:div(X::AlgebraElement, a::Number) + return MA.operate_to!(_preallocate_output(X, a, div), div, X, a) +end + +function Base.:+(X::AlgebraElement, Y::AlgebraElement) + return MA.operate_to!(_preallocate_output(X, Y, +), +, X, Y) +end +function Base.:-(X::AlgebraElement, Y::AlgebraElement) + return MA.operate_to!(_preallocate_output(X, Y, -), -, X, Y) +end +function Base.:*(X::AlgebraElement, Y::AlgebraElement) + return MA.operate_to!(_preallocate_output(X, Y, *), *, X, Y) +end Base.:^(a::AlgebraElement, p::Integer) = Base.power_by_squaring(a, p) -# mutable API; TODO: replace with MutableArithmetic +# mutable API -function zero!(a::AlgebraElement) - a.coeffs .= zero(eltype(coeffs(a))) +function MA.operate!(::typeof(zero), a::AlgebraElement) + MA.operate!(zero, coeffs(a)) return a end -function neg!(res::AlgebraElement, X::AlgebraElement) +function MA.operate_to!( + res::AlgebraElement, + ::typeof(*), + X::AlgebraElement, + a::Number, +) @assert parent(res) === parent(X) - res.coeffs .= .-coeffs(X) + MA.operate_to!(coeffs(res), *, coeffs(X), a) return res end -function add!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement) +function MA.operate_to!( + res::AlgebraElement, + ::typeof(div), + X::AlgebraElement, + a::Number, +) @assert parent(res) === parent(X) - @assert parent(X) === parent(Y) - if res === X - for (idx, y) in _nzpairs(coeffs(Y)) - res[idx] += y - end - elseif res === Y - for (idx, x) in _nzpairs(coeffs(X)) - res[idx] += x - end - else - zero!(res) - for (idx, x) in _nzpairs(coeffs(X)) - res[idx] += x - end - for (idx, y) in _nzpairs(coeffs(Y)) - res[idx] += y - end - end + MA.operate_to!(coeffs(res), div, coeffs(X), a) return res end -function sub!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement) - @assert parent(res) === parent(X) === parent(Y) - neg!(res, Y) - add!(res, res, X) - return res -end - -function mul!(res::AlgebraElement, X::AlgebraElement, a::Number) +function MA.operate_to!(res::AlgebraElement, ::typeof(-), X::AlgebraElement) @assert parent(res) === parent(X) - if res !== X - zero!(res) - end - for (idx, x) in _nzpairs(coeffs(X)) - res.coeffs[idx] = a * x - end + MA.operate_to!(coeffs(res), -, coeffs(X)) return res end -function mul!( - res::AbstractVector, - X::AbstractVector, - Y::AbstractVector, - ms::MultiplicativeStructure, +function MA.operate_to!( + res::AlgebraElement, + ::typeof(+), + X::AlgebraElement, + Y::AlgebraElement, ) - res = (res === X || res === Y) ? zero(res) : (res .= zero(eltype(res))) - return fmac!(res, X, Y, ms) -end - -function mul!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement) - res = (res === X || res === Y) ? zero(res) : zero!(res) - return fmac!(res, X, Y) -end - -function fmac!(res::AlgebraElement, X::AlgebraElement, Y::AlgebraElement) - @assert parent(res) === parent(X) === parent(Y) - fmac!(coeffs(res), coeffs(X), coeffs(Y), parent(res).mstructure) + @assert parent(res) === parent(X) + @assert parent(X) === parent(Y) + MA.operate_to!(coeffs(res), +, coeffs(X), coeffs(Y)) return res end -_nzpairs(v::AbstractVector) = pairs(v) -_nzpairs(v::AbstractSparseVector) = - zip(SparseArrays.nonzeroinds(v), SparseArrays.nonzeros(v)) - -function fmac!( - res::AbstractVector, - X::AbstractVector, - Y::AbstractVector, - mstr::MultiplicativeStructure, +function MA.operate_to!( + res::AlgebraElement, + ::typeof(-), + X::AlgebraElement, + Y::AlgebraElement, ) - @assert res !== X - @assert res !== Y - lx, ly = size(mstr) - @assert all(iszero, @view(X[lx+1:end])) - @assert all(iszero, @view(Y[ly+1:end])) - for iy in 1:ly - y = Y[iy] - iszero(y) && continue - for ix in 1:lx - x = X[ix] - iszero(x) && continue - res[mstr[ix, iy]] += X[ix] * y - end - end + @assert parent(res) === parent(X) === parent(Y) + MA.operate_to!(coeffs(res), -, coeffs(X), coeffs(Y)) return res end -function fmac!( - res::AbstractVector, - X::AbstractSparseVector, - Y::AbstractSparseVector, - mstr::MultiplicativeStructure, +function MA.operate_to!( + res::AlgebraElement, + ::typeof(*), + X::AlgebraElement, + Y::AlgebraElement, ) - @assert res !== X - @assert res !== Y - for iy in SparseArrays.nonzeroinds(Y) - y = Y[iy] - for ix in SparseArrays.nonzeroinds(X) - res[mstr[ix, iy]] += X[ix] * y - end - end + @assert parent(res) === parent(X) === parent(Y) + mstr = mstructure(basis(parent(res))) + MA.operate_to!(coeffs(res), mstr, coeffs(X), coeffs(Y)) return res end diff --git a/src/bases.jl b/src/bases.jl index b068f85..afd5951 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -1,26 +1,79 @@ -abstract type AbstractBasis{T,I} <: AbstractVector{T} end +""" + AbstractBasis{T,I} +Implements -struct Basis{T,I,A<:AbstractVector{T}} <: AbstractBasis{T,I} - basis::A - rbasis::Dict{T,I} -end + * `in(x, basis)` for basis elements + * `Base.getindex(b::AbstractBasis{T,I}, idx::I)::T` + * iteration protocol + +`AbstractBasis` from the outside may look like a vector with arbitrary indices, +however it may be possibly infinite (and therefore does not define its length). + +When `I<:Integer` one may store coefficients w.r.t. the basis in an ordinary +vector, however it's length has no relation to the size of the basis (which may +still be infinite). + +* [`ImplicitBasis`](@ref) is best suited for bases which can be enumerated but + are in princle infinite in size (e.g. monomial basis for the polynomial ring). + An example implementation of [`DiracBasis`](@ref) wrapping an iterator can be + used for this purpose. +* [`ExplicitBasis`](@ref) can be used for fixed, finite size bases, in particular + [`FixedBasis`](@ref) implements `AbstractVector`-type storage for elements. + This can be used e.g. for representing polynomials as (sparse) vectors of + coefficients w.r.t. a given `FixedBasis`. +expressed in it. +""" +abstract type AbstractBasis{T,I} end + +Base.eltype(::Type{<:AbstractBasis{T}}) where {T} = T +Base.eltype(b::AbstractBasis) = eltype(typeof(b)) +key_type(::Type{<:AbstractBasis{T,I}}) where {T,I} = I +key_type(b::AbstractBasis) = key_type(typeof(b)) + +""" + ImplicitBasis{T,I} +Implicit bases are not stored in memory and can be potentially infinite. +""" +abstract type ImplicitBasis{T,I} <: AbstractBasis{T,I} end -function Basis{I}(basis::AbstractVector) where {I} - length(basis) <= typemax(I) || - throw("index type $I is to small for basis of length $(length(basis))") - @assert !(eltype(basis) <: Integer) - return Basis(basis, Dict(b => I(idx) for (idx, b) in pairs(basis))) +function zero_coeffs(::Type{S}, ::ImplicitBasis{T}) where {S,T} + return SparseCoefficients(T[], S[]) end -Base.size(b::Basis) = size(b.basis) +""" + ExplicitBasis{T,I} +Explicit bases are stored e.g. in an `AbstractVector` and hence immutable +(in particular of well defined and fixed length). +""" +abstract type ExplicitBasis{T,I} <: AbstractBasis{T,I} end -Base.IndexStyle(::Type{<:Basis{T,I,A}}) where {T,I,A} = Base.IndexStyle(A) +Base.keys(ib::ExplicitBasis) = Base.OneTo(length(ib)) -Base.@propagate_inbounds Base.getindex(b::Basis, i::Integer) = b.basis[i] -Base.@propagate_inbounds Base.getindex(b::Basis{T}, g::T) where {T} = b.rbasis[g] +function zero_coeffs(::Type{S}, eb::ExplicitBasis{T,I}) where {S,T,I} + return spzeros(S, I, length(eb)) +end -Base.in(g, b::Basis) = haskey(b.rbasis, g) +function Base.getindex(eb::ExplicitBasis, range::AbstractRange{<:Integer}) + return [eb[i] for i in range] +end -## are these really that useful? -SparseArrays.indtype(::Type{<:Basis{T,I}}) where {T,I} = I -SparseArrays.indtype(b::Basis) = SparseArrays.indtype(typeof(b)) +""" + coeffs(cfs, source, target) +Translate coefficients `cfs` in `source::AbstractBasis` to basis +`target::AbstractBasis`. +""" +function coeffs(cfs, source::AbstractBasis, target::AbstractBasis) + source === target && return cfs + source == target && return cfs + res = zero_coeffs(valtype(cfs), target) + return coeffs!(res, cfs, source, target) +end + +function coeffs!(res, cfs, source::AbstractBasis, target::AbstractBasis) + MA.operate!(zero, res) + for (k, v) in nonzero_pairs(cfs) + x = source[k] + res[target[x]] += v + end + return res +end diff --git a/src/bases_dirac.jl b/src/bases_dirac.jl new file mode 100644 index 0000000..e4135ef --- /dev/null +++ b/src/bases_dirac.jl @@ -0,0 +1,30 @@ +mutable struct DiracBasis{T,I,S,M<:DiracMStructure} <: ImplicitBasis{T,I} + object::S # any iterable + moperation::M + + function DiracBasis{I}(itr, operation = *) where {I} + @assert !isempty(itr) + mstr = DiracMStructure(operation) + return new{eltype(itr),I,typeof(itr),typeof(mstr)}(itr, mstr) + end +end + +object(db::DiracBasis) = db.object +mstructure(db::DiracBasis{T}) where {T} = db.moperation + +function Base.IteratorSize(::Type{<:DiracBasis{T,I,S}}) where {T,I,S} + return Base.IteratorSize(S) +end +function Base.length(db::DiracBasis) + @assert Base.haslength(object(db)) + return length(object(db)) +end +Base.iterate(db::DiracBasis) = iterate(object(db)) +Base.iterate(db::DiracBasis, st) = iterate(object(db), st) + +Base.in(g, db::DiracBasis) = g in object(db) + +function Base.getindex(db::DiracBasis{T}, x::T) where {T} + @assert x in object(db) + return x +end diff --git a/src/bases_fixed.jl b/src/bases_fixed.jl new file mode 100644 index 0000000..ce302b7 --- /dev/null +++ b/src/bases_fixed.jl @@ -0,0 +1,41 @@ +mutable struct FixedBasis{T,I,V<:AbstractVector{T},M<:MTable{T,I}} <: + ExplicitBasis{T,I} + elts::V + table::M +end + +function FixedBasis(basis::AbstractBasis; n::Integer, mt::Integer = UInt32(0)) + @assert 0 ≤ mt ≤ n + elts = Iterators.take(basis, n) + return FixedBasis(collect(elts), mstructure(basis), (mt, mt)) +end + +function FixedBasis( + elts::AbstractVector, + mstr::MultiplicativeStructure, + dims::NTuple{2,I} = (UInt32(0), UInt32(0)), +) where {I<:Integer} + @assert 0 ≤ dims[1] ≤ length(elts) + @assert 0 ≤ dims[2] ≤ length(elts) + @assert !(eltype(elts) <: Integer) + return FixedBasis(elts, MTable(elts, mstr, dims)) +end + +mstructure(fb::FixedBasis) = fb.table + +Base.in(x, b::FixedBasis) = haskey(mstructure(b), x) +Base.getindex(b::FixedBasis{T}, x::T) where {T} = mstructure(b)[x] +Base.getindex(b::FixedBasis, i::Integer) = mstructure(b)[i] + +Base.IteratorSize(::Type{<:FixedBasis}) = Base.HasLength() +Base.length(b::FixedBasis) = length(b.elts) + +Base.iterate(b::FixedBasis) = iterate(b.elts) +Base.iterate(b::FixedBasis, state) = iterate(b.elts, state) +Base.IndexStyle(::Type{<:FixedBasis{T,I,V}}) where {T,I,V} = Base.IndexStyle(V) + +# To break ambiguity +Base.@propagate_inbounds Base.getindex( + b::FixedBasis{T,I}, + i::I, +) where {T,I<:Integer} = b.elts[i] diff --git a/src/coefficients.jl b/src/coefficients.jl new file mode 100644 index 0000000..02f40db --- /dev/null +++ b/src/coefficients.jl @@ -0,0 +1,243 @@ +""" + abstract type AbstractCoefficients{V,K} end +Everything that implements a fixed set of methods can be used as +`SparseCoefficients` without subtyping it. + +# "Read-only" coefficients +E.g. returned by calls to a `MultiplicativeStructure` need to implement + +* `Base.keys` +* `Base.values` +* `StarAlgebras.canonical` +* `StarAlgebras.star(b::AbstractBasis, ac)` + +For general types either all necessary arithmetic operations need to be +implemented, or fallbacks using the framework of `MutableArithmetics` are +provided based on random indexing. Additionally one needs to provide: + +* `Base.similar(ac, T::Type)` with the same semantics as the one for vectors +* `Base.getindex(ac, idx)` +* `Base.setindex!(ac, val, idx)` +* `MutableArithmetics.operate!(ms::UnsafeAddMul, ac, v::C, w::C) where C<:SA.AbstractCoefficients` + +""" +abstract type AbstractCoefficients{K,V} end + +key_type(::Type{<:AbstractCoefficients{K}}) where {K} = K +Base.valtype(::Type{<:AbstractCoefficients{K,V}}) where {K,V} = V +key_type(b::AbstractCoefficients) = key_type(typeof(b)) +Base.valtype(b::AbstractCoefficients) = valtype(typeof(b)) + +key_type(b) = keytype(b) +# `keytype(::Type{SparseVector{V,K}})` is not defined so it falls +# back to `keytype{::Type{<:AbstractArray})` which returns `Int`. +key_type(::Type{SparseArrays.SparseVector{V,K}}) where {V,K} = K +key_type(v::SparseArrays.SparseVector) = key_type(typeof(v)) + +Base.iszero(ac::AbstractCoefficients) = isempty(keys(ac)) + +Base.similar(ac::AbstractCoefficients) = similar(ac, valtype(ac)) + +""" + canonical(ac::AbstractCoefficients) +Compute the canonical form of `ac` (e.g. grouping coefficients together, etc). + +If `ac` can be brought to canonical form in-place one has to implement +* `MA.mutability(::Type{typeof(ac)}, canonical, ::Vararg{Type}) = MA.IsMutable()` +* `MA.operate!(canonical, ac)` that performs this canonicalization. + +otherwise `canonical(ac)` needs to be implemented. +""" +function canonical end +function MA.promote_operation(::typeof(canonical), ::Type{C}) where {C} + return C +end + +# example implementation for vectors +MA.operate!(::typeof(canonical), sv::SparseVector) = dropzeros!(sv) +MA.operate!(::typeof(canonical), v::Vector) = v + +function Base.:(==)(ac1::AbstractCoefficients, ac2::AbstractCoefficients) + MA.operate!(canonical, ac1) + MA.operate!(canonical, ac2) + all(x -> ==(x...), zip(keys(ac1), keys(ac2))) || return false + all(x -> ==(x...), zip(values(ac1), values(ac2))) || return false + return true +end + +function Base.hash(ac::AbstractCoefficients, h::UInt) + MA.operate!(canonical, ac) + return foldl((h, i) -> hash(i, h), nonzero_pairs(ac); init = h) +end + +""" + nonzero_pairs(ac::AbstractCoefficients) +Return an iterator over pairs `(k=>v)` of keys and values stored in `ac`. + +The iterator contains all pairs with `v` potentially non-zero. +""" +@inline function nonzero_pairs(ac) + return (k => v for (k, v) in zip(keys(ac), values(ac))) +end + +@inline nonzero_pairs(v::AbstractVector) = + (p for p in pairs(v) if !iszero(last(p))) +@inline function nonzero_pairs(v::AbstractSparseVector) + return zip(SparseArrays.nonzeroinds(v), SparseArrays.nonzeros(v)) +end + +function LinearAlgebra.norm(sc::AbstractCoefficients, p::Real) + isempty(values(sc)) && return (0^p)^(1 / p) + return sum(abs(v)^p for v in values(sc))^(1 / p) +end + +function LinearAlgebra.dot(ac::AbstractCoefficients, bc::AbstractCoefficients) + if isempty(values(ac)) || isempty(values(bc)) + return zero(MA.promote_sum_mul(valtype(ac), valtype(bc))) + else + return sum(c * star(bc[i]) for (i, c) in nonzero_pairs(ac)) + end +end + +function LinearAlgebra.dot(w::AbstractVector, ac::AbstractCoefficients) + @assert key_type(ac) <: Integer + if isempty(values(ac)) + return zero(MA.promote_sum_mul(eltype(w), valtype(ac))) + else + return sum(w[i] * star(v) for (i, v) in nonzero_pairs(ac)) + end +end + +function LinearAlgebra.dot(ac::AbstractCoefficients, w::AbstractVector) + @assert key_type(ac) <: Integer + if isempty(values(ac)) + return zero(MA.promote_sum_mul(eltype(w), valtype(ac))) + else + return sum(v * star(w[i]) for (i, v) in nonzero_pairs(ac)) + end +end + +# general mutable API +# why here? +MA.operate!(::typeof(zero), v::SparseVector) = (v .= 0; v) + +Base.zero(X::AbstractCoefficients) = MA.operate!(zero, similar(X)) +Base.:-(X::AbstractCoefficients) = MA.operate_to!(__prealloc(X, -1, *), -, X) +Base.:*(a::Number, X::AbstractCoefficients) = X * a +Base.:/(X::AbstractCoefficients, a::Number) = X * inv(a) +Base.://(X::AbstractCoefficients, a::Number) = X * 1 // a + +function Base.:*(X::AbstractCoefficients, a::Number) + return MA.operate_to!(__prealloc(X, a, *), *, X, a) +end +function Base.:div(X::AbstractCoefficients, a::Number) + return MA.operate_to!(__prealloc(X, a, div), div, X, a) +end +function Base.:+(X::AbstractCoefficients, Y::AbstractCoefficients) + return MA.operate_to!(__prealloc(X, Y, +), +, X, Y) +end +function Base.:-(X::AbstractCoefficients, Y::AbstractCoefficients) + return MA.operate_to!(__prealloc(X, Y, -), -, X, Y) +end + +## fallbacks; require Base.setindex! implemented for AbstractCoefficients + +function MA.operate!(::typeof(zero), X::AbstractCoefficients) + for (idx, x) in nonzero_pairs(X) + X[idx] = zero(x) + end + return X +end + +# unary - +function MA.operate_to!( + res::AbstractCoefficients, + ::typeof(-), + X::AbstractCoefficients, +) + if res !== X + MA.operate!(zero, res) + end + for (idx, x) in nonzero_pairs(X) + res[idx] = -x + end + return res +end + +# scalar +function MA.operate_to!( + res::AbstractCoefficients, + ::typeof(*), + X::AbstractCoefficients, + a::Number, +) + if res !== X + MA.operate!(zero, res) + end + for (idx, x) in nonzero_pairs(X) + res[idx] = a * x + end + return res +end + +function MA.operate_to!( + res::AbstractCoefficients, + ::typeof(div), + X::AbstractCoefficients, + a::Number, +) + if res !== X + MA.operate!(zero, res) + end + for (idx, x) in nonzero_pairs(X) + res[idx] = div(x, a) + end + return res +end + +# binary ops +function MA.operate_to!( + res::AbstractCoefficients, + ::typeof(+), + X::AbstractCoefficients, + Y::AbstractCoefficients, +) + if res === X + for (idx, y) in nonzero_pairs(Y) + res[idx] += y + end + elseif res === Y + for (idx, x) in nonzero_pairs(X) + res[idx] += x + end + else + MA.operate!(zero, res) + for (idx, x) in nonzero_pairs(X) + res[idx] += x + end + for (idx, y) in nonzero_pairs(Y) + res[idx] += y + end + end + return res +end + +function MA.operate_to!( + res::AbstractCoefficients, + ::typeof(-), + X::AbstractCoefficients, + Y::AbstractCoefficients, +) + if res === X + for (idx, y) in nonzero_pairs(Y) + X[idx] -= y + end + else + if res !== Y + MA.operate!(zero, res) + end + MA.operate_to!(res, -, Y) + MA.operate_to!(res, +, res, X) + end + return res +end diff --git a/src/diracs_augmented.jl b/src/diracs_augmented.jl new file mode 100644 index 0000000..e555472 --- /dev/null +++ b/src/diracs_augmented.jl @@ -0,0 +1,134 @@ +aug(cfs::Any) = sum(values(cfs)) +aug(a::AlgebraElement) = aug(coeffs(a)) + +function aug(ac::AbstractCoefficients) + isempty(keys(ac)) && return zero(valtype(ac)) + return sum(c * aug(x) for (x, c) in nonzero_pairs(ac)) +end + +struct Augmented{K} <: AbstractCoefficients{K,Int} + elt::K +end # corresponds to (elt - 1) + +function Base.getindex(aδ::Augmented{K}, i::K) where {K} + w = aug(aδ.elt) + isone(aδ.elt) && return zero(w) + i == aδ.elt && return w + isone(i) && return -w + return zero(w) +end + +MA.operate!(::typeof(canonical), aδ::Augmented) = aδ + +Base.keys(aδ::Augmented) = (k = keys(aδ.elt); (one(first(k)), first(k))) +function Base.values(aδ::Augmented) + (e, x) = keys(aδ) + return (aδ[e], aδ[x]) +end + +function star(basis::AbstractBasis, ad::Augmented) + return Augmented(star(basis, ad.elt)) +end + +function Base.show(io::IO, aδ::Augmented) + ioc = IOContext(io, :limit => true) + if iszero(aδ) + print(ioc, "(0)") + else + e, x = keys(aδ) + _, v = values(aδ) + print(ioc, '(', -v, '·', e, '+', v, '·', x, ')') + end +end + +Base.isless(ad1::Augmented, ad2::Augmented) = isless(ad1.elt, ad2.elt) + +aug(::Augmented) = 0 + +struct AugmentedBasis{T,I,A<:Augmented{T},B<:AbstractBasis{T,I}} <: + ImplicitBasis{A,I} + basis::B +end + +function AugmentedBasis(basis::DiracBasis{T,I}) where {T,I} + @assert one(object(basis)) in basis + return AugmentedBasis{T,I,Augmented{T},typeof(basis)}(basis) +end + +object(ab::AugmentedBasis) = object(ab.basis) + +function Base.IteratorSize(::Type{<:AugmentedBasis{T,A,I,B}}) where {T,A,I,B} + return Base.IteratorSize(B) +end +Base.haslength(ab::AugmentedBasis) = Base.haslength(ab.basis) + +function Base.length(ab::AugmentedBasis) + @assert Base.haslength(ab.basis) + return length(ab.basis) - 1 +end + +function Base.iterate(ab::AugmentedBasis) + (v, st) = iterate(object(ab)) + isone(v) && return iterate(ab, st) + return Augmented(v), st +end + +function Base.iterate(ab::AugmentedBasis, st) + (v, st) = let k = iterate(object(ab), st) + isnothing(k) && return nothing + k + end + isone(v) && return iterate(ab, st) + return Augmented(v), st +end + +Base.in(g, ab::AugmentedBasis) = false +Base.in(ad::Augmented, ab::AugmentedBasis) = ad.elt in ab.basis + +function Base.getindex(ab::AugmentedBasis{T,I,A}, x::A) where {T,I,A} + @assert x.elt in object(ab) + return x +end + +mstructure(db::AugmentedBasis) = AugmentedMStructure(mstructure(db.basis)) + +struct AugmentedMStructure{M<:DiracMStructure} <: MultiplicativeStructure + op::M +end + +function (mstr::AugmentedMStructure)(aδx::Augmented, aδy::Augmented) + δxy = first(keys(mstr.op(aδx.elt, aδy.elt))) + if isone(δxy) + return SparseCoefficients((aδx, aδy), (-1, -1)) + else + aδxy = Augmented(δxy) + #(x-1)*(y-1) = 1 - x - y + xy = -1·(x-1) - 1·(y-1) + 1·(xy-1) + return SparseCoefficients((aδx, aδy, aδxy), (-1, -1, 1)) + end +end + +function coeffs!( + res::AbstractCoefficients, + cfs::AbstractCoefficients, + source::DiracBasis, + target::AugmentedBasis, +) + s = aug(cfs) + if !iszero(s) + throw( + "Conversion to $target not possible due to non-zero augmentation: $s", + ) + end + for (k, v) in nonzero_pairs(cfs) + isone(k) && continue + x = source[k] + MA.operate!( + UnsafeAddMul(*), + res, + v, + SparseCoefficients((target[Augmented(x)],), (1,)), + ) + end + MA.operate!(canonical, res) + return res +end diff --git a/src/mstructures.jl b/src/mstructures.jl index c8e4532..2cf3e73 100644 --- a/src/mstructures.jl +++ b/src/mstructures.jl @@ -1,40 +1,90 @@ -abstract type MultiplicativeStructure{I} <: AbstractMatrix{I} end - -struct ProductNotDefined <: Exception +struct ProductNotWellDefined <: Exception i::Any j::Any msg::Any end -function Base.showerror(io::IO, ex::ProductNotDefined) - print(io, "Product of elements $(ex.i) and $(ex.j) is not defined on the basis") +function Base.showerror(io::IO, ex::ProductNotWellDefined) + print( + io, + "Product of elements $(ex.i) and $(ex.j) is not defined on the basis", + ) print(io, " or the multiplicative structure could not be completed") if isdefined(ex, :msg) print(io, ": $(ex.msg)") end - print(io, ".") + return print(io, ".") end -struct TrivialMStructure{I,B<:AbstractBasis} <: MultiplicativeStructure{I} - basis::B -end +""" + MultiplicativeStructure +Structure representing multiplication w.r.t its basis. + +Implements +* `basis(ms::MultiplicativeStructure{T}) → AbstractBasis{T}` +* `Basis.getindex(ms::MultiplicativeStructure{T}, i::T, j::T) → + Union{AbstractCoefficients, AbstractVector}` + the product of `i` and `j` represented by coefficients in `basis(ms)`. -TrivialMStructure(basis::AbstractBasis{T,I}) where {T,I} = - TrivialMStructure{I,typeof(basis)}(basis) +When the product is not representable faithfully, + `ProductNotWellDefined` exception should be thrown. +""" +abstract type MultiplicativeStructure end + +""" + struct UnsafeAddMul{M<:Union{typeof(*),MultiplicativeStructure}} + structure::M + end -basis(mstr::TrivialMStructure) = mstr.basis -Base.size(mstr::TrivialMStructure) = (l = length(basis(mstr)); (l, l)) -_get(mstr::TrivialMStructure, i) = i ≥ 0 ? i : (b = basis(mstr); b[star(b[-i])]) +The value of `(op::UnsafeAddMul)(a, b, c)` is `a + structure(b, c)` +where `a` is not expected to be canonicalized before the operation `+` +and should not be expected to be canonicalized after either. +""" +struct UnsafeAddMul{M<:Union{typeof(*),MultiplicativeStructure}} + structure::M +end + +function MA.operate_to!(res, ms::MultiplicativeStructure, v, w) + if res === v || res === w + throw(ArgumentError("No alias allowed")) + end + MA.operate!(zero, res) + MA.operate!(UnsafeAddMul(ms), res, v, w) + MA.operate!(canonical, res) + return res +end -Base.@propagate_inbounds function Base.getindex( - mstr::TrivialMStructure, - i::Integer, - j::Integer, +function MA.operate!( + ::UnsafeAddMul{typeof(*)}, + mc::SparseCoefficients, + val, + c::AbstractCoefficients, ) - b = basis(mstr) - i, j = _get(mstr, i), _get(mstr, j) - g, h = b[i], b[j] - gh = g * h - gh in b || throw(ProductNotDefined(i, j, "$g · $h = $gh")) - return basis(mstr)[gh] + append!(mc.basis_elements, keys(c)) + vals = values(c) + if vals isa AbstractVector + append!(mc.values, val .* vals) + else + append!(mc.values, val * collect(values(c))) + end + return mc +end + +function MA.operate!(ms::UnsafeAddMul, res, v, w) + for (kv, a) in nonzero_pairs(v) + for (kw, b) in nonzero_pairs(w) + c = ms.structure(kv, kw) + MA.operate!(UnsafeAddMul(*), res, a * b, c) + end + end + return res +end + +struct DiracMStructure{Op} <: MultiplicativeStructure + op::Op +end + +function (mstr::DiracMStructure)(x::T, y::T) where {T} + xy = mstr.op(x, y) + return SparseCoefficients((xy,), (1,)) end diff --git a/src/mtables.jl b/src/mtables.jl index 9410dd6..6b9fac0 100644 --- a/src/mtables.jl +++ b/src/mtables.jl @@ -1,148 +1,154 @@ -abstract type AbstractMTable{I} <: MultiplicativeStructure{I} end - -function _check(product_matrix::AbstractMatrix, basis::AbstractVector) - idx = findfirst(iszero, product_matrix) - if idx !== nothing - i, j = Tuple(idx) - x, y = basis[i], basis[j] - throw( - ProductNotDefined( - i, - j, - "$x · $y = $(x * y)", - ), - ) - end - return true +""" + MTable{T, I} <: MultiplicativeStructure{T} +Multiplicative table, stored explicitly as an AbstractMatrix{I}. + +!!! note + Accessing `mt[i,j]` with negative indices returns the product of + `star`ed basis elements, i.e. + ```julia + mt[-i, j] == b[star(b[i])*b[j]] + ``` +""" +struct MTable{T,I,V<:AbstractVector,M<:AbstractMatrix,Ms} <: + MultiplicativeStructure + elts::V + relts::Dict{T,I} + starof::Vector{I} + table::M + mstr::Ms + lock::Base.Threads.SpinLock end -Base.size(mt::AbstractMTable) = size(mt.table) -function _star_of(basis::AbstractBasis, len::Integer) - return [basis[star(basis[i])] for i in 1:len] +function MTable( + elts::AbstractVector, + mstr::MultiplicativeStructure, + dims::NTuple{2,I}, +) where {I<:Integer} + Base.require_one_based_indexing(elts) + @assert length(elts) ≥ first(dims) + @assert length(elts) ≥ last(dims) + + relts = Dict(b => I(idx) for (idx, b) in pairs(elts)) + starof = [relts[star(x)] for x in elts] + T = typeof(mstr(first(elts), first(elts))) + table = Matrix{T}(undef, dims) + @assert !isbitstype(T) || dims == (0, 0) + + return MTable(elts, relts, starof, table, mstr, Base.Threads.SpinLock()) end -## MTables -struct MTable{I,M<:AbstractMatrix{I}} <: AbstractMTable{I} - table::M - star_of::Vector{I} +Base.@propagate_inbounds function __absindex(mt::MTable, i::Integer) + return ifelse(i > 0, i, oftype(i, mt.starof[abs(i)])) end -function MTable(basis::AbstractBasis; table_size) - @assert length(table_size) == 2 - - @assert 1 <= first(table_size) <= length(basis) - @assert 1 <= last(table_size) <= length(basis) - - table = zeros(SparseArrays.indtype(basis), table_size) - - complete!(table, basis) - _check(table, basis) +Base.size(mt::MTable, i::Vararg) = size(mt.table, i...) +Base.haskey(mt::MTable, x) = haskey(mt.relts, x) +Base.getindex(mt::MTable{T}, x::T) where {T} = mt.relts[x] +Base.getindex(mt::MTable, i::Integer) = mt.elts[__absindex(mt, i)] - return MTable(table, _star_of(basis, max(table_size...))) +function __iscomputed(mt::MTable, i, j) + return isassigned(mt.table, i, j) end -function complete!(table::AbstractMatrix, basis::AbstractBasis, lck=Threads.SpinLock()) - Threads.@threads for j in axes(table, 2) - y = basis[j] - for i in axes(table, 1) - xy = basis[i] * y - lock(lck) do - table[i, j] = basis[xy] - end - end - end - return table +function (mt::MTable{T})(x::T, y::T) where {T} + i = __absindex(mt, mt[x]) + j = __absindex(mt, mt[y]) + return checkbounds(Bool, mt.table, i, j) ? mt(i, j) : mt.mstr(x, y) end -function complete!(table::Matrix, basis::AbstractBasis, lck=Threads.SpinLock()) - Threads.@threads for j in axes(table, 2) - y = basis[j] - for i in axes(table, 1) - x = basis[i] - table[i, j] = basis[x*y] +function (mt::MTable)(i::Integer, j::Integer) + if !checkbounds(Bool, mt.table, i, j) + x, y = mt[i], mt[j] + return mt.mstr(x, y) + end + if !__iscomputed(mt, i, j) + lock(mt.lock) do + return thread_unsafe_complete!(mt, i, j) end end - return table -end - -basis(mt::MTable) = throw("No basis is defined for a simple $(typeof(mt))") -Base.@propagate_inbounds _iscached(mt::MTable, i, j) = !iszero(mt.table[i, j]) -Base.@propagate_inbounds _get(cmt::MTable, i::Integer) = ifelse(i ≥ 0, i, cmt.star_of[abs(i)]) -Base.@propagate_inbounds function Base.getindex(m::MTable, i::Integer, j::Integer) - @boundscheck checkbounds(Bool, m, abs(i), abs(j)) || - throw(ProductNotDefined(i, j, "out of Mtable bounds")) - @boundscheck !_iscached(m, abs(i), abs(j)) && throw(ProductNotDefined(i, j, "product not stored")) - @inbounds begin - i = _get(m, i) - j = _get(m, j) - - return m.table[i, j] + res = mt.table[i, j] # load + while res != mt.table[i, j] # compare + res = mt.table[i, j] # and load again + yield() end -end - -## CachedMTables -struct CachedMTable{I,T,B<:AbstractBasis{T,I},M<:MTable} <: AbstractMTable{I} - basis::B - table::M - lock::Threads.SpinLock + return mt.table[i, j] # and load again end -function CachedMTable(basis::AbstractBasis{T,I}; table_size) where {T,I} - return CachedMTable(basis, zeros(I, table_size)) +function thread_unsafe_complete!(mt::MTable, i::Integer, j::Integer) + g, h = mt[i], mt[j] + gh_cfs = mt.mstr(g, h) + mt.table[i, j] = gh_cfs + return mt end -function CachedMTable( - basis::AbstractBasis{T,I}, - mt::AbstractMatrix{<:Integer}, -) where {T,I} - star_of = _star_of(basis, max(size(mt)...)) - mtable = MTable(mt, star_of) - return CachedMTable(basis, mtable, Threads.SpinLock()) +function complete!(mt::MTable) + Threads.@threads for j in axes(mt.table, 2) + for i in axes(mt.table, 1) + @inbounds thread_unsafe_complete!(mt, i, j) + end + end + return mt end -basis(m::CachedMTable) = m.basis -Base.@propagate_inbounds _iscached(cmt::CachedMTable, i, j) = _iscached(cmt.table, i, j) -Base.@propagate_inbounds _get(cmt::CachedMTable, i::Integer) = _get(cmt.table, i) - -Base.@propagate_inbounds function Base.getindex(cmt::CachedMTable, i::Integer, j::Integer) - @boundscheck checkbounds(Bool, cmt, abs(i), abs(j)) || - throw(ProductNotDefined(i, j, "out of Mtable bounds")) - @inbounds begin - i = _get(cmt, i) - j = _get(cmt, j) - if !_iscached(cmt, i, j) - cache!(cmt, i, j) +function MA.operate!( + ms::UnsafeAddMul{<:MTable}, + res::AbstractCoefficients, + v::AbstractCoefficients, + w::AbstractCoefficients, +) + for (kv, a) in nonzero_pairs(v) + for (kw, b) in nonzero_pairs(w) + c = ms.structure(kv, kw) + for (k, v) in nonzero_pairs(c) + res[ms.structure[k]] += v * a * b + end end - return cmt.table[i, j] end + return res end -Base.@propagate_inbounds function cache!(cmt::CachedMTable, i::Integer, j::Integer) - b = basis(cmt) - g, h = b[i], b[j] - gh = g * h - gh in b || throw(ProductNotDefined(i, j, "$g · $h = $gh")) - lock(cmt.lock) do - cmt.table.table[i, j] = b[gh] +function MA.operate!( + ms::UnsafeAddMul{<:MTable}, + res::AbstractSparseVector, + v::AbstractVector, + w::AbstractVector, +) + l1 = issparse(v) ? nnz(v) : 2^8 + l2 = issparse(w) ? nnz(w) : 2^8 + l = l1 * l2 + idcs = Vector{key_type(res)}() + vals = Vector{eltype(res)}() + sizehint!(idcs, l) + sizehint!(vals, l) + + for (kv, a) in nonzero_pairs(v) + for (kw, b) in nonzero_pairs(w) + c = ms.structure(kv, kw) + for (k, v) in nonzero_pairs(c) + push!(idcs, ms.structure[k]) + push!(vals, v * a * b) + end + end end - return cmt + res .+= sparsevec(idcs, vals, length(res)) + return res end -function cache!( - cmt::CachedMTable, - suppX::AbstractVector{<:Integer}, - suppY::AbstractVector{<:Integer}, +function MA.operate!( + ms::UnsafeAddMul{<:MTable}, + res::AbstractVector, + v::AbstractVector, + w::AbstractVector, ) - Threads.@threads for j in suppY - for i in suppX - if !_iscached(cmt, i, j) - cache!(cmt, i, j) + for (kv, a) in nonzero_pairs(v) + for (kw, b) in nonzero_pairs(w) + c = ms.structure(kv, kw) + for (k, v) in nonzero_pairs(c) + res[ms.structure[k]] += v * a * b end end end - return cmt + return res end - -complete!(cmt::CachedMTable) = cache!(cmt, 1:size(cmt, 1), 1:size(cmt, 2)) diff --git a/src/show.jl b/src/show.jl index 82301ed..f9b0afe 100644 --- a/src/show.jl +++ b/src/show.jl @@ -9,7 +9,7 @@ __needs_parens(::Any) = false __needs_parens(a::AlgebraElement) = true function _coeff_elt_print(io, c, elt) - print(io, c, "·") + print(io, c, '·') __needs_parens(elt) && print(io, '(') print(io, elt) __needs_parens(elt) && print(io, ')') @@ -19,18 +19,15 @@ end function Base.show(io::IO, a::AlgebraElement) A = parent(a) if iszero(a) - T = eltype(a) - if hasbasis(A) - _coeff_elt_print(io, zero(T), first(basis(A))) - else - print(io, zero(T)) - end - elseif hasbasis(A) - nzeros = findall(!iszero, coeffs(a)) - for (counter, idx) in enumerate(nzeros) - c, elt = coeffs(a)[idx], basis(A)[idx] - if counter == 1 + T = valtype(coeffs(a)) + _coeff_elt_print(io, zero(T), first(basis(A))) + else + _first = true + for (idx, value) in nonzero_pairs(coeffs(a)) + c, elt = value, basis(A)[idx] + if _first _coeff_elt_print(io, c, elt) + _first = false else if __prints_with_minus(c) print(io, ' ') @@ -40,14 +37,5 @@ function Base.show(io::IO, a::AlgebraElement) _coeff_elt_print(io, c, elt) end end - else - println(io, "algebra element without defined basis") - show(io, MIME("text/plain"), a.coeffs) end end - -function Base.show(io::IO, ::MIME"text/plain", mstr::TrivialMStructure) - l = length(basis(mstr)) - print(io, "TrivialMStructure over basis with $l elements") - return -end diff --git a/src/sparse_coeffs.jl b/src/sparse_coeffs.jl new file mode 100644 index 0000000..5f07251 --- /dev/null +++ b/src/sparse_coeffs.jl @@ -0,0 +1,157 @@ +struct SparseCoefficients{K,V,Vk,Vv} <: AbstractCoefficients{K,V} + basis_elements::Vk + values::Vv +end + +function SparseCoefficients(elts::Ks, vals::Vs) where {Ks,Vs} + return SparseCoefficients{eltype(elts),eltype(vals),Ks,Vs}(elts, vals) +end + +Base.keys(sc::SparseCoefficients) = sc.basis_elements +Base.values(sc::SparseCoefficients) = sc.values +function Base.copy(sc::SparseCoefficients) + return SparseCoefficients(copy(keys(sc)), copy(values(sc))) +end + +function Base.getindex(sc::SparseCoefficients{K}, key::K) where {K} + k = searchsortedfirst(sc.basis_elements, key; lt = comparable(K)) + if k in eachindex(sc.basis_elements) + v = sc.values[k] + return ifelse(sc.basis_elements[k] == key, v, zero(v)) + else + return zero(valtype(sc)) + end +end + +function Base.setindex!(sc::SparseCoefficients{K}, val, key::K) where {K} + k = searchsortedfirst(sc.basis_elements, key; lt = comparable(K)) + if k in eachindex(sc.basis_elements) && sc.basis_elements[k] == key + sc.values[k] += val + else + insert!(sc.basis_elements, k, key) + insert!(sc.values, k, val) + end + return sc +end + +function Base.zero(sc::SparseCoefficients) + return SparseCoefficients(empty(keys(sc)), empty(values(sc))) +end + +function Base.similar(s::SparseCoefficients, ::Type{T} = valtype(s)) where {T} + return SparseCoefficients(similar(s.basis_elements), similar(s.values, T)) +end + +function MA.mutability( + ::Type{<:SparseCoefficients{K,V,Vk,Vv}}, + ::typeof(canonical), + arg::Vararg{Type}, +) where {K,V,Vk,Vv} + # return MA.IsMutable() + return MA.mutability(Vk) +end + +### temporary convenience? how to handle this? +function __prealloc(X::SparseCoefficients, a::Number, op) + T = MA.promote_operation(op, valtype(X), typeof(a)) + return similar(X, T) +end + +function __prealloc(X::SparseCoefficients, Y::SparseCoefficients, op) + # this is not even correct for op = * + T = MA.promote_operation(op, valtype(X), valtype(Y)) + return similar(X, T) +end + +comparable(::Type) = isless +function MA.operate!(::typeof(canonical), res::SparseCoefficients) + return MA.operate!(canonical, res, comparable(key_type(res))) +end + +# `::C` is needed to force Julia specialize on the function type +# Otherwise, we get one allocation when we call `issorted` +# See https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing +function MA.operate!(::typeof(canonical), res::SparseCoefficients, cmp::C) where {C} + sorted = issorted(res.basis_elements; lt = cmp) + distinct = allunique(res.basis_elements) + if sorted && distinct && !any(iszero, res.values) + return res + end + + if !sorted + p = sortperm(res.basis_elements; lt = cmp) + permute!(res.basis_elements, p) + permute!(res.values, p) + end + + todelete = BitSet() + for i in firstindex(res.basis_elements):lastindex(res.basis_elements)-1 + if iszero(res.values[i]) + push!(todelete, i) + elseif res.basis_elements[i] == res.basis_elements[i+1] + res.values[i+1] += res.values[i] + push!(todelete, i) + end + end + if iszero(last(values(res))) + push!(todelete, lastindex(res.basis_elements)) + end + deleteat!(res.basis_elements, todelete) + deleteat!(res.values, todelete) + return res +end + +# arithmetic on coefficients; performance overloads + +function MA.operate!(::typeof(zero), s::SparseCoefficients) + empty!(s.basis_elements) + empty!(s.values) + return s +end + +function MA.operate_to!( + res::SparseCoefficients, + ::typeof(-), + X::SparseCoefficients, +) + return MA.operate_to!(res, *, X, -1) +end + +function MA.operate_to!( + res::SparseCoefficients, + ::typeof(*), + X::SparseCoefficients, + a::Number, +) + if res === X + res.values .*= a + else + resize!(res.basis_elements, length(X.basis_elements)) + resize!(res.values, length(res.basis_elements)) + res.basis_elements .= X.basis_elements + res.values .= a .* X.values + end + + return res +end + +function MA.operate_to!( + res::SparseCoefficients, + ::typeof(+), + X::SparseCoefficients, + Y::SparseCoefficients, +) + if res === X + append!(res.basis_elements, Y.basis_elements) + append!(res.values, Y.values) + elseif res === Y + append!(res.basis_elements, X.basis_elements) + append!(res.values, X.values) + else + MA.operate!(zero, res) + append!(res.basis_elements, X.basis_elements, Y.basis_elements) + append!(res.values, X.values, Y.values) + end + MA.operate!(canonical, res) + return res +end diff --git a/src/star.jl b/src/star.jl new file mode 100644 index 0000000..b0e971d --- /dev/null +++ b/src/star.jl @@ -0,0 +1,23 @@ +Base.adjoint(a::AlgebraElement) = star(a) +star(x::Any) = x' + +function star(X::AlgebraElement) + res = star(basis(parent(X)), coeffs(X)) + return AlgebraElement(res, parent(X)) +end + +star(::AbstractBasis, x) = star(x) + +function star(basis::AbstractBasis, d::SparseCoefficients) + k = star.(Ref(basis), keys(d)) + v = star.(values(d)) + return SparseCoefficients(k, v) +end + +function star(basis::FixedBasis, coeffs::SparseVector) + nzidcs = mstructure(basis).starof[SparseArrays.nonzeroinds(coeffs)] + nzvals = star.(SparseArrays.nonzeros(coeffs)) + + v = SparseVector(length(coeffs), nzidcs, nzvals) + return v +end diff --git a/src/types.jl b/src/types.jl index 1fc0039..11a37d5 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,132 +1,117 @@ abstract type AbstractStarAlgebra{O,T} end -struct StarAlgebra{O,T,M<:MultiplicativeStructure,B<:AbstractBasis{T}} <: - AbstractStarAlgebra{O,T} +mstructure(A::AbstractStarAlgebra) = mstructure(basis(A)) + +function _sanity_checks(coeffs, A::AbstractStarAlgebra) + @assert key_type(coeffs) == key_type(basis(A)) +end +function _sanity_checks(coeffs::AbstractVector, A::AbstractStarAlgebra) + @assert Base.haslength(basis(A)) + @assert length(coeffs) == length(basis(A)) +end + +# concrete implementation +struct StarAlgebra{O,T,B<:AbstractBasis{T}} <: AbstractStarAlgebra{O,T} object::O - mstructure::M basis::B - function StarAlgebra(obj, basis::AbstractBasis, mstr::MultiplicativeStructure) + function StarAlgebra(obj, basis::AbstractBasis) O = typeof(obj) T = eltype(basis) - M = typeof(mstr) B = typeof(basis) - return new{O,T,M,B}(obj, mstr, basis) - end - - function StarAlgebra(obj, mstr::MultiplicativeStructure) - O = typeof(obj) - T = eltype(obj) - M = typeof(mstr) - B = Basis{T,eltype(mstr)} - - return new{O,T,M,B}(obj, mstr) + return new{O,T,B}(obj, basis) end end -# TrivialMStructure: -function StarAlgebra(obj, basis::AbstractBasis) - mstr = TrivialMStructure(basis) - return StarAlgebra(obj, basis, mstr) -end - -# CachedMStructure: -function StarAlgebra( - obj, - basis::AbstractBasis, - cache_size::Tuple{<:Integer,Integer}; - precompute=false -) - mstr = CachedMTable(basis, table_size=cache_size) - precompute && complete!(mstr) - return StarAlgebra(obj, basis, mstr) -end - -hasbasis(A::StarAlgebra) = isdefined(A, :basis) - basis(A::StarAlgebra) = A.basis object(A::StarAlgebra) = A.object -# Base.eltype(A::StarAlgebra{O,B}) where {O,B} = eltype(B) -struct AlgebraElement{A,T,V<:AbstractVector{T}} +struct AlgebraElement{A,T,V} <: MA.AbstractMutable coeffs::V parent::A - - function AlgebraElement(coeffs::AbstractVector, A::AbstractStarAlgebra) - if hasbasis(A) - @assert length(coeffs) == length(basis(A)) - end - return new{typeof(A),eltype(coeffs),typeof(coeffs)}(coeffs, A) - end end -coeffs(a::AlgebraElement) = a.coeffs Base.parent(a::AlgebraElement) = a.parent -Base.eltype(a::AlgebraElement) = eltype(coeffs(a)) +Base.eltype(a::AlgebraElement) = valtype(coeffs(a)) +coeffs(a::AlgebraElement) = a.coeffs +function coeffs(x::AlgebraElement, b::AbstractBasis) + return coeffs(coeffs(x), basis(parent(x)), b) +end + +function AlgebraElement(coeffs, A::AbstractStarAlgebra) + _sanity_checks(coeffs, A) + return AlgebraElement{typeof(A),valtype(coeffs),typeof(coeffs)}(coeffs, A) +end + +function AlgebraElement( + coeffs::SparseCoefficients{T}, + A::AbstractStarAlgebra{O,T}, +) where {O,T} + return AlgebraElement{typeof(A),valtype(coeffs),typeof(coeffs)}(coeffs, A) +end ### constructing elements Base.zero(A::AbstractStarAlgebra) = zero(Int, A) function Base.zero(T::Type, A::AbstractStarAlgebra) - if hasbasis(A) - I = SparseArrays.indtype(basis(A)) - return AlgebraElement(sparsevec(I[], T[], length(basis(A))), A) - end - throw( - "Algebra without basis; use the `AlgebraElement` constructor directly.", - ) + cfs = zero_coeffs(T, basis(A)) + return AlgebraElement(cfs, A) end Base.one(A::AbstractStarAlgebra) = one(Int, A) function Base.one(T::Type, A::AbstractStarAlgebra) - if hasbasis(A) - b = basis(A) - i = b[one(object(A))] - return AlgebraElement(sparsevec([i], [one(T)], length(b)), A) + i = one(object(A)) + sc = SparseCoefficients([i], [one(T)]) + # TODO + # this is not correct, more thought is needed + if basis(A) isa DiracBasis + @assert i in basis(A) + return AlgebraElement(sc, A) + else + return AlgebraElement( + coeffs(sc, DiracBasis{UInt}(object(A)), basis(A)), + A, + ) end - throw( - "Algebra without basis; use the `AlgebraElement` constructor directly.", - ) end -Base.zero(a::AlgebraElement) = (b = similar(a); return zero!(b)) -Base.one(a::AlgebraElement) = one(parent(a)) +Base.zero(a::AlgebraElement) = (b = similar(a); return MA.operate!(zero, b)) +Base.one(a::AlgebraElement) = one(eltype(a), parent(a)) Base.iszero(a::AlgebraElement) = iszero(coeffs(a)) function Base.isone(a::AlgebraElement) - b = basis(parent(a)) - k = findfirst(!iszero, coeffs(a)) - k === nothing && return false - isone(a[k]) || return false - isone(b[k]) || return false - return isnothing(findnext(!iszero, coeffs(a), k + 1)) -end + c = coeffs(a) + A = parent(a) + cfs1 = SparseCoefficients((one(object(A)),), (1,)) -function (A::AbstractStarAlgebra{O,T})(elt::T) where {O,T} - if hasbasis(A) - b = basis(A) - i = b[elt] - return AlgebraElement(sparsevec([i], [1], length(b)), A) + if basis(A) isa DiracBasis + return c == cfs1 else - throw("Algebra without basis: cannot coerce $elt") + dc = coeffs(c, basis(parent(a)), DiracBasis{UInt}(object(parent(a)))) + return dc == cfs1 end end -function (A::AbstractStarAlgebra)(x::Number) - if hasbasis(A) - b = basis(A) - i = b[one(object(A))] - return AlgebraElement(sparsevec([i], [x], length(b)), A) - else - throw("Algebra without basis: cannot coerce $x") - end +function (A::AbstractStarAlgebra{O,T})(elt::T) where {O,T} + b = basis(A) + @assert elt in b + res = zero_coeffs(Int, b) + res[b[elt]] = 1 + return AlgebraElement(res, A) end -Base.similar(X::AlgebraElement, T=eltype(X)) = AlgebraElement(similar(coeffs(X), T), parent(X)) +(A::AbstractStarAlgebra)(x::Number) = x * one(A) + +function Base.similar(X::AlgebraElement, T = eltype(X)) + return AlgebraElement(similar(coeffs(X), T), parent(X)) +end function AlgebraElement{T}(X::AlgebraElement) where {T} v = coeffs(X) w = similar(v, T) - w .= v + MA.operate!(zero, w) + for (k, v) in nonzero_pairs(v) + w[k] = v + end return AlgebraElement(w, parent(X)) end diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 0000000..cc6317b --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,15 @@ +[deps] +Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557" +GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +Groups = "0.8" +GroupsCore = "0.5" +PermutationGroups = "0.6" diff --git a/test/abstract_coeffs.jl b/test/abstract_coeffs.jl new file mode 100644 index 0000000..efc6826 --- /dev/null +++ b/test/abstract_coeffs.jl @@ -0,0 +1,59 @@ +@testset "Abstract coefficients" begin + G = PermGroup(perm"(1,2,3)", perm"(1,2)") + RG = StarAlgebra(G, SA.DiracBasis{UInt8}(G)) + fRG = let RG = RG, n = length(basis(RG)) + fb = SA.FixedBasis(basis(RG); n = n) + StarAlgebra(SA.object(RG), fb) + end + + h = Permutation(perm"(2,3)", G) + + α, β = let a = fRG(1), b = fRG(h) + x = AlgebraElement(ACoeffs(coeffs(a, basis(fRG))), fRG) + y = AlgebraElement(ACoeffs(coeffs(b, basis(fRG))), fRG) + x, y + end + @test coeffs(2α) isa ACoeffs{Int} + @test coeffs(α - β) isa ACoeffs + @test coeffs(α - β // 3) isa ACoeffs{<:Rational} + @test 3(α // 3) == α + + l = length(basis(RG)) + P = sum(RG(g) for g in basis(RG)) // l + fP = AlgebraElement(ACoeffs(coeffs(P, basis(fRG))), fRG) + @test AlgebraElement(coeffs(fP, basis(RG)), RG) == P + @test fP * fP == fP + + fP2 = (α + β) // 2 + @test coeffs(fP2) isa ACoeffs{<:Rational} + + @test SA.key_type(coeffs(fP2)) == UInt32 + @test zero(coeffs(fP2)) == ACoeffs(coeffs(zero(P), basis(fRG))) + + P3 = 2 * sum(RG(g) for g in basis(RG) if sign(g) > 0) // l + fP3 = AlgebraElement(ACoeffs(coeffs(P3, basis(fRG))), fRG) + @test AlgebraElement(coeffs(fP3, basis(RG)), RG) == P3 + @test fP3 * fP3 == fP3 + + PAlt = sum(sign(g) * RG(g) for g in basis(RG)) // l + fPAlt = AlgebraElement(ACoeffs(coeffs(PAlt, basis(fRG))), fRG) + @test AlgebraElement(coeffs(fPAlt, basis(RG)), RG) == PAlt + @test fPAlt * fPAlt == fPAlt + + @test fP3 * fPAlt == fPAlt * fP3 + + @test fP2 * fP2 == fP2 + + @test fP2 * fP3 == fP3 * fP2 == fP + + P2m = (RG(1) - RG(h)) // 2 + fP2m = AlgebraElement(ACoeffs(coeffs(P2m, basis(fRG))), fRG) + @test fP2m * fP2m == fP2m + + @test fP2m * fP3 == fP3 * fP2m == fPAlt + @test iszero(fP2m * fP2) + + @test norm(fP2m) == norm(P2m) == norm(fP2m) + v = coeffs(P2m, basis(fRG)) # an honest vector + @test dot(fP2m, fP2m) == dot(coeffs(fP2m), v) == dot(v, coeffs(fP2m)) +end diff --git a/test/arithmetic.jl b/test/arithmetic.jl deleted file mode 100644 index 19b71e9..0000000 --- a/test/arithmetic.jl +++ /dev/null @@ -1,300 +0,0 @@ -@testset "Arithmetic" begin - G = PermGroup(perm"(1,2,3)", perm"(1,2)") - b = StarAlgebras.Basis{UInt8}(collect(G)) - l = length(b) - RG = StarAlgebra(G, b, (l, l)) - - @test contains(sprint(show, RG), "*-algebra of") - - @testset "Module structure" begin - a = AlgebraElement(ones(Int, order(G)), RG) - - @test -a isa AlgebraElement - @test coeffs(-a) == -coeffs(a) - - @test 2 * a isa AlgebraElement - @test eltype(2 * a) == typeof(2) - @test coeffs(2 * a) == 2coeffs(a) - - @test eltype(2a) == Int - y = div(2a, 2) - @test y == a - @test eltype(y) == Int - - @test 2.0 * a isa AlgebraElement - @test eltype(2.0 * a) == typeof(2.0) - @test coeffs(2.0 * a) == 2.0 * coeffs(a) - - @test coeffs(a / 2) == coeffs(a) / 2 - b = a / 2 - @test b isa AlgebraElement - @test eltype(b) == typeof(1 / 2) - @test coeffs(b / 2) == 0.25 * coeffs(a) - - c = a // 1 - - @test eltype(c) == Rational{Int} - @test c // 4 isa AlgebraElement - @test c // big(4) isa AlgebraElement - @test eltype(c // (big(4) // 1)) == Rational{BigInt} - - @test (1.0a) * 1 // 2 == (0.5a) == c // 2 - end - - @testset "Additive structure" begin - a = AlgebraElement(ones(Int, order(G)), RG) - b = sum(sign(g) * RG(g) for g in G) - - @test a == - AlgebraElement(ones(Int, order(G)), RG) == - sum(RG(g) for g in G) - - @test 1 / 2 * (coeffs(a + b)) == [1.0, 0.0, 1.0, 0.0, 1.0, 0.0] - - g, h = gens(G) - k = g * h - - a = RG(1) + RG(g) + RG(h) - b = RG(1) - RG(k) - RG(h) - - @test a - b == RG(g) + RG(k) + 2RG(h) - - @test a + b - 2a == b - a - - @test 1 // 2 * 2a == a - @test a + 2a == (3 // 1) * a - @test 2a - (1 // 1) * a == a - end - - @testset "Multiplicative structure" begin - for g in G, h in G - a = RG(g) - b = RG(h) - @test a * b == RG(g * h) - @test (a + b) * (a + b) == a * a + a * b + b * a + b * b - end - - for g in G - @test star(RG(g)) == RG(inv(g)) - @test (one(RG) - RG(g)) * star(one(RG) - RG(g)) == - 2 * one(RG) - RG(g) - RG(inv(g)) - @test aug(one(RG) - RG(g)) == 0 - end - - g, h = PermutationGroups.gens(G) - k = g * h - - a = RG(1) + RG(g) + RG(h) - b = RG(1) - RG(k) - RG(h) - - @test norm(a) == norm(b) - @test LinearAlgebra.dot(a, coeffs(a)) ≈ - norm(a)^2 ≈ - LinearAlgebra.dot(coeffs(a), a) - - @test a * b == StarAlgebras.mul!(a, a, b) - - @test aug(a) == 3 - @test aug(b) == -1 - @test aug(a) * aug(b) == aug(a * b) == aug(b * a) - - z = sum((one(RG) - RG(g)) * star(one(RG) - RG(g)) for g in G) - @test aug(z) == 0 - - @test supp(z) == StarAlgebras.basis(parent(z)) - @test supp(RG(1) + RG(g)) == [one(G), g] - @test supp(a) == [one(G), g, h] - - @testset "Projections in star algebras" begin - b = StarAlgebras.basis(RG) - l = length(b) - P = sum(RG(g) for g in b) // l - @test P * P == P - - P3 = 2 * sum(RG(g) for g in b if sign(g) > 0) // l - @test P3 * P3 == P3 - - PAlt = sum(sign(g) * RG(g) for g in b) // l - @test PAlt * PAlt == PAlt - - @test P3 * PAlt == PAlt * P3 - - P2 = (RG(1) + RG(b[2])) // 2 - @test P2 * P2 == P2 - - @test P2 * P3 == P3 * P2 == P - - P2m = (RG(1) - RG(b[2])) // 2 - @test P2m * P2m == P2m - - @test P2m * P3 == P3 * P2m == PAlt - @test iszero(P2m * P2) - end - end - - @testset "Mutable API and trivial mstructure" begin - A = [:a, :b, :c] - b = StarAlgebras.Basis{UInt16}(words(A, radius=8)) - l = findfirst(w -> length(w) > 4, b) - 1 - - RG = StarAlgebra(one(first(b)), b) - - @test basis(RG) === b - @test basis(RG.mstructure) === basis(RG) - - RGc = StarAlgebra(one(first(b)), b, (l, l)) - @test basis(RGc) === b - @test basis(RGc.mstructure) === basis(RGc) - - @test all(RG.mstructure[1:121, 1:121] .== RGc.mstructure) - - Z = zero(RG) - W = zero(RGc) - - let g = b[rand(1:121)] - X = RG(g) - Y = -RG(star(g)) - for i in 1:3 - X[b[rand(1:121)]] += rand(-3:3) - Y[b[rand(1:121)]] -= rand(3:3) - end - - Xc = AlgebraElement(coeffs(X), RGc) - Yc = AlgebraElement(coeffs(Y), RGc) - - @test coeffs(X * Y) == - coeffs(Xc * Yc) == - coeffs(StarAlgebras.mul!(Z, X, Y)) - - @test coeffs(X^2) == coeffs(Xc^2) == coeffs(X * X) - @test coeffs(Y^2) == coeffs(Yc^2) == coeffs(Y * Y) - - @test coeffs(Z) == StarAlgebras.mul!( - coeffs(W), - coeffs(X), - coeffs(Y), - RG.mstructure, - ) - @test coeffs(Z) == coeffs(W) - @test coeffs(Z) == StarAlgebras.mul!( - coeffs(W), - coeffs(X), - coeffs(Y), - RGc.mstructure, - ) - @test coeffs(Z) == coeffs(W) - - StarAlgebras.zero!(W) - StarAlgebras.fmac!(coeffs(W), coeffs(X), coeffs(Y), RG.mstructure) - - @test coeffs(2 * X * Y) == coeffs(StarAlgebras.mul!(W, W, 2)) - end - end - - @testset "mutable arithmetic" begin - A = [:a, :b, :c] - bas = StarAlgebras.Basis{UInt16}(words(A, radius=4)) - l = findfirst(w -> length(w) > 2, bas) - 1 - RG = StarAlgebra(one(first(bas)), bas, (l, l)) - - a = let c = rand(-3:3, l) - resize!(c, length(bas)) - c[l:end] .= 0 - AlgebraElement(c, RG) - end - b = let c = rand(-3:3, l) - resize!(c, length(bas)) - c[l:end] .= 0 - AlgebraElement(c, RG) - end - - let d = deepcopy(a) - StarAlgebras.zero!(d) - StarAlgebras.neg!(d, a) - - d = deepcopy(a) - @test !iszero(d) - @test @allocated(StarAlgebras.zero!(d)) == 0 - @test iszero(d) - - @test @allocated(StarAlgebras.neg!(d, a)) == 0 - @test d == -a - end - - let d = deepcopy(a) - StarAlgebras.add!(d, d, b) - StarAlgebras.add!(d, b, d) - StarAlgebras.add!(d, a, b) - - d = deepcopy(a) - @test @allocated(StarAlgebras.add!(d, d, b)) == 0 - @test d == a + b - - d = deepcopy(a) - @test @allocated(StarAlgebras.add!(d, b, d)) == 0 - @test d == a + b - - @test @allocated(StarAlgebras.add!(d, a, b)) == 0 - @test d == a + b - end - - let d = deepcopy(a) - StarAlgebras.mul!(d, d, 2) - StarAlgebras.mul!(d, a, 2) - StarAlgebras.mul!(d, a, b) - d = deepcopy(a) - StarAlgebras.mul!(d, d, b) - - d = deepcopy(a) - @test @allocated(StarAlgebras.mul!(d, d, 2)) == 0 - @test d == 2a - - @test @allocated(StarAlgebras.mul!(d, a, 2)) == 0 - @test d == 2a - - @test @allocated(StarAlgebras.mul!(d, a, b)) == 32 - @test d == a * b - - d = deepcopy(a) - @test @allocated(StarAlgebras.mul!(d, d, b)) != 0 - z = StarAlgebras.mul!(d, d, b) - @test z == a * b - @test z !== d - end - end -end - -@testset "Group Algebra caching" begin - A = [:a, :b, :c] - b = StarAlgebras.Basis{UInt8}(words(A, radius=4)) - k = findfirst(w -> length(w) == 3, b) - 1 - - RG = StarAlgebra(Word(A, Int[]), b, (k, k)) - @test RG isa StarAlgebra - - D = sum(RG(b[i]) for i in 1:k) - @test D isa AlgebraElement - g = one(RG) - @test isone(g) - - @test one(RG) == g - @test iszero(zero(RG)) - @test 0 * g == zero(RG) - @test iszero(0 * g) - - h = RG(b[3]) - - @test D * one(RG) == D - @test one(RG) * D == D - - @test supp(D) == b[1:k] - - @test_throws StarAlgebras.ProductNotDefined all(!iszero, RG.mstructure.table) - - @test D * D isa AlgebraElement - - @test all(!iszero, RG.mstructure.table) - - RG = StarAlgebra(Word(A, Int[]), b, (k, k), precompute=true) - @test all(!iszero, RG.mstructure.table) -end diff --git a/test/caching_allocations.jl b/test/caching_allocations.jl new file mode 100644 index 0000000..027dc0c --- /dev/null +++ b/test/caching_allocations.jl @@ -0,0 +1,70 @@ +# `@allocations` is not defined on Julia v1.6 +@static if v"1.10" ≤ VERSION + function _alloc_test(output, op, a, b, allocs) + MA.operate_to!(output, op, a, b) + expected = op(a, b) + @test @allocations(MA.operate_to!(output, op, a, b)) <= allocs + @test output == expected + end +end + +@testset "FixedBasis caching && allocations" begin + alph = [:a, :b, :c] + A★ = FreeWords(alph) + B = SA.DiracBasis{UInt16}(A★) + + fB = SA.FixedBasis(B; n = nwords(A★, 8), mt = UInt32(nwords(A★, 4))) + fRG = StarAlgebra(A★, fB) + + k = size(SA.mstructure(basis(fRG)), 1) + + y = spzeros(length(basis(fRG))) + y[1:k] .= 1 + Y = AlgebraElement(y, fRG) + @test Y == sum(fRG(basis(fRG)[i]) for i in 1:k) + + @test Y isa AlgebraElement + + @static if v"1.10" ≤ VERSION < v"1.11" + star(Y) + star(Y) + @test (@allocations star(Y)) ≤ 4 + end + + @test SA.supp(Y) == basis(fRG)[1:k] + + @test Y * one(fRG) == Y + @test one(fRG) * Y == Y + + @test_throws SA.UndefRefError all(!iszero, SA.mstructure(fRG).table) + mstr = deepcopy(SA.mstructure(fRG)) + SA.complete!(mstr) + @test all(!iszero(mstr.table)) + @test_throws SA.UndefRefError all(!iszero, SA.mstructure(fRG).table) + + @static if v"1.10" ≤ VERSION < v"1.11" + @test (@allocations Y * Y) > k^2 - 2 * k + @test Y * Y isa AlgebraElement + @test (@allocations Y * Y) ≤ 26 + else + k1 = @allocated Y * Y + @test Y * Y isa AlgebraElement + Y * Y + k2 = @allocated Y * Y + @test k2 / k1 < 0.7 + end + + @test all(!iszero, SA.mstructure(fRG).table) + + + @static if v"1.10" ≤ VERSION < v"1.11" + YY = deepcopy(Y) + _alloc_test(YY, *, Y, Y, 25) + _alloc_test(YY, +, Y, Y, 0) + _alloc_test(YY, -, Y, Y, 0) + + # SparseArrays calls `Base.unalias` which allocates: + _alloc_test(YY, +, YY, Y, 2) + _alloc_test(YY, +, Y, YY, 2) + end +end diff --git a/test/constructors.jl b/test/constructors.jl index ed019ee..9f5025f 100644 --- a/test/constructors.jl +++ b/test/constructors.jl @@ -1,78 +1,84 @@ @testset "Algebra and Elements" begin - A = [:a, :b, :c] - b = StarAlgebras.Basis{UInt8}(words(A, radius=2)) - l = length(b) + alph = [:a, :b, :c] + A★ = FreeWords(alph) + B = SA.DiracBasis{UInt16}(A★) + RG = StarAlgebra(A★, B) - RG = StarAlgebra(one(first(b)), b, (4, 4)) + @test typeof(zero(RG)) == typeof(RG(0)) + @test typeof(one(RG)) == typeof(RG(1)) - a = rand(l) + @test one(RG) == RG(1) + @test RG(1) == one(RG(0)) - @test AlgebraElement(a, RG) isa AlgebraElement - @test all(RG(g) isa AlgebraElement{typeof(RG)} for g in b) + @test zero(RG) == RG(0) + @test RG(0) == zero(RG(1)) - @test typeof(zero(RG)) == typeof(RG(0)) - @test typeof(one(RG)) == typeof(RG(1)) + x = let A = alph, n = 7 + elts = [Word(A, rand(1:length(A), rand(0:7))) for _ in 1:n] + vals = rand(-5:5, n) + coeffs = SA.SparseCoefficients(elts, vals) + MA.operate!(SA.canonical, coeffs) + end + + @test AlgebraElement(x, RG) isa AlgebraElement - @test_throws AssertionError AlgebraElement([1, 2, 3], RG) - @test AlgebraElement([1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], RG) isa AlgebraElement + X = AlgebraElement(x, RG) - x = AlgebraElement([1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], RG) - @test AlgebraElement{Float64}(x) isa AlgebraElement - y = AlgebraElement{Float64}(x) - @test parent(x) === parent(y) - @test x == y + @test AlgebraElement{Float64}(X) isa AlgebraElement + Y = AlgebraElement{Float64}(X) + @test parent(X) === parent(Y) + @test X == Y - p = Word(A, [2, 3]) + p = Word(alph, [2, 3]) a = RG(p) - @test coeffs(a)[b[p]] == 1 - @test coeffs(a) isa SparseVector - @test all(coeffs(a)[i] == 0 for i in 1:length(b) if i ≠ b[p]) + @test coeffs(a)[B[p]] == 1 + @test coeffs(a) isa SA.SparseCoefficients + @test length(collect(SA.nonzero_pairs(coeffs(a)))) == 1 @test a(p) == 1 - @test all(a(g) == 0 for g in b if g != p) + + let b = collect(Iterators.take(B, 121)) + @test all(coeffs(a)[x] == 0 for x in b if x ≠ p) + @test all(a(g) == 0 for g in b if g != p) + end @test sprint(show, zero(RG)) == "0·(id)" @test sprint(show, one(RG)) == "1·(id)" @test isone(one(a)) - @test !isone(x) + @test !isone(a) @test iszero(zero(a)) @test sprint(show, a) == "1·b·c" @test sprint(show, -a) == "-1·b·c" @test hash(a) == hash(one(RG) + RG(p) - one(RG)) - z = zeros(l) - z[b[p]] = 1 - @test AlgebraElement(z, RG) == a - - @test StarAlgebras.supp(a) == [p] - @test StarAlgebras.supp_ind(a) == [b[p]] + @test SA.supp(a) == [p] + @test collect(SA.nonzero_pairs(coeffs(a))) == [(p => 1)] - s = one(first(b)) + s = one(first(B)) @test a(s) == 0 a[s] = 2 - @test coeffs(a)[b[s]] == 2 - @test a[b[s]] == 2 + @test coeffs(a)[s] == 2 @test a(s) == 2 - dense_a = AlgebraElement(Vector(coeffs(a)), RG) - @test a == dense_a - @test hash(a) == hash(dense_a) + # dense_a = AlgebraElement(Vector(coeffs(a)), RG) + # @test a == dense_a + # @test hash(a) == hash(dense_a) - @test StarAlgebras.supp_ind(a) == [b[s], b[p]] == StarAlgebras.supp_ind(dense_a) - @test supp(a) == [s, p] == StarAlgebras.supp(dense_a) + # @test SA.supp_ind(a) == [b[s], b[p]] == SA.supp_ind(dense_a) + # @test SA.supp(a) == [s, p] == SA.supp(dense_a) aa = a - RG(p) - dense_aa = dense_a - RG(p) - @test StarAlgebras.supp_ind(aa) == [b[s]] == StarAlgebras.supp_ind(dense_aa) - @test supp(aa) == [s] == StarAlgebras.supp(dense_aa) + # dense_aa = dense_a - RG(p) + # @test SA.supp_ind(aa) == [b[s]] == SA.supp_ind(dense_aa) + @test SA.supp(aa) == [s] @test sprint(show, a) == "2·(id) +1·b·c" @test sprint(show, -a) == "-2·(id) -1·b·c" - z[b[s]] = 2 - @test AlgebraElement(z, RG) == a - @test sprint(show, AlgebraElement(z, RG)) == "2.0·(id) +1.0·b·c" + Z = AlgebraElement{Float64}(a) + @test Z == a + @test sprint(show, Z) == "2.0·(id) +1.0·b·c" @test sprint(show, 2one(RG) - RG(p)) == "2·(id) -1·b·c" @test LinearAlgebra.norm(a, 1) == 3 @@ -86,16 +92,4 @@ @test deepcopy(a) !== a @test coeffs(deepcopy(a)) !== coeffs(a) @test parent(deepcopy(a)) === parent(a) - - @testset "without basis" begin - O = one(first(b)) - RG2 = StarAlgebra(O, RG.mstructure) - @test_throws String zero(RG2) - @test_throws String one(RG2) - @test_throws String RG2(one(O)) - @test_throws String RG2(-5) - - @test sprint(show, AlgebraElement(rand(-2:2, 6), RG2)) isa String - @test sprint(show, AlgebraElement(zeros(Int, 6), RG2)) isa String - end end diff --git a/test/group_algebra.jl b/test/group_algebra.jl new file mode 100644 index 0000000..0dc2c25 --- /dev/null +++ b/test/group_algebra.jl @@ -0,0 +1,199 @@ +@testset "Permutation group algebra: arithmetic " begin + G = PermGroup(perm"(1,2,3)", perm"(1,2)") + b = SA.DiracBasis{UInt8}(G) + RG = StarAlgebra(G, b) + + @test contains(sprint(show, RG), "*-algebra of") + + @testset "Module structure" begin + sc = SA.SparseCoefficients( + collect(SA.object(RG)), + ones(Int, length(basis(RG))), + ) + a = AlgebraElement(sc, RG) + + @test -a isa AlgebraElement + @test coeffs(-a) == -coeffs(a) + @test zero(coeffs(a)) == coeffs(zero(a)) + + @test 2 * a isa AlgebraElement + @test eltype(2 * a) == typeof(2) + @test coeffs(2 * a) == 2coeffs(a) + + @test eltype(2a) == Int + y = div(2a, 2) + @test y == a + @test eltype(y) == Int + + @test 2.0 * a isa AlgebraElement + @test eltype(2.0 * a) == typeof(2.0) + @test coeffs(2.0 * a) == 2.0 * coeffs(a) + + @test coeffs(a / 2) == coeffs(a) / 2 + b = a / 2 + @test b isa AlgebraElement + @test eltype(b) == typeof(1 / 2) + @test coeffs(b / 2) == 0.25 * coeffs(a) + + c = a // 1 + + @test eltype(c) == Rational{Int} + @test c // 4 isa AlgebraElement + @test c // big(4) isa AlgebraElement + @test eltype(c // (big(4) // 1)) == Rational{BigInt} + + @test (1.0a) * 1 // 2 == (0.5a) == c // 2 + end + + @testset "Additive structure" begin + sc = SA.SparseCoefficients( + collect(SA.object(RG)), + ones(Int, length(basis(RG))), + ) + a = AlgebraElement(sc, RG) + b = sum(sign(g) * RG(g) for g in G) + + @test a == sum(RG(g) for g in G) + + @test 1 / 2 * (coeffs(a + b)) == SA.SparseCoefficients( + collect(SA.object(RG)), + [1.0, 0.0, 1.0, 0.0, 1.0, 0.0], + ) + + g, h = PermutationGroups.gens(G) + k = g * h + + a = RG(1) + RG(g) + RG(h) + b = RG(1) - RG(k) - RG(h) + + @test a - b == RG(g) + RG(k) + 2RG(h) + + @test a + b - 2a == b - a + @test (2a) // 2 == a + @test a + 2a == 3.0a + @test 2a - a // 1 == a + @test div(3a, 2) == a + + @test coeffs(a + b - 2a) == coeffs(a) + coeffs(b) - 2coeffs(a) + @test coeffs(2a // 2) == 2coeffs(a) // 2 == coeffs(a) + @test coeffs(3a) == 3.0coeffs(a) + @test coeffs(2a) - coeffs(a // 1) == + 2coeffs(a) - coeffs(a) // 1 == + coeffs(a) + @test div(3coeffs(a), 2) == coeffs(a) + end + + @testset "Multiplicative structure" begin + for g in G, h in G + a = RG(g) + b = RG(h) + @test a * b == RG(g * h) + @test (a + b) * (a + b) == a * a + a * b + b * a + b * b + end + + for g in G + @test star(RG(g)) == RG(inv(g)) + @test (one(RG) - RG(g)) * star(one(RG) - RG(g)) == + 2 * one(RG) - RG(g) - RG(inv(g)) + @test SA.aug(one(RG) - RG(g)) == 0 + end + + g, h = PermutationGroups.gens(G) + k = g * h + + a = RG(1) + RG(g) + RG(h) + b = RG(1) - RG(k) - RG(h) + + @test norm(a) == norm(b) + + @test a * b == MA.operate_to!(similar(a), *, a, b) + + @test SA.aug(a) == 3 + @test SA.aug(b) == -1 + @test SA.aug(a) * SA.aug(b) == SA.aug(a * b) == SA.aug(b * a) + + z = sum((one(RG) - RG(g)) * star(one(RG) - RG(g)) for g in G) + @test SA.aug(z) == 0 + + @test SA.supp(z) == sort(collect(basis(parent(z)))) + @test SA.supp(RG(1) + RG(g)) == [one(G), g] + @test SA.supp(a) == [one(G), h, g] + + @testset "Projections in star algebras" begin + b = basis(RG) + l = length(b) + P = sum(RG(g) for g in b) // l + @test P * P == P + + P3 = 2 * sum(RG(g) for g in b if sign(g) > 0) // l + @test P3 * P3 == P3 + + PAlt = sum(sign(g) * RG(g) for g in b) // l + @test PAlt * PAlt == PAlt + + @test P3 * PAlt == PAlt * P3 + + h = Permutation(perm"(2,3)", G) + P2 = (RG(1) + RG(h)) // 2 + @test P2 * P2 == P2 + + @test P2 * P3 == P3 * P2 == P + + P2m = (RG(1) - RG(h)) // 2 + @test P2m * P2m == P2m + + @test P2m * P3 == P3 * P2m == PAlt + @test iszero(P2m * P2) + + @testset "Dense projections and Fixed basis" begin + n = length(basis(RG)) + fb = SA.FixedBasis(basis(RG); n = n) + fRG = StarAlgebra(SA.object(RG), fb) + @testset "compatibility: sparse/dense" begin + fP = AlgebraElement(coeffs(P, basis(fRG)), fRG) + @test fP * fP == fP + dfP = AlgebraElement(fill(1 // n, n), fRG) + @test dfP == fP + @test dfP * dfP == dfP + @test fP * dfP == dfP + @test coeffs(fP * dfP) isa DenseArray + + @test 2dfP // 2 == dfP + + @test (dfP + fP) - dfP == dfP + end + + @testset "Projections" begin + RG = fRG + b = basis(RG) + l = length(b) + P = sum(RG(g) for g in b) // l + @test P * P == P + + P3 = 2 * sum(RG(g) for g in b if sign(g) > 0) // l + @test P3 * P3 == P3 + + PAlt = sum(sign(g) * RG(g) for g in b) // l + @test PAlt * PAlt == PAlt + + @test P3 * PAlt == PAlt * P3 + + h = Permutation(perm"(2,3)", G) + P2 = (RG(1) + RG(h)) // 2 + @test P2 * P2 == P2 + + @test P2 * P3 == P3 * P2 == P + + @test_broken -RG(h) == (-1) * RG(h) + @test !iszero(RG(1) - RG(h)) + + P2m = (RG(1) - RG(h)) // 2 + @test P2m * P2m == P2m + + @test P2m * P3 == P3 * P2m == PAlt + @test iszero(P2m * P2) + end + end + end + end +end diff --git a/test/monoid_algebra.jl b/test/monoid_algebra.jl new file mode 100644 index 0000000..92b22cc --- /dev/null +++ b/test/monoid_algebra.jl @@ -0,0 +1,136 @@ +@testset "Free monoid algebra" begin + alph = [:a, :b, :c] + A★ = FreeWords(alph) + B = SA.DiracBasis{UInt16}(A★) + RG = StarAlgebra(A★, B) + @test basis(RG) === B + + # no caching + fB = SA.FixedBasis(basis(RG); n = nwords(A★, 0, 8), mt = 0) + @test fB.table.elts === fB.elts + + fRG = StarAlgebra(A★, fB) + + g = one(fRG) + @test isone(g) + + @test one(fRG) == g + @test iszero(zero(fRG)) + @test zero(g) == zero(fRG) + @test iszero(0 * g) + + @testset "Translations between bases" begin + Z = zero(RG) + fZ = zero(fRG) + + @test coeffs(fZ) isa SparseVector + + @test coeffs(Z, fB) == coeffs(fZ) + @test coeffs(fZ, B) == coeffs(Z) + @test coeffs(coeffs(fZ, B), B, fB) == coeffs(fZ) + @test coeffs(coeffs(Z, fB), fB, B) == coeffs(Z) + + let g = basis(fRG)[rand(1:121)] + X = RG(g) + Y = -RG(star(g)) + for i in 1:3 + coeffs(X)[basis(fRG)[rand(1:121)]] += rand(-3:3) + coeffs(Y)[basis(fRG)[rand(1:121)]] -= rand(3:3) + end + + @test coeffs(X) == + coeffs(coeffs(X, basis(fRG)), basis(fRG), basis(RG)) + @test coeffs(Y) == + coeffs(coeffs(Y, basis(fRG)), basis(fRG), basis(RG)) + + fX = AlgebraElement(coeffs(X, basis(fRG)), fRG) + fY = AlgebraElement(coeffs(Y, basis(fRG)), fRG) + + @test dot(X, Y) == dot(fX, fY) == dot(coeffs(X), coeffs(Y)) + @test dot(fX, coeffs(fY)) == dot(coeffs(fX), fY) + + @test dot(X, X) ≈ norm(X)^2 ≈ dot(coeffs(X), coeffs(X)) + @test dot(fX, fX) ≈ norm(fX)^2 ≈ dot(coeffs(fX), coeffs(fX)) + + @test coeffs(fX) == + coeffs(coeffs(fX, basis(RG)), basis(RG), basis(fRG)) + @test coeffs(fY) == + coeffs(coeffs(fY, basis(RG)), basis(RG), basis(fRG)) + + @test coeffs(X * Y, basis(fRG)) == coeffs(fX * fY) + @test coeffs(X * X, basis(fRG)) == + coeffs(fX^2) == + coeffs(X^2, basis(fRG)) + + @test coeffs(Y * Y, basis(fRG)) == + coeffs(fY^2) == + coeffs(Y^2, basis(fRG)) + + MA.operate_to!(Z, *, X, Y) + MA.operate_to!(fZ, *, fX, fY) + + @test coeffs(Z) == coeffs(fZ, basis(RG)) + @test coeffs(fZ) == coeffs(Z, basis(fRG)) + + @test coeffs(2 * X * Y) == coeffs(MA.operate_to!(Z, *, Z, 2)) + @test coeffs(2 * fX * fY) == coeffs(MA.operate_to!(fZ, *, fZ, 2)) + end + end + @testset "mutable arithmetic" begin + a = let l = 12, R = 7 + support = [Word(alph, rand(1:length(alph), rand(0:R))) for _ in 1:l] + vals = rand(-3:3, l) + c = SA.SparseCoefficients(support, vals) + MA.operate!(SA.canonical, c) + AlgebraElement(c, RG) + end + b = let l = 7, R = 3 + support = [Word(alph, rand(1:length(alph), rand(0:R))) for _ in 1:l] + vals = rand(-3:3, l) + c = SA.SparseCoefficients(support, vals) + MA.operate!(SA.canonical, c) + AlgebraElement(c, RG) + end + + let d = deepcopy(a) + MA.operate!(zero, d) + MA.operate_to!(d, -, a) + + d = deepcopy(a) + @test !iszero(d) + @test @allocated(MA.operate!(zero, d)) == 0 + @test iszero(d) + + @test @allocated(MA.operate_to!(d, -, a)) == 0 + @test d == -a + end + + let d = zero(a) + MA.operate_to!(d, +, a, b) + @test d == a + b + MA.operate_to!(d, +, d, b) + @test d == a + 2b + MA.operate_to!(d, +, b, d) + @test d == a + 3b + MA.operate_to!(d, +, a, b) + @test d == a + b + end + + let d = deepcopy(a) + MA.operate_to!(d, *, a, 2) + @test d == 2a + MA.operate_to!(d, *, d, 2) + @test d == 4a + MA.operate_to!(d, *, a, b) + @test d == a * b + @test_throws ArgumentError MA.operate_to!(d, *, d, b) + + d = deepcopy(a) + @test @allocated(MA.operate_to!(d, *, d, 2)) == 0 + @test d == 2a + + @test @allocated(MA.operate_to!(d, *, a, 2)) == 0 + @test d == 2a + end + end +end diff --git a/test/mtables.jl b/test/mtables.jl index 97cc111..7c38f6a 100644 --- a/test/mtables.jl +++ b/test/mtables.jl @@ -1,5 +1,5 @@ @testset "TrivialMStructure" begin - b = StarAlgebras.Basis{UInt8}(words([:a, :b, :c], radius=4)) + b = StarAlgebras.Basis{UInt8}(words([:a, :b, :c]; radius = 4)) mstr = StarAlgebras.TrivialMStructure(b) @@ -36,9 +36,9 @@ end @testset "MTable" begin - b = StarAlgebras.Basis{UInt16}(words([:a, :b, :c, :d], radius=4)) + b = StarAlgebras.Basis{UInt16}(words([:a, :b, :c, :d]; radius = 4)) k = findfirst(w -> length(w) == 3, b) - 1 - mstr = StarAlgebras.MTable(b, table_size=(k, k)) + mstr = StarAlgebras.MTable(b; size = (k, k)) @test_throws String StarAlgebras.basis(mstr) @@ -47,25 +47,31 @@ end mstr[i, j] == b[b[i]*b[j]] for i in axes(mstr, 1) for j in axes(mstr, 2) ) @test all( - mstr[-i, j] == b[star(b[i])*b[j]] for i in axes(mstr, 1) for j in axes(mstr, 2) + mstr[-i, j] == b[star(b[i])*b[j]] for i in axes(mstr, 1) for + j in axes(mstr, 2) ) @test all( - mstr[i, -j] == b[b[i]*star(b[j])] for i in axes(mstr, 1) for j in axes(mstr, 2) + mstr[i, -j] == b[b[i]*star(b[j])] for i in axes(mstr, 1) for + j in axes(mstr, 2) ) @test all( - mstr[-i, -j] == b[star(b[i])*star(b[j])] for i in axes(mstr, 1) for j in axes(mstr, 2) + mstr[-i, -j] == b[star(b[i])*star(b[j])] for i in axes(mstr, 1) for + j in axes(mstr, 2) ) end @testset "CachedMTable" begin - b = StarAlgebras.Basis{UInt8}(words([:a, :b, :c], radius=4)) + b = StarAlgebras.Basis{UInt8}(words([:a, :b, :c]; radius = 4)) k = findfirst(w -> length(w) == 3, b) - 1 - mstr = StarAlgebras.CachedMTable(b, table_size=(k, k)) + mstr = StarAlgebras.CachedMTable(b; table_size = (k, k)) @test mstr isa StarAlgebras.CachedMTable{UInt8,Word{Symbol}} @test mstr.table.table isa Matrix{UInt8} - @test_throws StarAlgebras.ProductNotDefined StarAlgebras._check(mstr.table.table, StarAlgebras.basis(mstr)) + @test_throws StarAlgebras.ProductNotDefined StarAlgebras._check( + mstr.table.table, + StarAlgebras.basis(mstr), + ) StarAlgebras.complete!(mstr) @test all(!iszero, mstr.table) @@ -78,7 +84,7 @@ end @test mstr == mstr_sparse - mstr = StarAlgebras.CachedMTable(b, table_size=(k, k)) + mstr = StarAlgebras.CachedMTable(b; table_size = (k, k)) @test all(iszero, mstr.table.table) StarAlgebras.cache!(mstr, 1, 2) @@ -100,7 +106,7 @@ end @test_throws StarAlgebras.ProductNotDefined mstr[k+1, k] - mstr = StarAlgebras.CachedMTable(b, table_size=(k, k)) + mstr = StarAlgebras.CachedMTable(b; table_size = (k, k)) @test all(iszero, mstr.table.table) @test mstr[-1, 2] == 2 @test mstr[-2, 3] == b[star(b[2])*b[3]] diff --git a/test/perm_grp_algebra.jl b/test/perm_grp_algebra.jl new file mode 100644 index 0000000..12013a0 --- /dev/null +++ b/test/perm_grp_algebra.jl @@ -0,0 +1,98 @@ +@testset "POC: group algebra" begin + G = PermGroup(perm"(1,2,3,4,5,6)", perm"(1,2)") + g = Permutation(perm"(1,4,3,6)(2,5)", G) + h = Permutation(perm"(2,4,5,1)", G) + + db = SA.DiracBasis{UInt32}(G) + @test SA.mstructure(db) == SA.DiracMStructure(*) + @test SA.mstructure(db)(g, h) == SA.SparseCoefficients((g * h,), (1,)) + + @test db[g] == g + @test db[db[g]] == g + + xcfs = SA.SparseCoefficients([one(G), g], [1, -1]) + ycfs = SA.SparseCoefficients([one(G), inv(g)], [1, -1]) + xycfs = SA.SparseCoefficients([one(G), g, inv(g)], [2, -1, -1]) + + zcfs = SA.SparseCoefficients([one(G), h], [1, -1]) + xzcfs = SA.SparseCoefficients([one(G), g, h, g * h], [1, -1, -1, 1]) + + RG = SA.StarAlgebra(G, db) + + x = SA.AlgebraElement(xcfs, RG) + y = SA.AlgebraElement(ycfs, RG) + xy = SA.AlgebraElement(xycfs, RG) + @test x != y + @test x' == y + @test x * y == xy + + z = SA.AlgebraElement(zcfs, RG) + xz = SA.AlgebraElement(xzcfs, RG) + @test x * z == xz + + @testset "Augmented basis" begin + ad = SA.AugmentedBasis(db) + @test SA.mstructure(ad) == SA.AugmentedMStructure(SA.mstructure(db)) + @test ad[SA.Augmented(h)] isa SA.Augmented + @test sprint(show, ad[SA.Augmented(h)]) == "(-1·()+1·(1,2,4,5))" + + @test !(h in ad) + @test SA.Augmented(h) in ad + + IG = SA.StarAlgebra(G, ad) + + axcfs = SA.coeffs(x, basis(IG)) + aycfs = SA.coeffs(y, basis(IG)) + azcfs = SA.coeffs(z, basis(IG)) + ax = SA.AlgebraElement(axcfs, IG) + ay = SA.AlgebraElement(aycfs, IG) + az = SA.AlgebraElement(azcfs, IG) + + @test coeffs(ax * ay) == SA.coeffs(x * y, basis(IG)) + @test coeffs(ax * az) == SA.coeffs(x * z, basis(IG)) + @test SA.aug(ax) == 0 + @test star(ax) * star(ay) == star(ay) * star(ax) + + @test length(ad) == length(db) - 1 + @test Set(ad) == Set(SA.Augmented(g) for g in db if !isone(g)) + end + + @testset "Random elements" begin + rcfs = SA.SparseCoefficients(rand(G, 10), rand(-2:2, 10)) + r = SA.AlgebraElement(rcfs, RG) + scfs = SA.SparseCoefficients(rand(G, 10), rand(-2:2, 10)) + s = SA.AlgebraElement(scfs, RG) + + @test SA.aug(r * s) == SA.aug(r) * SA.aug(s) + end + @testset "Fixed Basis" begin + m = PermutationGroups.order(UInt16, G) + fb = SA.FixedBasis(collect(G), SA.DiracMStructure(*), (m, m)) + + @test fb[fb[g]] == g + + fRG = SA.StarAlgebra(G, fb) + + rcfs = SA.SparseCoefficients(rand(G, 10), rand(-2:2, 10)) + r = SA.AlgebraElement(rcfs, RG) + scfs = SA.SparseCoefficients(rand(G, 10), rand(-2:2, 10)) + s = SA.AlgebraElement(scfs, RG) + + @test coeffs(r, basis(fRG)) isa SparseVector + + fr = SA.AlgebraElement(coeffs(r, basis(fRG)), fRG) + fs = SA.AlgebraElement(coeffs(s, basis(fRG)), fRG) + + @test SA.aug(fr) == SA.aug(r) + @test SA.aug(fs) == SA.aug(s) + @test SA.aug(fr * fs) == SA.aug(fr) * SA.aug(fs) + + @test coeffs(r * s, basis(fRG)) isa AbstractVector + @test fr * fs == SA.AlgebraElement(coeffs(r * s, basis(fRG)), fRG) + + a, b = let mt = SA.mstructure(basis(fRG)).table + count(i -> isassigned(mt, i), eachindex(mt)), length(mt) + end + @test a ≤ b + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 4781da1..3d87d80 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,25 +1,48 @@ -using StarAlgebras using Test using LinearAlgebra using SparseArrays +using StarAlgebras +import StarAlgebras as SA +import MutableArithmetics as MA + +begin + # turning GroupsCore.GroupElement into "read-only" SparseCoefficients + import GroupsCore: GroupElement + Base.keys(g::GroupElement) = (g,) + Base.values(::GroupElement) = (1,) + SA.canonical(g::GroupElement) = g + # missing above: Base.isless for the SA.canonical + # TODO: implement sort by hashing for non-comparable elements? + + # generalities for meaningful *-algebras + SA.star(g::GroupElement) = inv(g) +end + +# implementations of the free monoid over an alphabet include("test_example_words.jl") +# example implementation of abstract coefficients +include("test_example_acoeffs.jl") + @testset "StarAlgebras" begin - include("mtables.jl") - include("constructors.jl") + # proof of concept + using PermutationGroups + include("perm_grp_algebra.jl") - using GroupsCore - StarAlgebras.star(g::GroupsCore.GroupElement) = inv(g) + include("constructors.jl") + include("group_algebra.jl") + include("abstract_coeffs.jl") - using SparseArrays - if VERSION < v"1.9" - Base.sum(v::SparseVector) = sum(nonzeros(v)) - end + # free monoid algebra + include("monoid_algebra.jl") - using PermutationGroups - include("arithmetic.jl") + include("caching_allocations.jl") + # some applications: using Groups + function Base.isless(g::Groups.FPGroupElement, h::Groups.FPGroupElement) + return isless(Groups.word(g), Groups.word(h)) + end include("sum_of_squares.jl") end diff --git a/test/sum_of_squares.jl b/test/sum_of_squares.jl index bfba90e..ed4df7b 100644 --- a/test/sum_of_squares.jl +++ b/test/sum_of_squares.jl @@ -1,99 +1,124 @@ @testset "sum of squares in FreeGroup *-algebra" begin - F = Groups.FreeGroup(4) - S = [Groups.gens(F); inv.(Groups.gens(F))] - - ID = one(F) - RADIUS = 3 - @time E_R, sizes = Groups.wlmetric_ball(S, ID, radius=2 * RADIUS) - @test sizes == [9, 65, 457, 3201, 22409, 156865] - - b = StarAlgebras.Basis{UInt32}(E_R) - - mstr = StarAlgebras.MTable(b, table_size=(sizes[RADIUS], sizes[RADIUS])) - - RG = StarAlgebra(F, b, mstr) - - g, h, k, l = S[1:4] - - G = (one(RG) - RG(g)) - @test G' == one(RG) - RG(inv(g)) - @test G' * G == StarAlgebras.mul!(zero(G), G', G) == 2one(RG) - RG(g) - RG(g)' - @test star(G * G) == G' * G' - - @testset "Sums of hermitian squares" begin - 𝕀 = one(RG) - - G = (𝕀 - RG(g)) - H = (𝕀 - RG(h)) - K = (𝕀 - RG(k)) - L = (𝕀 - RG(l)) - GH = (𝕀 - RG(g * h)) - KL = (𝕀 - RG(k * l)) - - X = (2𝕀 - star(RG(g)) - RG(h)) - Y = (2𝕀 - star(RG(g * h)) - RG(k)) - - @test -(2𝕀 - RG(g * h) - RG(g * h)') + 2G' * G + 2H' * H == X' * X - @test (2𝕀 - RG(g * h) - RG(g * h)') == GH' * GH - @test -(2𝕀 - RG(g * h * k) - RG(g * h * k)') + 2GH' * GH + 2K' * K == Y' * Y - @test -(2𝕀 - RG(g * h * k) - RG(g * h * k)') + - 2(GH' * GH - 2G' * G - 2H' * H) + - 4G' * G + - 4H' * H + - 2K' * K == Y' * Y - - @test GH' * GH - 2G' * G - 2H' * H == -X' * X - @test -(2𝕀 - RG(g * h * k) - RG(g * h * k)') + 4G' * G + 4H' * H + 2K' * K == 2X' * X + Y' * Y - - @test GH' * GH == 2G' * G + 2H' * H - (2𝕀 - RG(g)' - RG(h))' * (2𝕀 - RG(g)' - RG(h)) - @test KL' * KL == 2K' * K + 2L' * L - (2𝕀 - RG(k)' - RG(l))' * (2𝕀 - RG(k)' - RG(l)) - - @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + 2 * GH' * GH + 2 * KL' * KL == - (2𝕀 - RG(g * h)' - RG(k * l))' * (2𝕀 - RG(g * h)' - RG(k * l)) - - @test -(2𝕀 - star(RG(g * h * k * l)) - RG(g * h * k * l)) + - 2(2G' * G + 2H' * H - (2𝕀 - RG(g)' - RG(h))' * (2𝕀 - RG(g)' - RG(h))) + - 2(2K' * K + 2L' * L - (2𝕀 - RG(k)' - RG(l))' * (2𝕀 - RG(k)' - RG(l))) == - (2𝕀 - RG(g * h)' - RG(k * l))' * (2𝕀 - RG(g * h)' - RG(k * l)) - - @test -(2𝕀 - star(RG(g * h * k * l)) - RG(g * h * k * l)) + - 2(2G' * G + 2H' * H) + - 2(2K' * K + 2L' * L) == - (2𝕀 - RG(g * h)' - RG(k * l))' * (2𝕀 - RG(g * h)' - RG(k * l)) + - 2(2𝕀 - RG(g)' - RG(h))' * (2𝕀 - RG(g)' - RG(h)) + - 2(2𝕀 - RG(k)' - RG(l))' * (2𝕀 - RG(k)' - RG(l)) - - @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + - 2(2𝕀 - RG(g * h * k)' - RG(g * h * k)) + 2L' * L == - (2𝕀 - RG(g * h * k)' - RG(l))' * (2𝕀 - RG(g * h * k)' - RG(l)) - - @test 2𝕀 - RG(g * h * k)' - RG(g * h * k) == - 2GH' * GH + 2K' * K - (2𝕀 - star(RG(g * h)) - RG(k))' * (2𝕀 - star(RG(g * h)) - RG(k)) - - @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + - 2( - 2GH' * GH + 2K' * K - (2𝕀 - RG(g * h)' - RG(k))' * (2𝕀 - RG(g * h)' - RG(k)) - ) + 2L' * L == - (2𝕀 - RG(g * h * k)' - RG(l))' * (2𝕀 - RG(g * h * k)' - RG(l)) - - @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + 2(2GH' * GH + 2K' * K) + 2L' * L == - (2𝕀 - RG(g * h * k)' - RG(l))' * (2𝕀 - RG(g * h * k)' - RG(l)) + - 2(2𝕀 - RG(g * h)' - RG(k))' * (2𝕀 - RG(g * h)' - RG(k)) - - @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + - 8G' * G + 8H' * H + 4K' * K + 2L' * L == - (2𝕀 - RG(g * h * k)' - RG(l))' * (2𝕀 - RG(g * h * k)' - RG(l)) + - 2(2𝕀 - RG(g * h)' - RG(k))' * (2𝕀 - RG(g * h)' - RG(k)) + - 4(2𝕀 - RG(g)' - RG(h))' * (2𝕀 - RG(g)' - RG(h)) - - @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + 2GH' * GH + 2KL' * KL == - (2𝕀 - RG(g * h)' - RG(k * l))' * (2𝕀 - RG(g * h)' - RG(k * l)) - - @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + - 2(2G' * G + 2H' * H) + - 2(2K' * K + 2L' * L) == - (2𝕀 - RG(g * h)' - RG(k * l))' * (2𝕀 - RG(g * h)' - RG(k * l)) + - 2(2𝕀 - RG(k)' - RG(l))' * (2𝕀 - RG(k)' - RG(l)) + - 2(2𝕀 - RG(g)' - RG(h))' * (2𝕀 - RG(g)' - RG(h)) - end + F = Groups.FreeGroup(4) + S = [Groups.gens(F); inv.(Groups.gens(F))] + + ID = one(F) + RADIUS = 3 + E_R, sizes = Groups.wlmetric_ball(S, ID; radius = 2 * RADIUS) + @test sizes == [9, 65, 457, 3201, 22409, 156865] + + b = SA.DiracBasis{UInt32}(F) + + RG = StarAlgebra(F, b) + + g, h, k, l = S[1:4] + + G = (one(RG) - RG(g)) + @test G' == one(RG) - RG(inv(g)) + @test G' * G == + MA.operate_to!(zero(G), *, G', G) == + 2one(RG) - RG(g) - RG(g)' + @test star(G * G) == G' * G' + + @testset "Sums of hermitian squares" begin + 𝕀 = one(RG) + + G = (𝕀 - RG(g)) + H = (𝕀 - RG(h)) + K = (𝕀 - RG(k)) + L = (𝕀 - RG(l)) + GH = (𝕀 - RG(g * h)) + KL = (𝕀 - RG(k * l)) + + X = (2𝕀 - star(RG(g)) - RG(h)) + Y = (2𝕀 - star(RG(g * h)) - RG(k)) + + @test -(2𝕀 - RG(g * h) - RG(g * h)') + 2G' * G + 2H' * H == X' * X + @test (2𝕀 - RG(g * h) - RG(g * h)') == GH' * GH + @test -(2𝕀 - RG(g * h * k) - RG(g * h * k)') + 2GH' * GH + 2K' * K == + Y' * Y + @test -(2𝕀 - RG(g * h * k) - RG(g * h * k)') + + 2(GH' * GH - 2G' * G - 2H' * H) + + 4G' * G + + 4H' * H + + 2K' * K == Y' * Y + + @test GH' * GH - 2G' * G - 2H' * H == -X' * X + @test -(2𝕀 - RG(g * h * k) - RG(g * h * k)') + + 4G' * G + + 4H' * H + + 2K' * K == 2X' * X + Y' * Y + + @test GH' * GH == + 2G' * G + 2H' * H - (2𝕀 - RG(g)' - RG(h))' * (2𝕀 - RG(g)' - RG(h)) + @test KL' * KL == + 2K' * K + 2L' * L - (2𝕀 - RG(k)' - RG(l))' * (2𝕀 - RG(k)' - RG(l)) + + @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + + 2 * GH' * GH + + 2 * KL' * KL == + (2𝕀 - RG(g * h)' - RG(k * l))' * (2𝕀 - RG(g * h)' - RG(k * l)) + + @test -(2𝕀 - star(RG(g * h * k * l)) - RG(g * h * k * l)) + + 2( + 2G' * G + 2H' * H - + (2𝕀 - RG(g)' - RG(h))' * (2𝕀 - RG(g)' - RG(h)) + ) + + 2( + 2K' * K + 2L' * L - + (2𝕀 - RG(k)' - RG(l))' * (2𝕀 - RG(k)' - RG(l)) + ) == + (2𝕀 - RG(g * h)' - RG(k * l))' * (2𝕀 - RG(g * h)' - RG(k * l)) + + @test -(2𝕀 - star(RG(g * h * k * l)) - RG(g * h * k * l)) + + 2(2G' * G + 2H' * H) + + 2(2K' * K + 2L' * L) == + (2𝕀 - RG(g * h)' - RG(k * l))' * (2𝕀 - RG(g * h)' - RG(k * l)) + + 2(2𝕀 - RG(g)' - RG(h))' * (2𝕀 - RG(g)' - RG(h)) + + 2(2𝕀 - RG(k)' - RG(l))' * (2𝕀 - RG(k)' - RG(l)) + + @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + + 2(2𝕀 - RG(g * h * k)' - RG(g * h * k)) + + 2L' * L == + (2𝕀 - RG(g * h * k)' - RG(l))' * (2𝕀 - RG(g * h * k)' - RG(l)) + + @test 2𝕀 - RG(g * h * k)' - RG(g * h * k) == + 2GH' * GH + 2K' * K - + (2𝕀 - star(RG(g * h)) - RG(k))' * (2𝕀 - star(RG(g * h)) - RG(k)) + + @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + + 2( + 2GH' * GH + 2K' * K - + (2𝕀 - RG(g * h)' - RG(k))' * (2𝕀 - RG(g * h)' - RG(k)) + ) + + 2L' * L == + (2𝕀 - RG(g * h * k)' - RG(l))' * (2𝕀 - RG(g * h * k)' - RG(l)) + + @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + + 2(2GH' * GH + 2K' * K) + + 2L' * L == + (2𝕀 - RG(g * h * k)' - RG(l))' * (2𝕀 - RG(g * h * k)' - RG(l)) + + 2(2𝕀 - RG(g * h)' - RG(k))' * (2𝕀 - RG(g * h)' - RG(k)) + + @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + + 8G' * G + + 8H' * H + + 4K' * K + + 2L' * L == + (2𝕀 - RG(g * h * k)' - RG(l))' * (2𝕀 - RG(g * h * k)' - RG(l)) + + 2(2𝕀 - RG(g * h)' - RG(k))' * (2𝕀 - RG(g * h)' - RG(k)) + + 4(2𝕀 - RG(g)' - RG(h))' * (2𝕀 - RG(g)' - RG(h)) + + @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + + 2GH' * GH + + 2KL' * KL == + (2𝕀 - RG(g * h)' - RG(k * l))' * (2𝕀 - RG(g * h)' - RG(k * l)) + + @test -(2𝕀 - RG(g * h * k * l)' - RG(g * h * k * l)) + + 2(2G' * G + 2H' * H) + + 2(2K' * K + 2L' * L) == + (2𝕀 - RG(g * h)' - RG(k * l))' * (2𝕀 - RG(g * h)' - RG(k * l)) + + 2(2𝕀 - RG(k)' - RG(l))' * (2𝕀 - RG(k)' - RG(l)) + + 2(2𝕀 - RG(g)' - RG(h))' * (2𝕀 - RG(g)' - RG(h)) + end end diff --git a/test/test_example_acoeffs.jl b/test/test_example_acoeffs.jl new file mode 100644 index 0000000..2c9272f --- /dev/null +++ b/test/test_example_acoeffs.jl @@ -0,0 +1,27 @@ +struct ACoeffs{T} <: SA.AbstractCoefficients{UInt32,T} + vals::Vector{T} +end + +# convienience: +ACoeffs(v::AbstractVector) = ACoeffs{valtype(v)}(v) + +# AbstractCoefficients API + +## Basic API +Base.keys(ac::ACoeffs) = (k for (k, v) in pairs(ac.vals) if !iszero(v)) +Base.values(ac::ACoeffs) = (v for v in ac.vals if !iszero(v)) +MA.operate!(::typeof(SA.canonical), ac::ACoeffs) = ac +function SA.star(b::SA.AbstractBasis, ac::ACoeffs) + return ACoeffs([ac.vals[star(b, k)] for k in keys(ac)]) +end + +# for preallocation in arithmetic +function Base.similar(ac::ACoeffs, ::Type{T}) where {T} + vals = similar(ac.vals, T) + MA.operate!(zero, vals) + return ACoeffs(vals) +end + +# the default arithmetic implementation uses this access +Base.getindex(ac::ACoeffs, idx) = ac.vals[idx] +Base.setindex!(ac::ACoeffs, val, idx) = ac.vals[idx] = val diff --git a/test/test_example_words.jl b/test/test_example_words.jl index 276a20e..5093366 100644 --- a/test/test_example_words.jl +++ b/test/test_example_words.jl @@ -17,7 +17,9 @@ function Base.show(io::IO, w::Word) end end -Base.:(==)(w::Word, v::Word) = w.alphabet == v.alphabet && w.letters == v.letters +function Base.:(==)(w::Word, v::Word) + return w.alphabet == v.alphabet && w.letters == v.letters +end Base.hash(w::Word, h::UInt) = hash(w.alphabet, hash(w.letters, hash(Word, h))) Base.one(w::Word) = Word(w.alphabet, Int[]) @@ -39,19 +41,37 @@ function StarAlgebras.star(w::Word) return Word(w.alphabet, newletters) end -function words(alphabet; radius) +Base.isless(w::Word, v::Word) = w.letters < v.letters + +struct FreeWords{T} + alphabet::Vector{T} +end - words = [Word(alphabet, Int[])] # word identity +Base.one(fw::FreeWords) = Word(fw.alphabet, Int[]) - for r in 1:radius - append!( - words, - [ - Word(alphabet, collect(w)) for - w in Iterators.product(fill(1:length(alphabet), r)...) - ], - ) +Base.eltype(::Type{FreeWords{T}}) where {T} = Word{T} +Base.IteratorSize(::Type{<:FreeWords}) = Base.IsInfinite() +function Base.iterate(aw::FreeWords) + w = Word(aw.alphabet, Int[]) + stack = [w] + return w, (stack, 1) +end + +function Base.iterate(aw::FreeWords, state) + stack, l = state + if l > length(aw.alphabet) + popfirst!(stack) + l = 1 end + w = first(stack) + nw = Word(aw.alphabet, [w.letters; l]) + push!(stack, nw) + return nw, (stack, l + 1) +end - return words +nwords(M::FreeWords, maxl::Integer) = nwords(M, 0, maxl) +function nwords(M::FreeWords, minl::Integer, maxl::Integer) + maxl < minl && return zero(maxl) + k = oftype(maxl, length(M.alphabet)) + return sum(k^i for i in minl:maxl) end