JuliaAlgebra / StarAlgebras.jl

A package for computation in *-algebras with basis
MIT License
5 stars 2 forks source link

Refactoring #20

Open blegat opened 8 months ago

blegat commented 8 months ago

Here are the conclusion of our meeting with @kalmarek We should have two types of AbstractBasis

abstract type AbstractBasis{T,I} end
abstract type AbstractSubBasis{T,I} <: AbstractBasis{T,I} end
struct SparseBasis{T,I,B<:AbstractBasis{T,I}}
    parent::B
    sub::Vector{I}
end

For AbstractSubBasis, you can build the subbasis with

basis(::AbstractVector{T}, basis::AbstractBasis{T,I}) # -> AbstractSubBasis{T,I}
algebra_element(::AbstractVector{C}, ::AbstractVector{T}, basis::AbstractBasis} # -> AlgebraElement
kalmarek commented 8 months ago

how about implementing our own SparseVector (without length)

struct SparseCoeffs{Tv, Tb, Vv<:AbstractVector{Tv}, Vb<:AbstractVector{Tb}}
   coeffs::Vv
   basis_elts::Vb
end

(without length!) and then on top of this we could have

struct AlgebraElement{A,SC<:SparseCoeffs}
    parent::A
    coeffs::SC
end

Also: multiplicative structures for implicit (e.g. potentially infinite) bases should work on the level of basis elements and not integers. For explicit bases (i.e. finite dimensional) translation from the former to the latter is already implemented as SA.Basis.

blegat commented 8 months ago

For implicit bases, we have two issues with SparseArrays.SparseVector: the indices may e.g., be monomials and not integer and the length is infinite. We can create a infinite sparse vector as follows and use it so that the same interface can be used:

struct InfiniteSparseVector{T,I}
    nzind::Vector{I}
    nzval::Vector{T}
end
Base.IteratorSize(::Type{<:InfiniteSparseVector}) = Base.IsInfinite()
import SparseArrays
SparseArrays.nonzeroinds(v::InfiniteSparseVector) = v.nzidx
SparseArrays.nonzeros(v::InfiniteSparseVector) = v.nzval

The sophisticated part of the printing of SparseVector actually don't need the length so we can reuse it as follows:

Base.show(io::IO, x::InfiniteSparseVector) = show(convert(IOContext, io), x)

function Base.show(io::IO, ::MIME"text/plain", x::InfiniteSparseVector)
    xnnz = length(nonzeros(x))
    print(io, typeof(x), " with ", xnnz,
           " stored ", xnnz == 1 ? "entry" : "entries")
    if xnnz != 0
        println(io, ":")
        show(IOContext(io, :typeinfo => eltype(x)), x)
    end
end

import SparseArrays
function Base.show(io::IOContext, x::InfiniteSparseVector)
    show(io, SparseVector(last(x.nzind), x.nzind, x.nzval))
end
blegat commented 8 months ago

To handle the monomials efficiently, the InfiniteSparseVector should use https://juliaalgebra.github.io/MultivariatePolynomials.jl/stable/types/#MultivariatePolynomials.monomial_vector_type and not a Vector{I} actually.

blegat commented 8 months ago

The coefficient vector can also be https://github.com/jump-dev/MathOptInterface.jl/blob/77d6b6fffc7cbc59dfa1696f1905a88162e34dc6/src/Utilities/set_dot.jl#L173-L212 We can redefined it as

struct CanonicalVector{I,T} <: AbstractVector{T} # for explicit
    index::I
    length::I
end
struct InfiniteCanonicalVector{K} # for implicit
    key::K
end
blegat commented 8 months ago

For decompose, you won't be able to call it on a SubBasis, it only works on a Basis. For a caching of the multiplication structure, if the indices are Base.OneTo(n) you don't hash and if it is like monomials, you hash. For monomials, you can have an implicit basis indexed by integer and one by monomials and you can easily try each and use the one faster for you.