add experimental support for fermions and tensor contractions with super vector spaces

Author: Jutho

src/TensorKit.jl

@@ -13,7 +13,7 @@ export FusionStyle, UniqueFusion, MultipleFusion, MultiplicityFreeFusion,
         SimpleFusion, GenericFusion
 export BraidingStyle, SymmetricBraiding, Bosonic, Fermionic, Anyonic
 export Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep
-        #FermionParity, FermionNumber, FermionSpin # specific sectors
+export Fermion, FermionParity, FermionNumber, FermionSpin
 export FibonacciAnyon
 export IsingAnyon
 
@@ -45,6 +45,7 @@ export Trivial, ZNSpace, SU2Irrep, U1Irrep, CU1Irrep # Fermion
 # some unicode
 export ⊕, ⊗, ×, ⊠, ℂ, ℝ, ℤ, ←, →, ≾, ≿, ≅, ≺, ≻
 export ℤ₂, ℤ₃, ℤ₄, U₁, SU, SU₂, CU₁
+export fℤ₂, fU₁, fSU₂
 export ℤ₂Space, ℤ₃Space, ℤ₄Space, U₁Space, CU₁Space, SU₂Space
 
 # tensor maps

src/fusiontrees/manipulations.jl

@@ -875,7 +875,7 @@ function braid(f::FusionTree{I, N},
         for i = 1:N
             for j = 1:i-1
                 if p[j] > p[i]
-                    a, b = f.uncoupled[j], f.uncoupled[i]
+                    a, b = f.uncoupled[p[j]], f.uncoupled[p[i]]
                     coeff *= Rsymbol(a, b, first(a ⊗ b))
                 end
             end

src/sectors/fermions.jl

@@ -1 +1,92 @@
-# TODO
+struct Fermion{P,I<:Sector} <: Sector
+    sector::I
+    function Fermion{P,I}(sector::I) where {P, I<:Sector}
+        @assert BraidingStyle(I) isa Bosonic
+        return new{P,I}(sector)
+    end
+end
+Fermion{P}(sector::I) where {P, I<:Sector} = Fermion{P,I}(sector)
+Fermion{P,I}(sector) where {P, I<:Sector} = Fermion{P,I}(convert(I, sector))
+Base.convert(::Type{Fermion{P,I}}, a::Fermion{P,I}) where {P, I<:Sector} = a
+Base.convert(::Type{Fermion{P,I}}, a) where {P, I<:Sector} = Fermion{P,I}(convert(I, a))
+
+fermionparity(f::Fermion{P}) where P = P(f.sector)
+
+Base.IteratorSize(::Type{SectorValues{Fermion{P,I}}}) where {P, I<:Sector} =
+    Base.IteratorSize(SectorValues{I})
+Base.length(::SectorValues{Fermion{P,I}}) where {P, I<:Sector} = length(values(I))
+function Base.iterate(::SectorValues{Fermion{P, I}}) where {P, I<:Sector}
+    next = iterate(values(I))
+    @assert next !== nothing
+    value, state = next
+    return Fermion{P}(value), state
+end
+function Base.iterate(::SectorValues{Fermion{P, I}}, state) where {P, I<:Sector}
+    next = iterate(values(I), state)
+    if next === nothing
+        return nothing
+    else
+        value, state = next
+        return Fermion{P}(value), state
+    end
+end
+Base.getindex(::SectorValues{Fermion{P, I}}, i) where {P, I<:Sector} =
+    Fermion{P}(values(I)[i])
+findindex(::SectorValues{Fermion{P, I}}, f::Fermion{P, I}) where {P, I<:Sector} =
+    findindex(values(I), f.sector)
+
+Base.one(::Type{Fermion{P, I}}) where {P, I<:Sector} = Fermion{P}(one(I))
+Base.conj(f::Fermion{P}) where {P} = Fermion{P}(conj(f.sector))
+
+dim(f::Fermion) = dim(f.sector)
+
+FusionStyle(::Type{<:Fermion{<:Any,I}}) where {I<:Sector} = FusionStyle(I)
+BraidingStyle(::Type{<:Fermion}) = Fermionic()
+Base.isreal(::Type{Fermion{<:Any,I}}) where {I<:Sector} = isreal(I)
+
+⊗(a::F, b::F) where {F<:Fermion} = SectorSet{F}(a.sector ⊗ b.sector)
+
+Nsymbol(a::F, b::F, c::F) where {F<:Fermion} = Nsymbol(a.sector, b.sector, c.sector)
+
+Fsymbol(a::F, b::F, c::F, d::F, e::F, f::F) where {F<:Fermion} =
+    Fsymbol(a.sector, b.sector, c.sector, d.sector, e.sector, f.sector)
+
+function Rsymbol(a::F, b::F, c::F) where {F<:Fermion}
+    if fermionparity(a) && fermionparity(b)
+        return -Rsymbol(a.sector, b.sector, c.sector)
+    else
+        return +Rsymbol(a.sector, b.sector, c.sector)
+    end
+end
+
+twist(a::Fermion) = ifelse(fermionparity(a), -1, +1)*twist(a.sector)
+
+type_repr(::Type{Fermion{P,I}}) where {P, I<:Sector} = "Fermion{$P, " * type_repr(I) * "}"
+
+function Base.show(io::IO, a::Fermion{P, I}) where {P, I<:Sector}
+    if get(io, :typeinfo, nothing) !== Fermion{P, I}
+        print(io, type_repr(typeof(a)), "(")
+    end
+    print(IOContext(io, :typeinfo => I), a.sector)
+    if get(io, :typeinfo, nothing) !== Fermion{P, I}
+        print(io, ")")
+    end
+end
+
+Base.hash(f::Fermion, h::UInt) = hash(f.sector, h)
+Base.isless(a::F, b::F) where {F<:Fermion} = isless(a.sector, b.sector)
+
+_fermionparity(a::Z2Irrep) = isodd(a.n)
+_fermionnumber(a::U1Irrep) = isodd(convert(Int, a.charge))
+_fermionspin(a::SU2Irrep) = isodd(twice(a.j))
+
+const FermionParity = Fermion{_fermionparity, Z2Irrep}
+const FermionNumber = Fermion{_fermionnumber, U1Irrep}
+const FermionSpin = Fermion{_fermionspin, SU2Irrep}
+const fℤ₂ = FermionParity
+const fU₁ = FermionNumber
+const fSU₂ = FermionSpin
+
+type_repr(::Type{FermionParity}) = "FermionParity"
+type_repr(::Type{FermionNumber}) = "FermionNumber"
+type_repr(::Type{FermionSpin}) = "FermionSpin"

src/sectors/irreps.jl

@@ -57,6 +57,7 @@ Nsymbol(a::I, b::I, c::I) where {I<:AbelianIrrep} = c == first(a ⊗ b)
 Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I<:AbelianIrrep} =
     Int(Nsymbol(a, b, e)*Nsymbol(e, c, d)*Nsymbol(b, c, f)*Nsymbol(a, f, d))
 frobeniusschur(a::AbelianIrrep) = 1
