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?

[lkdvos]

The @planar macro does not actually know that T'` has trivial codomain, since it runs at parse time. You have to actually tell it that this is the shape:

@planar TT[ld lu ru rd] := T[ld u rd] * T'[(); lu u ru]

Note that it can sometimes be more convenient to write things in terms of conj, since that would keep the indices in the same place:

@planar TT[ld lu ru rd] := T[ld u rd] * conj(T[lu u ru])

[Chenqitrg]

Thanks for your answer.

A = Rep[ℤ₃](0=>2, 1=>3, 2=>1)
M = rand(ComplexF64, A←A)
@planar cont = M[a; a]

I am trying to perform the tensor trace for the planar tensor. The @tensor contract reports no error. But for @plannar, here is the error information:

ERROR: UndefVarError: `p` not defined
Stacktrace:
 [1] planartrace!(C::Tensor{…}, A::TensorMap{…}, ::Tuple{…}, ::Tuple{…}, α::VectorInterface.One, β::VectorInterface.Zero, backend::TensorOperations.DefaultBackend, allocator::TensorOperations.DefaultAllocator)
   @ TensorKit ~/.julia/packages/TensorKit/esxcj/src/planar/planaroperations.jl:71
 [2] planartrace!
   @ ~/.julia/packages/TensorKit/esxcj/src/planar/planaroperations.jl:43 [inlined]
 [3] planartrace!
   @ ~/.julia/packages/TensorKit/esxcj/src/planar/planaroperations.jl:39 [inlined]

Edit:

I have changed the p and q to (p1,p2) and (q1,q2). Then it is correct.

[lkdvos]

Thanks for reporting that, this is definitely a bug.

[Chenqitrg]

Does there exist a combiner tensor to combining several tensor legs into a single one?

[Jutho]

Not yet; but you can easily do this yourself with

Combiner = unitary(fuse(V1 * V2 * V3), V1 * V2 * V3)

This will create a specific TensorMap that can be used to fuse and unfuse several indices by contracting with it, or its dagger. There is some overhead, as you are contracting with essentially an identity matrix, so at some point there should be some functionality for doing this more efficiently, but often, this is not the cost-dominating step of an algorithm.

[Chenqitrg]

I am wondering for the VecG category, how can I use setindex? For example, suppose I have defined the VecD6 sector correctly. How should I set the index by sector?

e = D6(0,0)
r = D6(0,1)
r2 = D6(0,2)
s = D6(1,0)
sr = D6(1,1)
sr2 = D6(1,2)
A = Vect[D6](e=>1, r=>1, sr=>1, sr2=>2)
B = Vect[D6](e=>1, r=>1, sr=>2, sr2=>1)
C = Vect[D6](e=>1, r=>2, sr=>1, sr2=>1)
D = Vect[D6](e=>2, r=>1, sr=>1, sr2=>1)
T = rand(ComplexF64, A⊗B←C⊗D)

@show T[(e,e,e,e)] # This getindex function has no problem
T[(e,e,e,e)] = rand(1,1,1,2) #This reports an error.

Thanks Zhengyuan Yue for answering. Writing

T[(e,e,e,e)][:] = rand(1,1,1,2)

would be fine.

[Jutho]

Or T[(e,e,e,e)] .= rand(1,1,1,2)

[Chenqitrg]

Hi I noticed that for a fusion category without braiding structure, doing tsvd will always need to perform the permute function, which requires symmetric braiding. Is there a way to let tsvd using transpose to permute legs?

[Jutho]

If you just use tsvd(tensor) without additional permutation arguments, it shouldn't require any braiding structure. If you first want to transpose the tensor into a specific (p1, p2) structure, you will have to do this manually beforehand.