Non monomial basis

[kalmarek]

This is not even a working code, but a place where we can discuss more concretely

What still needs to be implemented:

• arithmetic on coefficient vectors
• decomposition in bases
• how to encode whether MStructure is of "Dirac type"? (traits?)
• how to interact with Implicit basis integer indexing seems to be hard in general -- maybe only optional? (again hidden behind a trait? IndexType?)
• … (and many more)

@blegat

[blegat]

What's failing now are the MTables change which I'm less confident with so I hand over the reins back to you

[kalmarek]

@blegat I'm not trying to get tests running now. (given the number of planned changes it is a bit pointless., I'm trying to get the implicit/explicit Bases running on my examples;

e.g. with

G = PermGroup(perm"(1,2,3,4,5,6,7,8)", perm"(1,2)")
db = SA.DiracBasis{UInt32}(G) # <: ImplicitBasis{UInt32, eltype(G)} for an iterable G
mstr = SA.LazyMStructure(db)
g, h = rand(G, 2);

and

function Base.show(io::IO, sc::SA.AbstractCoefficients)
    if iszero(sc)
        print(io, zero(valtype(sc)))
    else
        for (i, (k, v)) in enumerate(pairs(sc))
            print(io, v, '·', k)
            i == length(pairs(sc)) && break
            print(io, " + ")
        end
    end
end

I can do now:

julia> mstr[g, h] == SA.DiracDelta(g * h)
true

julia> x = SA.SparseCoefficients([one(G), g], [1, -1])
1·() + -1·(1,6,5,7,2,8,3,4)

julia> y = sc = SA.SparseCoefficients([one(G), inv(g)], [1, -1])
1·() + -1·(1,4,3,8,2,7,5,6)

julia> res = zero(x)
0

julia> SA.mul!(mstr, res, x, y); res
┌ Error: in __canonicalize!: Not implemented yet
└ @ StarAlgebras ~/.julia/dev/StarAlgebras/src/coefficients.jl:71
1·() + -1·(1,4,3,8,2,7,5,6) + -1·(1,6,5,7,2,8,3,4) + 1·()

🎉 Looks good to me even without implementing integer access to DiracBasis!

What I don't quite like at the moment is that I had to implement a specialized

Base.getindex(::LazyMStructure{I,<:DiracBasis{T}}, x::T, y::T) where {I,T}

It doesn't feel right...

---

Also: The concept of I as indextype for SA.MultiplicativeStructure{I} maybe should be tossed out.

[blegat]

It's working with chebyshev polynomials now, see JuliaAlgebra/MultivariateBases.jl#23
The only issues are

• I kind of worked around the mtables by directly implementing fmac! for the AblebraElement. I need to think about it more
• Didn't do canonicalize!. It's actually nontrivial to do this inplace. If you call sortperm on the basis_elements then you still need to reorder inplace which isn't so easy. It's best to call sort! instead but here we have two separate vectors so if we sort one they won't stay in sync. What I do in MultivariatePolynomials is to call sort! on a vector of terms. Here, We should create some sort of AbstractVector{<:Pair} that we give to sort! and for each operation used by sort!, we update both vectors, which is what https://github.com/JuliaArrays/StructArrays.jl is doing. They even have some sorting things: https://github.com/JuliaArrays/StructArrays.jl/blob/master/src/sort.jl

[kalmarek]

@blegat Sure, we could switch to storing pairs (key => value) along each other and the whole interface to using pairs.

[kalmarek]

I got both implicit (DiracBasis) and explicit (FixedBasis) to work;

with some preparation

using Revise
using StarAlgebras
import StarAlgebras as SA
using PermutationGroups
import PermutationGroups.AP as AP

SA.star(x::Number) = x'
SA.star(g::AP.AbstractPermutation) = inv(g)

struct Lex <: Base.Order.Ordering end
function Base.lt(::Lex, p::AP.AbstractPermutation, q::AP.AbstractPermutation)
    for i in Base.OneTo(min(AP.degree(p), AP.degree(q)))
        ip = i^p
        iq = i^q
        if ip < iq
            return true
        elseif ip > iq
            return false
        end
    end
    return AP.degree(p) < AP.degree(q)
end
Base.isless(p::AP.AbstractPermutation, q::AP.AbstractPermutation) =
    Base.lt(Lex(), p, q)

G = PermGroup(perm"(1,2,3,4,5,6,7,8)", perm"(1,2)") # group of order 40320

one can run now:

db = SA.DiracBasis{UInt32}(G)
RG = SA.StarAlgebra(G, db)

g, h = rand(G, 2)

x = SA.SparseCoefficients([one(G), g], [1, -1])
X = SA.AlgebraElement(x, RG)

y = SA.SparseCoefficients([one(G), inv(g)], [1, -1])
Y = SA.AlgebraElement(y, RG)

xy = SA.SparseCoefficients([one(G), g, inv(g)], [2, -1, -1])
XY = SA.AlgebraElement(xy, RG)

@assert X != Y
@assert X' == Y
@assert RG.mstructure[g, h] == SA.DiracDelta(g * h)
@assert X*Y == XY

But also this:

RGf = let G = G, g = g, h = h
    fb = let (g, h) = (g, h)
        # somehow generate enough elements of G to have
        # a non-trivial multiplication table
        for a in (g, h, inv(g), inv(h), one(g))
            for b in (g, h, inv(g), inv(h), one(g))
                SA.__cache(db, a * b)
            end
        end
        # elts must be star-invariant
        elts = union(db.basis, SA.star.(db.basis))
        SA.FixedBasis{UInt32}(elts)
    end

    k = max(fb[g], fb[h], fb[inv(g)], fb[inv(h)])
    mtbl = SA.MTable(fb, size=(k, k))
    StarAlgebra(G, fb, mtbl)
end

let fb = SA.basis(RGf), g = g
    k = fb[g]
    @assert isone(fb[k] * fb[-signed(k)])
end

using SparseArrays
xf = sparsevec([fb[one(G)], fb[g]], [1, -1], length(fb))
Xf = AlgebraElement(xf, RGf)

yf = sparsevec([fb[one(G)], fb[inv(g)]], [1, -1], length(fb))

Yf = AlgebraElement(yf, RGf)

@assert Xf == Yf'
@assert (Xf')' == Xf
@assert Xf * Xf' == Xf' * Xf

[kalmarek]

Ok, more input/experimentation here. The following fits my usecase (where AlgebraElements have their parent objects, but maybe we could tweak it to the Polynomial case as well?)

f = # polynomial, AlgebraElement, etc, possibly written in different basis
g = # polynomial, AlgebraElement, etc, possibly written in different basis

function +(f, g)
    A = parent(f)
    res = coeffs(f) + coeffs(g, basis(A))
    return AlgebraElement(res, A)
end

function *(f, g)
    A = parent(f)
    mstr = mstructure(A)
    res = mul(mstr, coeffs(f), coeffs(g, basis(A)))
    return AlgebraElement(res, A)
end

To spell out assumptions:

• objects f have their parent object which comes with Explicit/Implicit basis
• coeffs(g, basis(A))::AbstractCoefficients decomposes coeffs(g) into (possibly sparse) AbstractCoefficients written in basis of A.
• there is some basic arithmetic on AbstractCoefficients over the same basis (but without checking this assumption, by design)
• mul(mstr, cf, cg) returns AbstractCoefficients which will depend very much on what cf, cg and mstr are. The assumption is (again not checked) that cf and cg are already written in the basis over which mstr is defined.

@blegat I've implemented both implicit, fixed, and augmented (consisting of (1-x)) basis and it seems to be a valid design. what do you think?

[blegat]

What if mul(mstr, cf, cg) becomes mstr(cf, cg), then we'll also have operate!(mstr, cf, cg) etc... and it would play nicely with MutableArithmetics and MultivariatePolynomials

[blegat]

function *(f, g)
    A = parent(f)
    mop = moperator(A)
    res = mop(coeffs(f), coeffs(g, basis(A)))
    return AlgebraElement(res, A)
end

[blegat]

Suggested change

	mutable struct DiracBasis{T,I,S,St} <: ImplicitBasis{T,I}
	mutable struct Basis{T,I,S,St} <: ImplicitBasis{T,I}

[blegat]

Add a field which is * by default and contains the multiplicative structure

[blegat]

mstructure(::T, ::T) should be defined

[kalmarek]

Basis(itr, moperator=LazyMStructure(*))

(lms::LazyMStructure)(x,y) = lms.op(x,y)::?

lms.op(x,y)::AbstractCoefficients

[blegat]

Basis(itr) = Basis(itr, *)

[kalmarek]

things considered as AbstractCoefficients:

MPoly
AugmentedDirac (1-x)
SparseVector # fixed length
SparseCoefficients # arbitrary length

[blegat]

struct AddMul{M}
    mstr::M
end
function MA.operate!(::AddMul{Cheby}, accumulator, a, b) # the new unsafe_append!
    MA.operate!(MA.add_mul, accumulator, 1/2, ...) # e.g. in Chebyshev basis we add two terms
    MA.operate!(MA.add_mul, accumulator, -1/2, ...)
    return result
end
struct BuildContext{C}
    inner::C
end
function MA.operate!(::MA.add_mul, accumulator::BuildContext, a, b) # the new unsafe_append!
    MA.operate!(UnsafeAdd(), accumulator.inner, a, b)
    return result
end
function MA.operate!(::UnsafeAdd, coef::SparseCoefficients, a, b) # the new unsafe_append!
    push!(coef.values, a)
    push!(coef.basis_elements, b)
    return coef
end

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 93.79433% with 35 lines in your changes are missing coverage. Please review.

Project coverage is 93.51%. Comparing base (f1b98a0) to head (fc2e9b2).
Report is 4 commits behind head on main.

Files	Patch %	Lines
src/coefficients.jl	87.38%	14 Missing ⚠️
src/mtables.jl	91.25%	7 Missing ⚠️
src/diracs_augmented.jl	91.78%	6 Missing ⚠️
src/mstructures.jl	86.66%	4 Missing ⚠️
src/sparse_coeffs.jl	97.59%	2 Missing ⚠️
src/bases.jl	95.65%	1 Missing ⚠️
src/bases_fixed.jl	95.23%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #21      +/-   ##
==========================================
- Coverage   98.60%   93.51%   -5.10%     
==========================================
  Files           8       14       +6     
  Lines         358      632     +274     
==========================================
+ Hits          353      591     +238     
- Misses          5       41      +36     

Flag	Coverage Δ
unittests	93.51% <93.79%> (-5.10%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[kalmarek]

I'm reasonably happy with this coverage.
There are 6 more broken tests @blegat, please have a look

[kalmarek]

🎉