+Asymbol(a::I, b::I, c::I) where {I<:AbelianIrrep} = Int(Nsymbol(a, b, c))
 Bsymbol(a::I, b::I, c::I) where {I<:AbelianIrrep} = Int(Nsymbol(a, b, c))
 Rsymbol(a::I, b::I, c::I) where {I<:AbelianIrrep} = Int(Nsymbol(a, b, c))
 

src/sectors/sectors.jl

@@ -98,7 +98,11 @@ which case it is complex).
 """
 function Base.isreal(I::Type{<:Sector})
     u = one(I)
-    return (eltype(Fsymbol(u, u, u, u, u, u))<:Real) && (eltype(Rsymbol(u, u, u))<:Real)
+    if BraidingStyle(I) isa HasBraiding
+        return (eltype(Fsymbol(u, u, u, u, u, u))<:Real) && (eltype(Rsymbol(u, u, u))<:Real)
+    else
+        return (eltype(Fsymbol(u, u, u, u, u, u))<:Real)
+    end
 end
 Base.isreal(::Type{Trivial}) = true
 
@@ -379,6 +383,6 @@ end
 # possible sectors
 include("groups.jl")
 include("irreps.jl") # irreps of symmetry groups, with bosonic braiding
-# include("fermions.jl") # irreps with defined fermionparity and fermionic braiding
+include("fermions.jl") # irreps with defined fermionparity and fermionic braiding
 include("anyons.jl") # non-group sectors
 include("product.jl") # direct product of different sectors

src/tensors/indexmanipulations.jl

@@ -11,13 +11,12 @@ codomain or range of the map, and indices in `p2` indicating the domain.
 To permute into an existing `tdst`, see [`add!`](@ref)
 """
 function permute(t::TensorMap{S},
-                    p1::IndexTuple, p2::IndexTuple=();
-                    copy::Bool = false) where {S}
-    cod = ProductSpace{S}(map(n->space(t, n), p1))
-    dom = ProductSpace{S}(map(n->dual(space(t, n)), p2))
-
+                    p1::IndexTuple{N₁},  p2::IndexTuple{N₂}=();
+                    copy::Bool = false) where {S, N₁, N₂}
+    cod = ProductSpace{S, N₁}(map(n->space(t, n), p1))
+    dom = ProductSpace{S, N₂}(map(n->dual(space(t, n)), p2))
+    # share data if possible
     if !copy
