Can the category VecG be supported?

[Chenqitrg]

Hi thanks for your guys' good package.

I am writing a code, whose tensors are morphisms in VecG, the category of G-graded vector spaces. It is not a braided category so there is only the fusion structure.

The following is the code to define the VecD8 category structure following your documentation.

using TensorKit

struct D8 <: Sector
    s::Int
    r::Int
end

Base.one(::Type{D8}) = D8(0,0)
Base.conj(g::D8) = D8(mod(-g.s,2), mod((-1)^(g.s + 1) * g.r,4))
dim(g::D8) = 1
Base.isreal(::Type{D8}) = true
TensorKit.FusionStyle(::Type{D8}) = UniqueFusion()
TensorKit.Nsymbol(a::D8, b::D8, c::D8) = (mod(a.s + b.s, 2), mod((-1)^b.s * a.r + b.r, 4)) == (c.s, c.r) ? 1 : 0
⊗(a::D8, b::D8) = D8(mod(a.s + b.s, 2), mod((-1)^b.s * a.r + b.r, 4))
Base.isless(a::D8, b::D8) = a.s < b.s || (a.s == b.s && a.r < b.r)
function Base.iterate(::TensorKit.SectorValues{D8},s=0,r=0)
    if s==0 && r < 4
        return (D8(s,r), (s, r+1))
    elseif s==0 && r==4
        return (D8(0,4), (1,0))
    elseif s==1 && r < 4
        return (D8(s,r), (s, r+1))
    else
        return nothing
    end
end
TensorKit.Fsymbol(a::D8, b::D8, c::D8, d::D8, e::D8, f::D8) = 1
Base.length(::TensorKit.SectorValues{D8}) = 8
Base.getindex(::TensorKit.SectorValues{D8}, i::Int) = D8((i-1)÷4, (i-1)%4)
TensorKit.findindex(::TensorKit.SectorValues{D8}, c::D8) = c.r + 4 * c.s + 1
V1 = Vect[D8](D8(0,0) => 1, D8(1,2) => 2)
V2 = Vect[D8](D8(0,0) => 1, D8(1,2) => 2)
V3 = Vect[D8](D8(0,0) => 1, D8(1,2) => 2)
V4 = Vect[D8](D8(0,0) => 1, D8(1,2) => 2)

function ⊗(V1::GradedSpace{D8, NTuple{8, Int64}}, V2::GradedSpace{D8, NTuple{8, Int64}})
    G = TensorKit.SectorValues{D8}()
    multiplicity = zeros(Int, 8)
    for i in 1:8, j in 1:8
        delta = V1.dims[i] * V2.dims[j]
        g = G[i] ⊗ G[j]
        multiplicity[TensorKit.findindex(G, g)] += delta
    end
    return GradedSpace{D8, NTuple{8, Int64}}(Tuple(multiplicity), false)
end

@show D8(1,2) ⊗ D8(0,0)
@show V1 ⊗ V4
rand(V1 ⊗ V4 ← V2 ⊗ V3)

It reports error:

ERROR: ArgumentError: Union{} does not have elements
Stacktrace:
 [1] eltype(::Type{Union{}})
   @ Base ./abstractarray.jl:240
 [2] sectorscalartype(::Type{D8})
   @ TensorKitSectors ~/.julia/packages/TensorKitSectors/cqTdc/src/sectors.jl:112
 [3] TensorMap{Float64, GradedSpace{…}, 1, 1, Vector{…}}(data::Vector{Float64}, space::TensorMapSpace{GradedSpace{…}, 1, 1})
   @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/tensors/tensor.jl:34
 [4] TensorMap{Float64, GradedSpace{…}, 1, 1, Vector{…}}(::UndefInitializer, space::TensorMapSpace{GradedSpace{…}, 1, 1})
   @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/tensors/tensor.jl:24
 [5] (TensorMap{…} where {…})(::UndefInitializer, V::TensorMapSpace{…})
   @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/tensors/tensor.jl:90
 [6] rand(rng::Random.TaskLocalRNG, ::Type{Float64}, V::TensorMapSpace{GradedSpace{D8, NTuple{8, Int64}}, 1, 1})
   @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/tensors/tensor.jl:259
 [7] rand(::Type{Float64}, V::TensorMapSpace{GradedSpace{D8, NTuple{8, Int64}}, 1, 1})
   @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/tensors/tensor.jl:252
 [8] rand(V::TensorMapSpace{GradedSpace{D8, NTuple{8, Int64}}, 1, 1})
   @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/tensors/tensor.jl:245

I am wondering whether it can be supported with only a small modification? I am also writing a code, but I am only trying to tackle VecG case.

