[NDTensors] replace UnitRangeDual with GradedUnitRangeDual

[ogauthe]

This PR aims to simplify the use of GradedAxes.dual, especially to be used in a BlockSparseArray. It proposes to remove the UnitRangeDual type:

• a generic AbstractUnitRange is now self dual
• a labelled GradedUnitRange or GradedOneTo are cast to GradedUnitRangeDual{<:LabelledInteger,<:AbstractGradedUnitRange} <: AbstractGradedUnitRange

This should simplify slicing a BlockSparseArray with dual axes and allows to preserve labels or take their dual.

This branch is similar to #1529 but is based on main.

Currently, this branche solves:

• preserving labels in blocklengths(dual(...))

g = gradedrange([U1(1) => 2, U1(2) => 2])
@test GradedAxes.blocklengths(dual(g)) isa Vector{<:LabelledInteger}  #    # broken in main and in non-abelian_fusion

• preserving GradedUnitRange behavior when slicing:

blocklengths(dual(gradedrange([U1(0) => 2, U1(1) => 2])[1:3])) isa Vector  # broken in main and in non-abelian_fusion

• manipulating arrays with DualGradedUnitRange

r = gradedrange([U1(0) => 2, U1(1) => 2])[1:3]
a = BlockSparseArray{Float64}(r, r)
b = BlockSparseArray{Float64}(dual(r),dual(r))
@test a[:, :] isa BlockSparseArray   # broken in main and in non-abelian_fusion
@test b[:, :] isa BlockSparseArray   # broken in main and in non-abelian_fusion

Meaning one can create and manipulate BlockSparseArray with GradedUnitRange or its dual (not just BlockedOneTo). However axes(a[:,:]) are still Base.IdentityUnitRange(GradedUnitRangeDual{...}), which creates other issues.

It also adds many tests, including broken ones.

[ogauthe]

Overall I still find there are too many kind of UnitRange. I wonder if we could avoid the doublet GradedOneTo and GradedUnitRange. gradedrange creates a GradedOneTo, but when slicing it one gets a GradedUnitRange.

For instance, blockfirsts(gradedrange(["x" => 2, "y" => 3])[1:5]) currently fails due to ambiguity.

[mtfishman]

For dispatch it is helpful to have special ranges that start at one since that is used as axes by most arrays, like how Base has UnitRange and Base.OneTo. Handling both shouldn't be too difficult if we write good generic code.

[mtfishman]

To be more specific, if we write most code in terms of AbstractGradedUnitRange and AbstractUnitRange, that should allow us to have a lot of subtypes (such as GradedOneTo, GradedUnitRage, GradedUnitRangeDual, etc.) without adding too much code complexity. Though perhaps one question is, is there a need for a GradedOneToDual? That does start to sound excessive and it is best to try to design it without that for now.

One thing to consider is that dual only really makes sense for graded unit ranges, so maybe we could consider dropping UnitRangeDual, and just define dual(r::AbstractUnitRange) = r as a generic fallback (i.e. if the unit range isn't graded, the dual is itself).

[ogauthe]

I have no strong opinion whether there should be a UnitRangeDual type or not. This is mostly a design question for BlockSparseArray, about how much complexity you want there. In the specific case of FusionTensor, I will always use a graded axis. My objective here is to make dual(GradedUnitRange) an AbstractGradedUnitRange to deal with labels in a generic way.

However I believe dual and not dual axes should be handled on equal footing: deciding who is dual and who is not is a pure convention and they should always appear with a 1:1 ratio for contraction to be possible. If there is some gain in having both GradedOneTo and GradedUnitRangeDual, then I think there should be a GradedOneToDual.

[ogauthe]

I realize that with the current design, there are already GradedUnitRangeDual{<:LabelledInteger, BlockArrays.BlockedOneTo and GradedUnitRangeDual{<:LabelledInteger, BlockArrays.BlockedUnitRange}, so there should be no need for an explicit GradedOneToDual.

I will remove UnitRangeDual and use self-dual axes as the default.

[mtfishman]

I realize that with the current design, there are already GradedUnitRangeDual{<:LabelledInteger, BlockArrays.BlockedOneTo and GradedUnitRangeDual{<:LabelledInteger, BlockArrays.BlockedUnitRange}, so there should be no need for an explicit GradedOneToDual.

I will remove UnitRangeDual and use self-dual axes as the default.

Sounds good, that's what I was going to suggest (i.e. GradedOneToDual can be a type alias, if needed).

[ogauthe]

I think this is now ready for GradedAxes. I will focus on BlockSparseArrays in another PR.

[mtfishman]

@ogauthe I realized that one reason to keep UnitRangeDual is that it may be a reasonable return type when slicing blocks of a GradedUnitRangeDual, i.e. it might make sense for dual(gradedrange([U1(0) => 2, U1(1) => 2]))[Block(1)] to return a UnitRangeDual, to preserve the fact that the block was taken from a dual space.

That could be helpful since then it could make it easier to create a dual graded unit range out of the blocks of another dual graded unit range (say if you have a dual graded unit range and want to remove one of the blocks, it could make it easier to write that kind of code).

[ogauthe]

Indeed slicing a GradedUnitRangeDual with a Block was an issue I spotted (this is associated with the slicing issues in #1336). Currently it yields a LabelledUnitRange and the dual information is lost.

Now, how generic UnitRangeDual should be? I guess it should be a subtype of AbstractUnitRange and should be valid axis for a BlockSparseArray. Should we come back to "any AbstractUnitRange can be dualed and become a UnitRangeDual" or should it be more specific, as "obtained from slicing a GradedUnitRangeDual"?

I favor the second as it would make it simpler to use a BlockSparseArray without caring about mathematical details, but I think both options carry the same amount of code complexity (mostly coming from being valid BlockSparseArray axes).

[mtfishman]

I agree with your assessment that UnitRangeDual should be a subtype of AbstractUnitRange, and also that we should be thoughtful about which operations create UnitRangeDual. I'll try to summarize that view:

I think it should be valid to make a BlockSparseArray with any kind of unit range as long as it satisfies the BlockArrays.jl axis interface, and maybe some additional interface functions we require in BlockSparseArrays, so axes may be dual or non-dual axes of whatever kind.

The question to me more comes down to how and when dual axes get constructed, depending on the unit range type and operation. So to be concrete, if you have a block sparse matrix a, in general a' should have axes: axes(a') == dual.(reverse(axes(a))), so that operations like a' * a make sense (i.e. involve a contraction of a dual axis with a non-dual axis).

However, I think it would be surprising and unnecessarily complicated if a started with non-graded axes, say OneTo or BlockOneTo, and a' had dual axes, say UnitRangeDual{<:Any,<:OneTo} or UnitRangeDual{<:Any,<:BlockOneTo}, since they are self-dual. So that leads me to think that we should define dual(r::OneTo) = r and dual(r::BlockOneTo) = r (since adjoint(::BlockSparseMatrix) will be defined in terms of dual to take the dual of the axes where appropriate), and if a user wants to create a BlockSparseArray with dual axes of type UnitRangeDual{<:Any,<:BlockOneTo}, they'll have to do that manually in another way besides the dual(...) function.

We may want to define that dual(r::LabelledUnitRange) outputs a UnitRangeDual(r), though in that case it may be better to rename LabelledUnitRange to something else to indicate that it is more specifically the block of a GradedUnitRange.