-        # share data if possible
         if p1 === codomainind(t) && p2 === domainind(t)
             return t
         elseif has_shared_permute(t, p1, p2)

src/tensors/tensor.jl

@@ -79,7 +79,7 @@ function TensorMap(data::DenseArray, codom::ProductSpace{S,N₁}, dom::ProductSp
         size(data) == (dims(codom)..., dims(dom)...))
         throw(DimensionMismatch())
     end
-    if sectortype(S) === Trivial # For now, we can only accept array data for Trivial sectortype
+    if sectortype(S) === Trivial
         data2 = reshape(data, (d1, d2))
         A = typeof(data2)
         return TensorMap{S, N₁, N₂, Trivial, A, Nothing, Nothing}(data2, codom, dom)

src/tensors/tensoroperations.jl

@@ -1,17 +1,13 @@
 function cached_permute(sym::Symbol, t::TensorMap{S},
-                            p1::IndexTuple{N₁},  p2::IndexTuple{N₂}=()) where {S, N₁, N₂}
+                            p1::IndexTuple{N₁},  p2::IndexTuple{N₂}=();
+                            copy::Bool = false) where {S, N₁, N₂}
     cod = ProductSpace{S, N₁}(map(n->space(t, n), p1))
     dom = ProductSpace{S, N₂}(map(n->dual(space(t, n)), p2))
-
     # share data if possible
-    if p1 === codomainind(t) && p2 === domainind(t)
-        return t
-    elseif isa(t, TensorMap) && sectortype(S) === Trivial
-        stridet = i->stride(t[], i)
-        sizet = i->size(t[], i)
-        canfuse1, d1, s1 = TensorOperations._canfuse(sizet.(p1), stridet.(p1))
-        canfuse2, d2, s2 = TensorOperations._canfuse(sizet.(p2), stridet.(p2))
-        if canfuse1 && canfuse2 && s1 == 1 && (d2 == 1 || s2 == d1)
+    if !copy
+        if p1 === codomainind(t) && p2 === domainind(t)
+            return t
+        elseif has_shared_permute(t, p1, p2)
             return TensorMap(reshape(t.data, dim(cod), dim(dom)), cod, dom)
         end
     end
@@ -23,11 +19,11 @@ function cached_permute(sym::Symbol, t::TensorMap{S},
 end
 
 function cached_permute(sym::Symbol, t::AdjointTensorMap{S},
-                            p1::IndexTuple{N₁},  p2::IndexTuple{N₂}=()) where {S, N₁, N₂}
-
+                            p1::IndexTuple,  p2::IndexTuple=();
+                            copy::Bool = false) where {S, N₁, N₂}
     p1′ = adjointtensorindices(t, p2)
     p2′ = adjointtensorindices(t, p1)
-    adjoint(cached_permute(sym, adjoint(t), p1′, p2′))
+    adjoint(cached_permute(sym, adjoint(t), p1′, p2′; copy = copy))
 end
 
 scalar(t::AbstractTensorMap{S}) where {S<:IndexSpace} =
@@ -254,21 +250,21 @@ function contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S},
     memcost4 += dB*(!hsp(B, oindB, cindB′′)) +
                 dA*(!hsp(A, cindA′′, oindA))
 
