Implement SubBasis

src/bases.jl

@@ -24,7 +24,6 @@ still be infinite).
   [`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
 
@@ -34,8 +33,11 @@ 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
+
+Implicit bases are bases that contains the product of all its elements.
+This makes these bases particularly useful to work with [`AlgebraElement`](@ref)s with supports that can not be reasonably bounded.
+Note that these bases may not explictly store its elements in memory as they may be potentially infinite.
 """
 abstract type ImplicitBasis{T,I} <: AbstractBasis{T,I} end
 
@@ -44,21 +46,22 @@ function zero_coeffs(::Type{S}, ::ImplicitBasis{T,I}) where {S,T,I}
 end
 
 """
-    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<:Integer} <: AbstractBasis{T,I} end
+
+Explicit bases are bases of finite length for which the keys are integers.
 """
-abstract type ExplicitBasis{T,I} <: AbstractBasis{T,I} end
+abstract type ExplicitBasis{T,I<:Integer} <: AbstractBasis{T,I} end
 
+Base.IteratorSize(::Type{<:ExplicitBasis}) = Base.HasLength()
 Base.keys(eb::ExplicitBasis) = Base.OneTo(length(eb))
 Base.size(eb::ExplicitBasis) = (length(eb),)
 
 function zero_coeffs(::Type{S}, eb::ExplicitBasis{T,I}) where {S,T,I}
     return spzeros(S, I, length(eb))
 end
 
-function Base.getindex(eb::ExplicitBasis, range::AbstractRange{<:Integer})
-    return [eb[i] for i in range]
+function Base.getindex(eb::ExplicitBasis{T}, range::AbstractRange{<:Integer}) where {T}
+    return T[eb[i] for i in range]
 end
 
 """

src/bases_fixed.jl

@@ -8,14 +8,66 @@ mutable struct FixedBasis{T,I,V<:AbstractVector{T}} <:
     starof::Vector{I}
 end
 
+"""
+    FixedBasis(elts::AbstractVector)
+
+Represents a linear basis which is fixed of length (given by `elts`).
+The elements of `elts` are assumed to be _linearly independent_ and (as a set) invariant under `star`.
+Due to its finiteness, `b::FixedBasis` can be indexed like a vector (by positive integers) and therefore elements expressed in `b` can be represented by simple (sparse) Vectors.
+
+## Examples
+
+```julia
+julia> import StarAlgebras as SA;
+
+julia> SA.star(s::String) = reverse(s) # identity is also fine
+
+julia> b = SA.FixedBasis(["", "a", "b", "aa", "ab", "ba", "bb"]);
+
+julia> b["a"]
+1
+
+julia> b["ab"]
+4
+
+julia> b[4]
+"ab"
+
+julia> b["abc"]
+ERROR: KeyError: key "abc" not found
+[...]
+
+julia> A = SA.StarAlgebra(String, SA.DiracMStructure(b, *))
+*-algebra of String
+
+julia> c = A("a") + A("b")
+1·"a" + 1·"b"
+
+julia> c*A("a")
+1·"aa" + 1·"ba"
+
+julia> A("a")*c
+1·"aa" + 1·"ab"
+
+julia> c*c
+1·"aa" + 1·"ab" + 1·"ba" + 1·"bb"
+
+julia> SA.coeffs(c*c)
+7-element SparseArrays.SparseVector{Int64, Int64} with 4 stored entries:
+  [4]  =  1
+  [5]  =  1
+  [6]  =  1
+  [7]  =  1
+```
+"""
+FixedBasis(elts::AbstractVector{T}) where {T} = FixedBasis{T,keytype(elts)}(elts)
+
 function FixedBasis{T,I}(elts::AbstractVector{T}) where {T,I}
     relts = Dict(b => I(idx) for (idx, b) in pairs(elts))
     starof = [relts[star(x)] for x in elts]
     return FixedBasis{T,I,typeof(elts)}(elts, relts, starof)
 end
 
-FixedBasis(elts::AbstractVector{T}) where {T} = FixedBasis{T,keytype(elts)}(elts)
-
 function FixedBasis{T,I}(basis::AbstractBasis{T}; n::Integer) where {T,I}
     return FixedBasis{T,I}(collect(Iterators.take(basis, n)))
 end
@@ -26,15 +78,71 @@ Base.in(x, b::FixedBasis) = haskey(b.relts, x)
 Base.getindex(b::FixedBasis{T}, x::T) where {T} = b.relts[x]
 Base.getindex(b::FixedBasis, i::Integer) = b.elts[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]
+"""
+    struct SubBasis{T,I,K,B<:AbstractBasis{T,K},V<:AbstractVector{K}} <:
+        ExplicitBasis{T,I}
+        parent_basis::B
+        keys::V
+        is_sorted::Bool
+    end
+
+Represents a sub-basis of a given basis, where `keys` is a vector of keys
+representing the sub-basis elements, `parent_basis` is the parent basis from
+which the sub-basis is derived, and `is_sorted` indicates whether the keys are
+sorted.
+"""
+struct SubBasis{T,I,K,B<:AbstractBasis{T,K},V<:AbstractVector{K}} <:
+       ExplicitBasis{T,I}
+    parent_basis::B
+    keys::V
+    is_sorted::Bool
+    function SubBasis(parent_basis::AbstractBasis{T,K}, keys::AbstractVector{K}) where {T,K}
+        return new{T,keytype(keys),K,typeof(parent_basis),typeof(keys)}(parent_basis, keys, issorted(keys))
+    end
+end
+
+Base.parent(sub::SubBasis) = sub.parent_basis
+
+Base.length(b::SubBasis) = length(b.keys)
+function _iterate(b::SubBasis, elem_state)
+    if isnothing(elem_state)
+        return
+    end
+    elem, state = elem_state
+    return parent(b)[elem], state
+end
+Base.iterate(b::SubBasis) = _iterate(b, iterate(b.keys))
+Base.iterate(b::SubBasis, st) = _iterate(b, iterate(b.keys, st))
+
+function Base.get(b::SubBasis{T,I}, x::T, default) where {T,I}
+    key = b.parent_basis[x]
+    if b.is_sorted
+        i = searchsortedfirst(b.keys, key)
+        if i in eachindex(b.keys) && b.keys[i] == key
+            return convert(I, i)
+        end
+    else
+        i = findfirst(isequal(key), b.keys)
+        if !isnothing(i)
+            return convert(I, i)
+        end
+    end
+    return default
+end
+
+Base.in(x::T, b::SubBasis{T}) where T = !isnothing(get(b, x, nothing))
+Base.haskey(b::SubBasis, i::Integer) = i in eachindex(b.keys)
+
+Base.getindex(b::SubBasis, i::Integer) = parent(b)[b.keys[i]]
+function Base.getindex(b::SubBasis{T,I}, x::T) where {T,I}
+    i = get(b, x, nothing)
+    if isnothing(i)
+        throw(KeyError(x))
+    end
+    return i::I
+end

test/perm_grp_algebra.jl

@@ -1,6 +1,10 @@
 # This file is a part of StarAlgebras.jl. License is MIT: https://github.com/JuliaAlgebra/StarAlgebras.jl/blob/main/LICENSE
 # Copyright (c) 2021-2025: Marek Kaluba, Benoît Legat
 
+using Test
+using PermutationGroups
+import Random
+import StarAlgebras as SA
 @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)
@@ -99,4 +103,42 @@
         end
         @test a ≤ b
     end
+
+    @testset "SubBasis" begin
+        # If we're unlucky, `a^3` might belong to the basis.
+        # We fix the seed to be sure we are never in that case.
+        Random.seed!(0)
+        S1 = unique!(rand(G, 7))
+        S = unique!([S1; [a * b for a in S1 for b in S1]])
+        subb = SA.SubBasis(db, S)
+        a = S1[1]
+        @test subb[a] == 1
+        @test a in subb
+        @test isnothing(get(subb, a^3, nothing))
+        @test_throws KeyError(a^3) subb[a^3]
+        @test !(a^3 in subb)
+        @test collect(subb) == S
+        smstr = SA.DiracMStructure(subb, *)
+        @test only(smstr(1, 2).basis_elements) == subb[subb[1] * subb[2]]
+        @test only(smstr(1, 2, eltype(subb)).basis_elements) == subb[1] * subb[2]
+
+        sbRG = SA.StarAlgebra(G, subb)
+
+        x = let z = zeros(Int, length(SA.basis(sbRG)))
+            z[1:length(S1)] .= rand(-1:1, length(S1))
+            SA.AlgebraElement(z, sbRG)
+        end
+
+        y = let z = zeros(Int, length(SA.basis(sbRG)))
+            z[1:length(S1)] .= rand(-1:1, length(S1))
+            SA.AlgebraElement(z, sbRG)
+        end
+
+        dx = SA.AlgebraElement(SA.coeffs(x, SA.basis(RG)), RG)
+        dy = SA.AlgebraElement(SA.coeffs(y, SA.basis(RG)), RG)
+
+        @test dx + dy == SA.AlgebraElement(SA.coeffs(x + y, SA.basis(RG)), RG)
+
+        @test dx * dy == SA.AlgebraElement(SA.coeffs(x * y, SA.basis(RG)), RG)
+    end
 end

test/quadratic_form.jl

@@ -74,19 +74,6 @@
     @test A(Q) == π * b[3] * b[2] + b[4] * b[3]
 end
 
-struct SubBasis{T,I,V<:AbstractVector{I},B<:SA.ImplicitBasis{T,I}} <: SA.ExplicitBasis{Float64,Int}
-    implicit::B
-    indices::V
-end
-Base.length(b::SubBasis) = length(b.indices)
-function Base.iterate(b::SubBasis, args...)
-    elem_state = iterate(b.indices, args...)
-    if isnothing(elem_state)
-        return
-    end
-    return b.implicit[elem_state[1]], elem_state[2]
-end
-
 @testset "Int -> Float basis" begin
     limited = SA.MappedBasis(1:2, float, Int)
     @test !(3.0 in limited)
@@ -95,7 +82,10 @@ end
     implicit = SA.MappedBasis(NaturalNumbers(), float, Int)
     @test 3.0 in implicit
     @test haskey(implicit, 3)
-    explicit = SubBasis(implicit, 1:3)
+    explicit = SA.SubBasis(implicit, 1:3)
+    @test 3.0 in explicit
+    @test collect(explicit) == [1.0, 2.0, 3.0]
+    @test haskey(explicit, 3)
     m = Bool[
         true  false true
         false true  false
@@ -119,7 +109,7 @@ end
     implicit = cheby_basis()
     mstr = ChebyMStruct(implicit)
     mt = SA.MTable(implicit, mstr, (0, 0))
-    sub = SubBasis(implicit, 1:3)
+    sub = SA.SubBasis(implicit, 1:3)
     fixed = SA.FixedBasis(implicit; n = 3)
     a = ChebyPoly(2)
     b = ChebyPoly(3)