[WIP] GradedArraysNext

[lkdvos]

This is a large WIP of my progress for the GradedArraysNext refactor.
The main problem to simply start merging is that the testsuite needs either updating or marking things as broken, and I don't think I can find the time to get that done before tomorrow, so for now I'll just leave this as is.

[lkdvos]

@mtfishman I'm going to be honest, I'm kind of stuck here. I've trying to figure out what is going on and what codepaths the original implementation was going through, but I'm a bit at a loss. Any chance you have an idea what is going on, and what is going wrong? You can trigger an error here by running the test_tensoralgebraext.jl testset, the first matricize call hits this and I'm not entirely sure what is missing to make it work.

[mtfishman]

I can look into it, I guess this calls out to some of the BlockSparseArrays slicing code which is pretty complicated right now (we should be able to simplify it based on improvements to BlockArrays but I was waiting for this PR before doing that).

[mtfishman]

It's a bit involved, it looks like there is something a bit funny going on in the slicing code in BlockSparseArrays that was ok for dense blocks but causes issues for Kronecker/sector array blocks (at some level it erroneously switches over to doing some of these slicing operations elementwise instead of blockwise, it will take some investigating to figure out why it is doing that). For now I'm ok with matricize not sorting and merging sectors (and we can skip tests relying on that), I'll look into fixing this on the BlockSparseArrays side. Hopefully I can fix it for this PR, otherwise I'm fine with having it fixed in a followup PR.

[mtfishman]

It doesn't fully fix the issue, but it looks like it helps to overload matricize for SectorDelta:

using GradedArrays: SectorDelta, flip, tensor_product
using TensorAlgebra: TensorAlgebra as TA, trivialbiperm, blocks
struct SectorDeltaFusion <: TA.FusionStyle end
TA.FusionStyle(::Type{<:SectorDelta}) = SectorDeltaFusion()
function TA.matricize(style::SectorDeltaFusion, a::AbstractArray, ndims_codomain::Val)
    biperm = trivialbiperm(ndims_codomain, Val(ndims(a)))
    ax_codomain, ax_domain = blocks(axes(a)[biperm])
    ax_codomain = tensor_product(ax_codomain...)
    ax_domain = flip(tensor_product(ax_domain...))
    return SectorDelta{eltype(a)}((ax_codomain, ax_domain))
end

and probably we also need an overload of unmatricize.

[lkdvos]

Thanks for looking into this, I really was struggling and this does seem to (partially) resolve the issue. We can discuss this in a bit, but I am still not getting the required results (I think), which should also merge equal sectors, and I'm not sure what kind of indexing behavior would be used to get that to happen

[mtfishman]

No problem, it is definitely a bit hard to follow. The idea is that the block permutation and merging is done with these kinds of slice operations:

using BlockArrays: Block
using BlockSparseArrays: blocksparsezeros
a = blocksparsezeros(Float64, [2, 2, 2, 2], [2, 2, 2, 2])
for i in 1:4
    a[Block(i, i)] = randn(2, 2)
end
I = [Block(4), Block(3), Block(2), Block(1)]
a[I, I]
J = BlockedVector(I, [1, 2, 1])
a[J, J]

which gives:

julia> a
4×4-blocked 8×8 BlockSparseMatrix{Float64, Matrix{Float64}, …, …}:
 -0.338878  -0.407518  │    ⋅          ⋅        │   ⋅          ⋅       │    ⋅         ⋅      
  1.02568    1.25297   │    ⋅          ⋅        │   ⋅          ⋅       │    ⋅         ⋅      
 ──────────────────────┼────────────────────────┼──────────────────────┼─────────────────────
   ⋅          ⋅        │  -0.50356   -1.23022   │   ⋅          ⋅       │    ⋅         ⋅      
   ⋅          ⋅        │  -0.772989  -0.670701  │   ⋅          ⋅       │    ⋅         ⋅      
 ──────────────────────┼────────────────────────┼──────────────────────┼─────────────────────
   ⋅          ⋅        │    ⋅          ⋅        │  1.1336    -0.60612  │    ⋅         ⋅      
   ⋅          ⋅        │    ⋅          ⋅        │  0.321213  -1.15558  │    ⋅         ⋅      
 ──────────────────────┼────────────────────────┼──────────────────────┼─────────────────────
   ⋅          ⋅        │    ⋅          ⋅        │   ⋅          ⋅       │  -2.31704   0.937104
   ⋅          ⋅        │    ⋅          ⋅        │   ⋅          ⋅       │   1.09776  -0.30121 

julia> I = [Block(4), Block(3), Block(2), Block(1)]
4-element Vector{Block{1, Int64}}:
 Block(4)
 Block(3)
 Block(2)
 Block(1)

julia> a[I, I]
4×4-blocked 8×8 BlockSparseMatrix{Float64, Matrix{Float64}, …, …}:
 -2.31704   0.937104  │   ⋅          ⋅       │    ⋅          ⋅        │    ⋅          ⋅      
  1.09776  -0.30121   │   ⋅          ⋅       │    ⋅          ⋅        │    ⋅          ⋅      
 ─────────────────────┼──────────────────────┼────────────────────────┼──────────────────────
   ⋅         ⋅        │  1.1336    -0.60612  │    ⋅          ⋅        │    ⋅          ⋅      
   ⋅         ⋅        │  0.321213  -1.15558  │    ⋅          ⋅        │    ⋅          ⋅      
 ─────────────────────┼──────────────────────┼────────────────────────┼──────────────────────
   ⋅         ⋅        │   ⋅          ⋅       │  -0.50356   -1.23022   │    ⋅          ⋅      
   ⋅         ⋅        │   ⋅          ⋅       │  -0.772989  -0.670701  │    ⋅          ⋅      
 ─────────────────────┼──────────────────────┼────────────────────────┼──────────────────────
   ⋅         ⋅        │   ⋅          ⋅       │    ⋅          ⋅        │  -0.338878  -0.407518
   ⋅         ⋅        │   ⋅          ⋅       │    ⋅          ⋅        │   1.02568    1.25297 