-    if min(memcost1, memcost2) <= min(memcost3, memcost4)
+    # if min(memcost1, memcost2) <= min(memcost3, memcost4)
         if memcost1 <= memcost2
             return _contract!(α, A, B, β, C, oindA, cindA′, oindB, cindB′, p1, p2, syms)
         else
             return _contract!(α, A, B, β, C, oindA, cindA′′, oindB, cindB′′, p1, p2, syms)
         end
-    else
-        p1′ = map(n->ifelse(n>N₁, n-N₁, n+N₂), p1)
-        p2′ = map(n->ifelse(n>N₁, n-N₁, n+N₂), p2)
-        if memcost3 <= memcost4
-            return _contract!(α, B, A, β, C, oindB, cindB′, oindA, cindA′, p1′, p2′, syms)
-        else
-            return _contract!(α, B, A, β, C, oindB, cindB′′, oindA, cindA′′, p1′, p2′, syms)
-        end
-    end
+    # else
+    #     p1′ = map(n->ifelse(n>N₁, n-N₁, n+N₂), p1)
+    #     p2′ = map(n->ifelse(n>N₁, n-N₁, n+N₂), p2)
+    #     if memcost3 <= memcost4
+    #         return _contract!(α, B, A, β, C, oindB, cindB′, oindA, cindA′, p1′, p2′, syms)
+    #     else
+    #         return _contract!(α, B, A, β, C, oindB, cindB′′, oindA, cindA′′, p1′, p2′, syms)
+    #     end
+    # end
 end
 
 function _contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S},
@@ -278,13 +274,31 @@ function _contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S},
                     p1::IndexTuple, p2::IndexTuple,
                     syms::Union{Nothing, NTuple{3, Symbol}} = nothing) where {S, N₁, N₂}
 
+    if !(BraidingStyle(sectortype(S)) isa SymmetricBraiding)
+        throw(SectorMismatch("only tensors with symmetric braiding rules can be contracted; try `@planar` instead"))
+    end
+    copyA = false
+    if BraidingStyle(sectortype(S)) isa Fermionic
+        for i in cindA
+            if !isdual(space(A, i))
+                copyA = true
+            end
+        end
+    end
     if syms === nothing
-        A′ = permute(A, oindA, cindA)
+        A′ = permute(A, oindA, cindA; copy = copyA)
         B′ = permute(B, cindB, oindB)
     else
-        A′ = cached_permute(syms[1], A, oindA, cindA)
+        A′ = cached_permute(syms[1], A, oindA, cindA; copy = copyA)
         B′ = cached_permute(syms[2], B, cindB, oindB)
     end
+    if BraidingStyle(sectortype(S)) isa Fermionic
+        for i in domainind(A′)
+            if !isdual(space(A′, i))
+                A′ = twist!(A′, i)
+            end
+        end
+    end
     ipC = TupleTools.invperm((p1..., p2...))
     oindAinC = TupleTools.getindices(ipC, ntuple(n->n, N₁))
     oindBinC = TupleTools.getindices(ipC, ntuple(n->n+N₁, N₂))

test/runtests.jl

@@ -19,6 +19,7 @@ const TK = TensorKit
 Random.seed!(1234)
 
 smallset(::Type{I}) where {I<:Sector} = take(values(I), 5)
+smallset(::Type{FermionNumber}) = FermionNumber.((0, +1, -1, +2, -2))
 function smallset(::Type{ProductSector{Tuple{I1,I2}}}) where {I1,I2}
     iter = product(smallset(I1), smallset(I2))
     s = collect(i ⊠ j for (i,j) in iter if dim(i)*dim(j) <= 6)
@@ -51,9 +52,10 @@ function hasfusiontensor(I::Type{<:Sector})
 end
 
 sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewSU2Irrep, SU3Irrep,
