Inconsistent handling of empty tensors when using multifusion sectors

[leburgel]

There is some inconsistency in how tensor maps without any valid fusion channels are handled when using multifusion sectors. Depending on how the spaces of such a tensor are partitioned between the domain and codomain, construction and manipulation eithers returns an empty tensor with norm zero, or throws an ArgumentError.

For example:

using TensorKit, TensorKitSectors

C0, C1 = IsingBimodule(1, 1, 0), IsingBimodule(1, 1, 1)
D0, D1 = IsingBimodule(2, 2, 0), IsingBimodule(2, 2, 1)
M, Mop = IsingBimodule(1, 2, 0), IsingBimodule(2, 1, 0)

V = (
    Vect[IsingBimodule](C0 => 2, C1 => 1), Vect[IsingBimodule](Mop => 1)',
    Vect[IsingBimodule](D0 => 2, D1 => 2), Vect[IsingBimodule](M => 2)',
    Vect[IsingBimodule](C0 => 3, C1 => 2),
)
t = randn(Float64, V[1] ⊗ V[2] ← V[3] ⊗ V[4] ⊗ V[5])
# t is fine: norm(t) == 0.0, blocksectors(t) == IsingBimodule[]
@show norm(t)
@show blocksectors(t)

repartition(t, 0)   # ok, norm = 0.0
repartition(t, 1)   # ok, norm = 0.0
repartition(t, 2)   # ok, norm = 0.0
repartition(t, 3)   # ArgumentError: invalid fusion channel
repartition(t, 4)   # ArgumentError: invalid fusion channel
repartition(t, 5)   # ok, norm = 0.0

The issue seems to be that blocksectors calls corresponding to an empty set sometimes return and empty set and sometimes throw, depening on the specific arguments:

codom = Vect[IsingBimodule](C0 => 1) ⊗ Vect[IsingBimodule](C0 => 1)
blocksectors(codom ← Vect[IsingBimodule](D0 => 1))   # ArgumentError: invalid fusion channel
blocksectors(codom ← Vect[IsingBimodule](C1 => 1))   # IsingBimodule[]

I'm not sure whether it's better to always throw or never throw. @lkdvos, @borisdevos.

[borisdevos]

Just tagging my response in the PR where this issue originated, so we can discuss here: #513 (comment)

[lkdvos]

Looking at the last example of the blocksectors, that to me looks like the behavior we would expect, in the sense that the bottom line is an empty set of sectors because fusion is allowed but there is no overlap between domain and codomain, which is equivalent to U1Irrep(0) otimes U1Irrep(0) -> U1Irrep(1) behavior.
On the other hand, the first is not because this truly is an invalid fusion diagram, connecting regions with different colors together, which should error.

Given this, it might be the case of investigating a bit more in detail what is going on and being slightly more consistent with that, since I think the original tensor you construct is of the allowed but no overlap kind, so presumably all repartitions should indeed work and behave equivalently to something like t = rand(U1Space(0 => 1) * U1Space(0 => 1), U1Space(1 => 1)). (Which we might want to test works as well)

[Jutho]

But even for the initial tensor, the coloring is problematic for the double fusion tree. The domain tree has color 2 on the left and color 1 on the right, and the codomain tree has the opposite.

Are the spaces

V = (
    Vect[IsingBimodule](C0 => 2, C1 => 1), Vect[IsingBimodule](Mop => 1)',
    Vect[IsingBimodule](D0 => 2, D1 => 2), Vect[IsingBimodule](M => 2)',
    Vect[IsingBimodule](C0 => 3, C1 => 2),
)

those from (some of) the tests? I know that I at some point a few months ago redesigned the test spaces such that (V[1] ⊗ V[2] ⊗ V[3] ⊗ V[4] ⊗ V[5]) should be a valid product space, and then the test tensors should be constructed as either

t = randn(Float64, V[1] ⊗ V[2] ⊗ V[3] ⊗ V[4] ⊗ V[5])

or anything that is isomorphic to this, such as

t = randn(Float64, V[1] ⊗ V[2] ← (V[3] ⊗ V[4] ⊗ V[5])')

I think there is even a comment about this somewhere, if I recall correctly (and it survived further changes).

None of this is to say that we do not need fix certain issues in the code though.

[borisdevos]

Are the spaces... those from (some of) the tests?

The spaces are indeed from the tests, and are valid in the way you intended it to be. The way it was tested where it errored, namely V[1] ⊗ V[2] ← V[3] ⊗ V[4] ⊗ V[5] would indeed mismatch color-wise.

My point is more that this can error at the fusion tree/sector level by checking colors in a different way to checking N-symbols. We currently shield potential coloring errors at Nsymbol, which currently isn't a choice in the sense of, we may want to return 0 if the fusion is forbidden, but we can't right now. So if we can check forbidden fusions differently to Nsymbol where colorings matter (through ⊗, which returns an empty collection), then the expected behavior will appear.

At least for this specific example with blocksectors, replacing Nsymbol(a, b, c) calls with c ∈/∉ (a ⊗ b) depending on the context fixes it. It turns out that this kind of guard appears in some form everywhere else in the code as far as I can tell by just searching for Nsymbol, so it seems the way to go.