Make better use of BlockedTuple in contract logic to track codomain and domain

[ogauthe]

solves #22
see also ITensor/BlockSparseArrays.jl#55

This is a draft to systematically preserve domain and codomain information in contract. This PR makes several changes in contract, fusedims and splitdims. It preserves the current high level interface, adds new interface using BlockedTuple and AbstractBlockPermutation, and rewrite internals using these AbstractBlockTuple types.

I tried to relax BlockedPermutation to AbstractBlockPermutation wherever I could. We may consider more lax signatures using BlockedTuple{Parameters}.

In the details:

• contract imposes AbstractBlockPermutation{2} at lower levels of contract. It implies that final contraction is always done between matrices, adding trailing 1 dimension if needed. In particular, it means dispatch on _mul! to handle specific cases with 0 or 1 dim arrays are no more needed. The code is simpler, but some specificities of low dimensions are lost. This is how NumPy handles matmul, so I think this should not affect performances much, although the devil is in the details and I may be wrong.
• fusedims stack is written as

fusedims(a::AbstractArray, labels)  # labels=((1,),(2,3)) or tuplemortar(((1,),(2,3))))
    blockperms(labels)
    |-> fusedims(a::AbstractArray, bp::AbstractBlockPermutation)  # expect FusionTensors to define this one
            permutedims(a, Tuple(bp)); return fusedims(permuted, trivialperm(bp))
           |-> fusedims(a::AbstractArray, bp::TrivialBlockPermutation)
                new_axes = fuseaxes(axes(a), bp))
                |-> fusedims(a::AbstractArray, new_axes::AbstractUnitRanges...)
                     |-> fusedims(FusionStyle(a), new_axes)

-splidims is written using BlockedTuple{AbstractUnitRange} at its lowest level.

The associated BlockSparseArrays PR ITensor/BlockSparseArrays.jl#55 allows tests to pass for GradedUnitRanges axes.

TBD: get rid of _permutedims and _permutedims and replace them with dispatch on permutedims(a::AbstractArray, bp::BlockedPermutation

[codecov]

Codecov Report

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

Project coverage is 0.00%. Comparing base (9407722) to head (a8d078b).

Files with missing lines	Patch %	Lines
src/splitdims.jl	0.00%	16 Missing ⚠️
src/factorizations.jl	0.00%	6 Missing ⚠️
src/fusedims.jl	0.00%	5 Missing ⚠️
src/contract/blockedperms.jl	0.00%	3 Missing ⚠️
src/contract/contract_matricize/contract.jl	0.00%	3 Missing ⚠️
src/blockedpermutation.jl	0.00%	2 Missing ⚠️

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

HEAD has 7 uploads less than BASE

Flag	BASE (9407722)	HEAD (a8d078b)
docs	2	1
	6	0

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #33       +/-   ##
==========================================
- Coverage   93.95%   0.00%   -93.96%     
==========================================
  Files          15      13        -2     
  Lines         364     335       -29     
==========================================
- Hits          342       0      -342     
- Misses         22     335      +313     

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.

[mtfishman]

It definitely would be nice to get rid of the various _mul! functions handling combinations of matrices and vectors. I do remember considering the alternative of adding trailing singleton dimensions, however the reason I didn't go that route was that I thought it might be tricky to handle that for certain tensor types like block sparse or symmetric tensors (obviously not impossible, but I didn't want to make the ability to add/remove singleton dimensions a requirement of backend tensors, and it may not be easy to implement that in more complicated cases without copying data). How does this version handle that?

[ogauthe]

How does this version handle that?

The tests do not cover this case except for the native Array type. Currently, it is dealt with
⊗() = Base.OneTo(1) that fusedims calls. There are indeed issues with BlockSparseArrays. I will investigate and add more tests.

[ogauthe]

There is a real difficulty here. Creating a trivial axis ex nihilo is problematic as it may not fit with other axes used in the code or with the array. In the case of FusionTensors.jl, it was solved by designing TrivialSector with specific fusion rules.

• we may have a dedicated trivial_axis logic in TensorAlgebra, that may return Base.OneTo(1), blockedrange([1]) or gradedrange([TrivialSector() => 1]) (defined in package extensions) depending on input. We may use GradedUnitRanges.OneToOne for this role. In any case, this trivial axis should only use in intermediate steps and never appear in the ouput
• or we may avoid this burden by preserving filter(!isempty(permblocks. Then we are back to the current issue: a FusionTensor requires both domain and codomain and requires dedicated function. We may move the filter(!isempty logic deeper in the call stack, but FusionTensors.jl would need to use a custom path at some point.

[mtfishman]

• we may have a dedicated trivial_axis logic in TensorAlgebra, that may return Base.OneTo(1), blockedrange([1]) or gradedrange([TrivialSector() => 1]) (defined in package extensions) depending on input. We may use GradedUnitRanges.OneToOne for this role. In any case, this trivial axis should only use in intermediate steps and never appear in the ouput

That seems reasonable to me, though I don't have all of the context so I don't know all of the tradeoffs between the different choices. But I would call it singleton_axis (dimensions with length 1 are generally referred to as "singleton dimensions" in Julia, see for example https://docs.julialang.org/en/v1.13-dev/base/arrays/#Base.insertdims).