changes for SymbolicWedderburn

Project.toml

@@ -9,6 +9,6 @@ MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-MutableArithmetics = "1"
+MutableArithmetics = "1.4.3"
 SparseArrays = "1"
 julia = "1.6"

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,98 +44,112 @@ 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
 
 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.

src/StarAlgebras.jl

@@ -6,13 +6,12 @@ import LinearAlgebra
 import MutableArithmetics as MA
 
 export StarAlgebra, AlgebraElement
-export aug, basis, coeffs, star, supp
+export basis, coeffs, star
 
 # AbstractCoefficients
 ## abstract definitions
 include("coefficients.jl")
-## concrete implementations
-include("diracs_augmented.jl")
+## concrete implementation
 include("sparse_coeffs.jl")
 
 # MultiplicativeStructures
@@ -34,4 +33,7 @@ include("star.jl")
 include("arithmetic.jl")
 include("show.jl")
 
+# augmented basis implementation
+include("diracs_augmented.jl")
+
 end

src/algebra_elts.jl

@@ -15,8 +15,8 @@ function Base.deepcopy_internal(a::AlgebraElement, id::IdDict)
 end
 
 # call overload:
-(a::AlgebraElement)(x) = coeffs(a)[x]
-Base.setindex!(a::AlgebraElement, v, idx) = a.coeffs[idx] = v
+(a::AlgebraElement)(x) = coeffs(a)[basis(parent(a))[x]]
+Base.setindex!(a::AlgebraElement, v, idx) = a.coeffs[basis(parent(a))[idx]] = v
 
 # AlgebraElement specific functions
 
@@ -25,8 +25,6 @@ function supp(a::AlgebraElement)
     return [b[i] for (i, _) in nonzero_pairs(coeffs(a))]
 end
 
-aug(a::AlgebraElement) = aug(coeffs(a))
-
 function LinearAlgebra.norm(a::AlgebraElement, p::Real)
     return LinearAlgebra.norm(coeffs(a), p)
 end