Merge pull request #18 from kalmarek/enh/rewrite_mstructures

remove twisted from mstructures

Project.toml

@@ -1,7 +1,7 @@
 name = "StarAlgebras"
 uuid = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 authors = ["Marek Kaluba <[email protected]>"]
-version = "0.1.8"
+version = "0.2.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

README.md

@@ -15,19 +15,21 @@ julia> using StarAlgebras
 julia> using PermutationGroups
 
 julia> G = PermGroup(perm"(1,2)", perm"(1,2,3)")
-Full symmetric group over 3 elements
+Permutation group on 2 generators generated by
+ (1,2)
+ (1,2,3)
 
 julia> b = StarAlgebras.Basis{UInt8}(collect(G))
-6-element StarAlgebras.Basis{Permutation{...}, UInt8, Vector{...}}:
+6-element StarAlgebras.Basis{Permutation{Int64, …}, UInt8, Vector{Permutation{Int64, …}}}:
  ()
+ (2,3)
  (1,2)
  (1,3,2)
- (2,3)
- (1,2,3)
  (1,3)
+ (1,2,3)
 
- julia> RG = StarAlgebra(G, b)
- *-algebra of Permutation group on 2 generators of order 6
+julia> RG = StarAlgebra(G, b)
+*-algebra of Permutation group on 2 generators of order 6
 
 ```
 
@@ -66,6 +68,9 @@ julia> StarAlgebras.coeffs(f)
 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> f' # the same
+1·() -1·(2,3) +3·(1,2) +3·(1,3,2) -3·(1,3) -1·(1,2,3)
+
 julia> g = rand(G); g
 (1,2,3)
 
@@ -132,6 +137,24 @@ 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):
+```julia
+julia> x = RG(G(perm"(1,2,3)"))
+1·(1,2,3)
+
+julia> X = one(RG) - x
+1·() -1·(1,2,3)
+
+julia> X'
+1·() -1·(1,3,2)
+
+julia> X'*X
+2·() -1·(1,3,2) -1·(1,2,3)
+
+julia> @test X'*X == star(X)*X == 2one(X) - x - star(x)
+Test Passed
+
+```
 
 ### More advanced use
 
@@ -194,53 +217,30 @@ Finally, if the group is infinite (or just too large), but we need specific prod
 
 ### Even more advanced use (for experts only)
 
-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 `*`-algebra.
-
-To compute such sums of squares one may either sprinkle the code with `star`s, or define
+For low-level usage `MultiplicativeStructures` follow the sign convention:
 ```julia
-julia> tcmt = StarAlgebras.CachedMTable{true}(b, table_size=(length(b), length(b)));
+julia> mt = StarAlgebras.CachedMTable(b, table_size=(length(b), length(b)));
 
-```
-This multiplicative structure is **twisted** in the sense that `tcmt[i,j]` does not compute the product of `i`-th and `j`-th elements of the basis, but rather the `star` of `i`-th and `j`-th. An example
 ```julia
-julia> k = tcmt[i,j]
-0x03
+julia> k = mt[-i,j]
+0x06
 
 julia> @test star(b[i])*b[j] == b[k]
 Test Passed
 
 ```
-you should only use it with extreme care, as this "product" is no longer associative! Observe:
+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> @test tcmt[tcmt[3, 5], 4] == 0x06 # star(star(b[3])*b[5])*b[4] = star(b[5])*b[3]*b[4]
-Test Passed
-
-julia> @test tcmt[3, tcmt[5, 4]] == 0x01 # star(b[3])*star(b[5])*b[4]
+julia> @test mt[-signed(mt[-3, 5]), 4] == 0x06 # star(star(b[3])*b[5])*b[4] = star(b[5])*b[3]*b[4]
 Test Passed
-```
 
-However, writing sums of heritian squares is a breeze:
-```julia
-julia> tRG = StarAlgebra(G, b, tcmt)
-*-Algebra of Full symmetric group over 3 elements
-
-julia> x = tRG(G(perm"(1,2,3)"))
-1·(1,2,3)
-
-julia> X = one(tRG) - x
-1·() -1·(1,2,3)
-
-julia> @test X^2 == X*X == 2one(tRG) - x - star(x)
+julia> @test mt[-3, mt[-5, 4]] == 0x01 # star(b[3])*star(b[5])*b[4]
 Test Passed
-
 ```
 
-#### WARNING!
-
-Before using this mode you should consult the code (and potentially its author :) and understand very precisely what and where is happening!
 
 
 -----

src/algebra_elts.jl

@@ -30,19 +30,31 @@ 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
+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] = b[star(b[idx])]
+        idcs[i] = star(parent(X), idx)
         vals[i] = X[idx]
     end
     return AlgebraElement(sparsevec(idcs, vals, length(b)), A)
 end
 
+Base.adjoint(a::AlgebraElement) = star(a)
+
 LinearAlgebra.norm(a::AlgebraElement, p::Real) =
     LinearAlgebra.norm(coeffs(a), p)
 aug(a::AlgebraElement) = sum(coeffs(a))
@@ -54,5 +66,5 @@ 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(deepcopy(coeffs(a)), parent(a))
+    return AlgebraElement(Base.deepcopy_internal(coeffs(a), stackdict), parent(a))
 end

src/mstructures.jl

@@ -1,10 +1,4 @@
-abstract type MultiplicativeStructure{Twisted,I} <: AbstractMatrix{I} end
-
-_istwisted(::MultiplicativeStructure{T}) where {T} = T
-_product(ms::MultiplicativeStructure, g, h) = _product(Val(_istwisted(ms)), g, h)
-
-_product(::Val{false}, g, h) = g * h
-_product(::Val{true}, g, h) = star(g) * h
+abstract type MultiplicativeStructure{I} <: AbstractMatrix{I} end
 
 struct ProductNotDefined <: Exception
     i::Any
@@ -21,25 +15,26 @@ function Base.showerror(io::IO, ex::ProductNotDefined)
     print(io, ".")
 end
 
-struct TrivialMStructure{Tw,I,B<:AbstractBasis} <: MultiplicativeStructure{Tw,I}
+struct TrivialMStructure{I,B<:AbstractBasis} <: MultiplicativeStructure{I}
     basis::B
 end
 
-TrivialMStructure{Tw}(basis::AbstractBasis{T,I}) where {Tw,T,I} =
-    TrivialMStructure{Tw,I,typeof(basis)}(basis)
+TrivialMStructure(basis::AbstractBasis{T,I}) where {T,I} =
+    TrivialMStructure{I,typeof(basis)}(basis)
 
 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])])
 
 Base.@propagate_inbounds function Base.getindex(
     mstr::TrivialMStructure,
     i::Integer,
     j::Integer,
 )
-    @boundscheck checkbounds(mstr, i, j)
     b = basis(mstr)
+    i, j = _get(mstr, i), _get(mstr, j)
     g, h = b[i], b[j]
-    gh = _product(mstr, g, h)
+    gh = g * h
     gh in b || throw(ProductNotDefined(i, j, "$g · $h = $gh"))
-    return b[gh]
+    return basis(mstr)[gh]
 end