[Chenqitrg]

Oh writing

TensorKit.BraidingStyle(::Type{D8}) = NoBraiding()

will be good.

[lkdvos]

Thanks for reporting this, this actually looks like a small bug in TensorKitSectors, where our automatic fallback is failing. I'll look into fixing that asap.
You can fix this very easily in the meantime by defining TensorKitSectors.sectorscalartype(::Type{D8}) = Int, which should reflect the scalartype of the Fsymbol.

[edit]: since Fsymbol returns Int in your case it should actually be Int.

[Chenqitrg]

A related question: the current tensor product of Graded spaces is written by my own. Directly using the tensor product without definition will report an error:

MethodError: no method matching ⊗(::GradedSpace{D8, NTuple{8, Int64}}, ::GradedSpace{D8, NTuple{8, Int64}})

And my way of realizing the tensor product may not in a proper way, since the correct output should be a ProductSpace.

[lkdvos]

I might be mistaken, but could it be that you are defining ⊗(::D8, ::D8) instead of TensorKitSectors.⊗(::D8, ::D8)?

[Chenqitrg]

Thanks for answering.
This is indeed not following the standard one in Sectors . If I write

Base.:⊗(a::D8, b::D8) = D8(mod(a.s + b.s, 2), mod((-1)^b.s * a.r + b.r, 4))

as in documentation, then there is a error

UndefVarError: `⊗` not defined

Edit: OK writting

TensorKitSectors.:⊗(a::D8, b::D8) = D8(mod(a.s + b.s, 2), mod((-1)^b.s * a.r + b.r, 4))

would be good.

Edit:
Now A = rand(V1 ← V1 ⊗ V2) can pass. However, A = rand(V1 ⊗ V2 ← V1 ⊗ V2) reports an error

ERROR: MethodError: no method matching iterate(::D8)

Closest candidates are:
  iterate(::Pkg.Resolve.NodePerm, Any...)
   @ Pkg ~/.julia/juliaup/julia-1.10.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/Pkg/src/Resolve/maxsum.jl:240
  iterate(::TensorKitSectors.CU1ProdIterator)
   @ TensorKitSectors ~/.julia/packages/TensorKitSectors/cqTdc/src/irreps.jl:311
  iterate(::TensorKitSectors.CU1ProdIterator, ::Int64)
   @ TensorKitSectors ~/.julia/packages/TensorKitSectors/cqTdc/src/irreps.jl:311
  ...

Stacktrace:
  [1] blocksectors(P::ProductSpace{GradedSpace{D8, NTuple{8, Int64}}, 2})
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/spaces/productspace.jl:158
  [2] blocksectors
    @ ~/.julia/packages/TensorKit/esxcj/src/spaces/homspace.jl:93 [inlined]
  [3] fusionblockstructure(W::TensorMapSpace{GradedSpace{D8, NTuple{8, Int64}}, 2, 2}, ::TensorKit.NoCache)
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/spaces/homspace.jl:273
  [4] #150
    @ ~/.julia/packages/TensorKit/esxcj/src/spaces/homspace.jl:353 [inlined]
  [5] get!(default::TensorKit.var"#150#151"{…}, lru::LRUCache.LRU{…}, key::TensorMapSpace{…})
    @ LRUCache ~/.julia/packages/LRUCache/ctUcD/src/LRUCache.jl:169
  [6] fusionblockstructure(W::TensorMapSpace{GradedSpace{D8, NTuple{8, Int64}}, 2, 2}, ::TensorKit.GlobalLRUCache)
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/spaces/homspace.jl:352
  [7] fusionblockstructure
    @ ~/.julia/packages/TensorKit/esxcj/src/spaces/homspace.jl:251 [inlined]
  [8] TensorMap
    @ ~/.julia/packages/TensorKit/esxcj/src/tensors/tensor.jl:19 [inlined]
  [9] TensorMap
    @ ~/.julia/packages/TensorKit/esxcj/src/tensors/tensor.jl:90 [inlined]
 [10] rand
    @ ~/.julia/packages/TensorKit/esxcj/src/tensors/tensor.jl:259 [inlined]
 [11] rand
    @ ~/.julia/packages/TensorKit/esxcj/src/tensors/tensor.jl:252 [inlined]
 [12] rand(V::TensorMapSpace{GradedSpace{D8, NTuple{8, Int64}}, 2, 2})
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/tensors/tensor.jl:245

[Chenqitrg]

I isolated the bug function and deleted the irrelevant lines, it seems that ⊗(s...) can not be used as an iterator because it is uniquefusion. Changing the definition to be a tuple would be ok:

TensorKitSectors.:⊗(a::D8, b::D8) = (D8(mod(a.s + b.s, 2), mod((-1)^b.s * a.r + b.r, 4)),)

[Chenqitrg]

Hi, please forgive me for my many naive questions.
Since the category VecD8 does not have braiding, its permutation may be problematic. It seems that in the permute function, there is an assertion of symmetric braiding.

However, one can still perform the plannar permutation for VecD8 tensors, that is, it is allowed to write permute(T, (4,1), (2,3)). Can I do that?

[Jutho]

Yes, there is the transpose function for that: transpose only uses the duality structure without using braiding and supports "cyclic permutations" (graphically cyclic, due the way we label the indices it doesn't look very cyclic when you look at how the permutation is specified). It will error if the "permutation" is not planar/cyclic.

[Chenqitrg]

In function tsvd, can I simultaneously give the truncation dimension and truncation error, such that the truncation is by the minimum of these two?

[lkdvos]

yes, these truncation schemes support the following syntax: truncdim(maxdim) & truncbelow(minval).
(see also the documentation page, where this is mentioned)

[Chenqitrg]

For the convenience of future generalization, I am trying to use the @planar to contract tensors. The planar tensor I am trying to get is of the form

--lu--T†-ru--
      |
      u
      |      
--ld--T--rd--

A = Rep[ℤ₃](0=>2, 1=>3, 2=>1)
B = Rep[ℤ₃](0=>4, 1=>5, 2=>2)
C = Rep[ℤ₃](0=>1, 1=>2, 2=>3)
T = rand(ComplexF64, A⊗B⊗C←one(Rep[ℤ₃]))
@planar TT[ld lu; rd ru] := T[ld u rd] * T'[lu u ru]

It reports error

ERROR: LoadError: ArgumentError: not a planar diagram expression: TT[ld lu; rd ru] := T[ld u rd] * T'[lu u ru]
Stacktrace:
  [1] not_planar_err(ex::Any)
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/planar/preprocessors.jl:2
  [2] _check_planarity(ex::Any)
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/planar/preprocessors.jl:398
  [3] (::TensorKit.var"#353#355")(a::Any)
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/planar/preprocessors.jl:402
  [4] foreach(f::TensorKit.var"#353#355", itr::Vector{Any})
    @ Base ./abstractarray.jl:3097
  [5] _check_planarity(ex::Any)
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/planar/preprocessors.jl:401
  [6] (::TensorKit.var"#353#355")(a::Any)
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/planar/preprocessors.jl:402
  [7] foreach(f::TensorKit.var"#353#355", itr::Vector{Any})
    @ Base ./abstractarray.jl:3097
  [8] _check_planarity(ex::Any)
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/planar/preprocessors.jl:401
  [9] (::TensorKit.var"#377#388")(ex::Any)
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/planar/macros.jl:85
 [10] (::TensorOperations.TensorParser)(ex::Expr)
    @ TensorOperations ~/.julia/packages/TensorOperations/Pw3eO/src/indexnotation/parser.jl:28
 [11] var"@planar"(__source__::LineNumberNode, __module__::Module, args::Vararg{Expr})
    @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/planar/macros.jl:9
 [12] eval
    @ ./boot.jl:385 [inlined]
 [13] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:2076

But it seems to be a plannar diagram.

[lkdvos]

Your drawing does not correspond to the network. You have to draw a tensor t[l1 l2 l3; r1 r2 r3] (also use ; for the input tensors) like this:

     ┌───────┐      
     │       │      
l1 ──┼       ├─── r1
     │       │      
     │       │      
l2 ──┼       ├─── r2
     │       │      
     │       │      
l3 ──┼       ├─── r3
     │       │      
     └───────┘      

[Chenqitrg]

Thanks for your answer. I am still confused, in what case will the duality be automatically applied? Or do I need to apply the transpose first? The following are ok, except the third one.

A = Rep[ℤ₃](0=>2, 1=>3, 2=>1)
B = Rep[ℤ₃](0=>4, 1=>5, 2=>2)
C = Rep[ℤ₃](0=>1, 1=>2, 2=>3)
T = rand(ComplexF64, B⊗A⊗C←one(Rep[ℤ₃]))
F = rand(ComplexF64, B⊗ A'←B)
@planar K[a b; c] := T[a b c]
@planar G[1 3; 5 2] := F[1 a; 5] * T[2 a 3]
@planar TT[ld lu ru rd] := T[ld u rd] * T'[lu u ru]

Note that the T has trivial domain and T' has trivial codomain. So it should be a planar graph. Perhaps when domain or codomain are empty, the automatic check will treat it as the empty domain rather than codomain?