-              FibonacciAnyon, IsingAnyon, Z3Irrep ⊠ Z4Irrep, U1Irrep ⊠ SU2Irrep, SU2Irrep ⊠
-              SU2Irrep, NewSU2Irrep ⊠ NewSU2Irrep, NewSU2Irrep ⊠ SU2Irrep, SU2Irrep ⊠
-              NewSU2Irrep, Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon)
+              FibonacciAnyon, IsingAnyon, FermionParity, FermionNumber, FermionSpin,
+              FermionParity ⊠ FermionParity, Z3Irrep ⊠ Z4Irrep, FermionNumber ⊠ SU2Irrep,
+              FermionSpin ⊠ SU2Irrep, NewSU2Irrep ⊠ NewSU2Irrep, NewSU2Irrep ⊠ SU2Irrep,
+              FermionSpin ⊠ NewSU2Irrep, Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon)
 
 Ti = time()
 include("sectors.jl")

test/tensors.jl

@@ -8,6 +8,11 @@ Vℤ₂ = (ℂ[Z2Irrep](0=>1, 1=>1),
         ℂ[Z2Irrep](0=>3, 1=>2)',
         ℂ[Z2Irrep](0=>2, 1=>3),
         ℂ[Z2Irrep](0=>2, 1=>5))
+Vfℤ₂ = (ℂ[FermionParity](0=>1, 1=>1),
+        ℂ[FermionParity](0=>1, 1=>2)',
+        ℂ[FermionParity](0=>3, 1=>2)',
+        ℂ[FermionParity](0=>2, 1=>3),
+        ℂ[FermionParity](0=>2, 1=>5))
 Vℤ₃ = (ℂ[Z3Irrep](0=>1, 1=>2, 2=>2),
         ℂ[Z3Irrep](0=>3, 1=>1, 2=>1),
         ℂ[Z3Irrep](0=>2, 1=>2, 2=>1)',
@@ -18,6 +23,11 @@ VU₁ = (ℂ[U1Irrep](0=>1, 1=>2, -1=>2),
         ℂ[U1Irrep](0=>2, 1=>2, -1=>1)',
         ℂ[U1Irrep](0=>1, 1=>2, -1=>3),
         ℂ[U1Irrep](0=>1, 1=>3, -1=>3)')
+VfU₁ = (ℂ[FermionNumber](0=>1, 1=>2, -1=>2),
+        ℂ[FermionNumber](0=>3, 1=>1, -1=>1),
+        ℂ[FermionNumber](0=>2, 1=>2, -1=>1)',
+        ℂ[FermionNumber](0=>1, 1=>2, -1=>3),
+        ℂ[FermionNumber](0=>1, 1=>3, -1=>3)')
 VCU₁ = (ℂ[CU1Irrep]((0,0)=>1, (0,1)=>2, 1=>1),
         ℂ[CU1Irrep]((0,0)=>3, (0,1)=>0, 1=>1),
         ℂ[CU1Irrep]((0,0)=>1, (0,1)=>0, 1=>2)',
@@ -28,13 +38,18 @@ VSU₂ = (ℂ[SU2Irrep](0=>3, 1//2=>1),
         ℂ[SU2Irrep](1//2=>1, 1=>1)',
         ℂ[SU2Irrep](0=>2, 1//2=>2),
         ℂ[SU2Irrep](0=>1, 1//2=>1, 3//2=>1)')
+VfSU₂ = (ℂ[FermionSpin](0=>3, 1//2=>1),
+            ℂ[FermionSpin](0=>2, 1=>1),
+            ℂ[FermionSpin](1//2=>1, 1=>1)',
+            ℂ[FermionSpin](0=>2, 1//2=>2),
+            ℂ[FermionSpin](0=>1, 1//2=>1, 3//2=>1)')
 VSU₃ = (ℂ[SU3Irrep]((0,0,0)=>3, (1,0,0)=>1),
                ℂ[SU3Irrep]((0,0,0)=>3, (2,0,0)=>1)',
                ℂ[SU3Irrep]((1,1,0)=>1, (2,1,0)=>1),
                ℂ[SU3Irrep]((1,0,0)=>1, (2,0,0)=>1),
                ℂ[SU3Irrep]((0,0,0)=>1, (1,0,0)=>1, (1,1,0)=>1)')
 
-for V in (Vtr, Vℤ₂, Vℤ₃, VU₁, VCU₁, VSU₂)#, VSU₃)
+for V in (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂)#, VSU₃)
     V1, V2, V3, V4, V5 = V
     @assert V3 * V4 * V2 ≿ V1' * V5' # necessary for leftorth tests
     @assert V3 * V4 ≾ V1' * V2' * V5' # necessary for rightorth tests