julia> J = BlockedVector(I, [1, 2, 1])
3-blocked 4-element BlockedVector{Block{1, Int64}}:
 Block(4)
 ────────
 Block(3)
 Block(2)
 ────────
 Block(1)

julia> a[J, J]
3×3-blocked 8×8 BlockSparseMatrix{Float64, Matrix{Float64}, …, …}:
 -2.31704   0.937104  │   ⋅          ⋅         ⋅          ⋅        │    ⋅          ⋅      
  1.09776  -0.30121   │   ⋅          ⋅         ⋅          ⋅        │    ⋅          ⋅      
 ─────────────────────┼────────────────────────────────────────────┼──────────────────────
   ⋅         ⋅        │  1.1336    -0.60612   0.0        0.0       │    ⋅          ⋅      
   ⋅         ⋅        │  0.321213  -1.15558   0.0        0.0       │    ⋅          ⋅      
   ⋅         ⋅        │  0.0        0.0      -0.50356   -1.23022   │    ⋅          ⋅      
   ⋅         ⋅        │  0.0        0.0      -0.772989  -0.670701  │    ⋅          ⋅      
 ─────────────────────┼────────────────────────────────────────────┼──────────────────────
   ⋅         ⋅        │   ⋅          ⋅         ⋅          ⋅        │  -0.338878  -0.407518
   ⋅         ⋅        │   ⋅          ⋅         ⋅          ⋅        │   1.02568    1.25297 

so block permutations are performed by slicing with a list of blocks, and then if you block the block list that indicates which blocks should then be merged. You can see that right now those operations work for block sparse arrays with dense blocks but some step in the generic slicing code in BlockSparseArrays fails for Kronecker/sector blocks.

[lkdvos]

Ok that's really helpful (and I can now also see how it makes sense), let me see if I can try and figure out where things are going wrong, or where to intercept

[mtfishman]

Thanks, though again, if it turns out to be too tricky to figure out, I can look into it (I started to look into it and could see where things might be going wrong but it required some more careful investigating to figure out). I don't think this should necessarily hold up this PR since I think it is more of a BlockSparseArrays.jl issue. It involves some of the more complicated code in BlockSparseArrays.jl which I'm not so happy about anyway so I don't want to subject you too much to it unnecessarily. Also I'm happy to discuss/walk through the code that is being called for these kinds of slicing operations, it may help to talk through it since that may lead to some insights into how to simplify it (and in the process make it easier to debug/fix issues like this that pop up).

[codecov]

Codecov Report

❌ Patch coverage is 76.08696% with 55 lines in your changes missing coverage. Please review.
✅ Project coverage is 59.66%. Comparing base (b5565ba) to head (016d1d0).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/sectorarray.jl	70.68%	34 Missing ⚠️
src/tensoralgebra.jl	16.66%	10 Missing ⚠️
src/gradedarray.jl	85.18%	8 Missing ⚠️
src/sectorrange.jl	88.23%	2 Missing ⚠️
src/sectorproduct.jl	88.88%	1 Missing ⚠️

❗ There is a different number of reports uploaded between BASE (b5565ba) and HEAD (016d1d0). Click for more details.

HEAD has 2 uploads less than BASE

Flag	BASE (b5565ba)	HEAD (016d1d0)
docs	3	1

Additional details and impacted files

@@             Coverage Diff             @@
##             main      #79       +/-   ##
===========================================
- Coverage   93.64%   59.66%   -33.98%     
===========================================
  Files          13       10        -3     
  Lines         928      833       -95     
===========================================
- Hits          869      497      -372     
- Misses         59      336      +277     

Flag	Coverage Δ
docs	0.00% <0.00%> (ø)

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.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[lkdvos]

So, from what I can tell I now marked everything either as broken or fixed, this might be a good starting point for a review

[mtfishman]

So, from what I can tell I now marked everything either as broken or fixed, this might be a good starting point for a review

Great, thanks! I'll take a look.

[mtfishman]

gradedrange(start, [U1(1) => 2, U1(2) => 3, U1(3) => 2]) was syntax for adding an offset/start to the graded range (in the test here it is set to start at 1 so is kind of trivial but internally the type is still different, i.e. GradedUnitRange instead of GradedOneTo). It is analogous to the behavior of BlockArrays.blockedrange:

julia> using BlockArrays

julia> blockedrange([2, 3])
2-blocked 5-element BlockedOneTo{Int64, Vector{Int64}}:
 1
 2
 ─
 3
 4
 5

julia> blockedrange(1, [2, 3])
2-blocked 5-element BlockedUnitRange{Int64, Vector{Int64}}:
 1
 2
 ─
 3
 4
 5

julia> blockedrange(2, [2, 3])
2-blocked 5-element BlockedUnitRange{Int64, Vector{Int64}}:
 2
 3
 ─
 4
 5
 6

But, we can leave that syntax out for now and just use broadcasting syntax to add offsets, so gradedrange(1, [U1(1) => 2, U1(2) => 3, U1(3) => 2]) would be equivalent to gradedrange([U1(1) => 2, U1(2) => 3, U1(3) => 2]) .+ 0.

[mtfishman]

@lkdvos thanks for going through all of this and getting it working. This looks like a good starting point for future work, ready to merge from your end?

[lkdvos]

Yes, I don't think I have anything immediate to add!