From c61b531a3f9c68b4184ff3f619509f593c21f7a5 Mon Sep 17 00:00:00 2001 From: maarten Date: Thu, 12 Aug 2021 14:43:52 +0200 Subject: [PATCH 01/15] can hang --- src/fusiontrees/manipulations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fusiontrees/manipulations.jl b/src/fusiontrees/manipulations.jl index 9f67c1306..3a17e2405 100644 --- a/src/fusiontrees/manipulations.jl +++ b/src/fusiontrees/manipulations.jl @@ -541,7 +541,7 @@ function _transpose((f1, f2, p1, p2)::TransposeKey{I,N₁,N₂}) where {I<:Secto newtrees = newtrees′ i1 -= 1 end - while Nhalf < i1 + while Nhalf < i1 && Nhalf > 0 local newtrees′ for ((f1a, f2a), coeffa) in newtrees for ((f1b, f2b), coeffb) in cycleclockwise(f1a, f2a) From 71474964ba04c82f6b0151608e8a3c077f8d5219 Mon Sep 17 00:00:00 2001 From: maarten Date: Thu, 12 Aug 2021 14:44:08 +0200 Subject: [PATCH 02/15] guessed --- src/tensors/braidingtensor.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index ac75b14a5..87e1442a5 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -65,7 +65,6 @@ function fusiontrees(b::BraidingTensor) end end dim2 = offset2 - push!(data, c=>f((dim1, dim2))) push!(rowr, c=>rowrc) push!(colr, c=>colrc) end @@ -78,8 +77,8 @@ end V1, V2 = domain(b) @boundscheck begin c == f2.coupled || throw(SectorMismatch()) - ((f1.uncoupled[1] ∈ V2) && (f2.uncoupled[1] ∈ V1)) || throw(SectorMismatch()) - ((f1.uncoupled[2] ∈ V1) && (f2.uncoupled[2] ∈ V2)) || throw(SectorMismatch()) + ((f1.uncoupled[1] ∈ sectors(V2)) && (f2.uncoupled[1] ∈ sectors(V1))) || throw(SectorMismatch()) + ((f1.uncoupled[2] ∈ sectors(V1)) && (f2.uncoupled[2] ∈ sectors(V2))) || throw(SectorMismatch()) end @inbounds begin d = (dims(V2 ⊗ V1, f1.uncoupled)..., dims(V1 ⊗ V2, f2.uncoupled)...) @@ -89,10 +88,11 @@ end data = fill!(storagetype(b)(undef, (n1, n2)), zero(eltype(b))) if f1.uncoupled == (a2, a1) braiddict = artin_braid(f2, 1; inv = b.adjoint) - r = get(dict, f1, zero(valtype(braiddict))) + r = get(braiddict, f1, zero(valtype(braiddict))) data[1:(n1+1):end] .= r # set diagonal end - return permutedims(sreshape(StridedView(data), d), (1,2,4,3)) + + return sreshape(StridedView(data), d) end end From 7ea1b35450177c8cae08ec8b74ef16177d120d83 Mon Sep 17 00:00:00 2001 From: maarten Date: Thu, 12 Aug 2021 14:44:18 +0200 Subject: [PATCH 03/15] istensor comes from tensoroperations --- src/tensors/planar.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tensors/planar.jl b/src/tensors/planar.jl index 30e865558..df11a4626 100644 --- a/src/tensors/planar.jl +++ b/src/tensors/planar.jl @@ -409,7 +409,7 @@ function _construct_braidingtensors(ex::Expr) lhs, rhs = TO.getlhs(ex), TO.getrhs(ex) if TO.istensorexpr(rhs) list = TO.gettensors(rhs) - if TO.isassignment(ex) && istensor(lhs) + if TO.isassignment(ex) && TO.istensor(lhs) obj, l, r = TO.decomposetensor(lhs) lhs_as_rhs = Expr(:typed_vcat, Expr(TO.prime, obj), Expr(:tuple, r...), Expr(:tuple, l...)) push!(list, lhs_as_rhs) From 7c2ee7df70bfbc3314d55ad4c1bd930bf45983e8 Mon Sep 17 00:00:00 2001 From: maarten Date: Thu, 12 Aug 2021 16:45:50 +0200 Subject: [PATCH 04/15] Trivial braidingtensors are difficult to get right --- src/tensors/braidingtensor.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index 87e1442a5..adc903fdd 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -71,6 +71,12 @@ function fusiontrees(b::BraidingTensor) return TensorKeyIterator(rowr, colr) end +function Base.getindex(b::BraidingTensor{S}) where S + sectortype(S) == Trivial || throw(SectorMismatch()) + (V1,V2) = domain(b); + return reshape(storagetype(b)(LinearAlgebra.I, dim(V1)*dim(V2), dim(V1)*dim(V2)),(dim(V2),dim(V1),dim(V1),dim(V2))); +end + @inline function Base.getindex(b::BraidingTensor, f1::FusionTree{I,2}, f2::FusionTree{I,2}) where {I<:Sector} I == sectortype(b) || throw(SectorMismatch()) c = f1.coupled @@ -101,6 +107,11 @@ function Base.copy!(t::TensorMap, b::BraidingTensor) space(t) == space(b) || throw(SectorMismatch()) fill!(t, zero(eltype(t))) for (f1, f2) in fusiontrees(t) + if f1 == nothing || f2 == nothing + copy!(t.data,one(t.data)) + return t + end + a1, a2 = f2.uncoupled c = f2.coupled f1.uncoupled == (a2, a1) || continue From 53b7a2b3e19a6e564bd8414a47c6b8d6e944d306 Mon Sep 17 00:00:00 2001 From: maarten Date: Thu, 12 Aug 2021 19:01:34 +0200 Subject: [PATCH 05/15] special case artin braiding for trivial legs --- src/fusiontrees/manipulations.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/fusiontrees/manipulations.jl b/src/fusiontrees/manipulations.jl index 3a17e2405..a73e3287f 100644 --- a/src/fusiontrees/manipulations.jl +++ b/src/fusiontrees/manipulations.jl @@ -766,7 +766,32 @@ function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I<:Sector inner = f.innerlines vertices = f.vertices u = one(I) + + if u in (uncoupled[i],uncoupled[i+1]) # the braid simplifies drastically, we are braiding a trivial sector + a, b = uncoupled[i], uncoupled[i+1] + uncoupled′ = TupleTools.setindex(uncoupled,b,i); + uncoupled′ = TupleTools.setindex(uncoupled′,a,i+1); + vertices′ = vertices; + + if i > 1 #we also need to alter innerlines and vertices + incharges = (uncoupled[1],inner...,coupled′); + + vertices′ = TupleTools.setindex(vertices′,vertices[i],i-1) + vertices′ = TupleTools.setindex(vertices′,vertices[i-1],i) + + if a == u + inner = TupleTools.setindex(inner,incharges[i+1],i-1); + else + inner = TupleTools.setindex(inner,incharges[i-1],i-1); + end + end + + f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices′) + return fusiontreedict(I)(f′ => true) + end + oneT = one(eltype(Rsymbol(u,u,u))) * one(eltype(Fsymbol(u,u,u,u,u,u))) + if i == 1 a, b = uncoupled[1], uncoupled[2] c = N > 2 ? inner[1] : coupled′ From cae2d8ff4d5b0762d3fa35e62d45026c86fc699d Mon Sep 17 00:00:00 2001 From: maarten Date: Thu, 12 Aug 2021 19:20:29 +0200 Subject: [PATCH 06/15] 2 special cases --- src/tensors/braidingtensor.jl | 93 +++++++++++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index adc903fdd..18174560e 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -139,3 +139,96 @@ function block(b::BraidingTensor, s::Sector) # a1, a2 = f2.uncoupled # end + +function planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S}, + β, C::AbstractTensorMap{S}, + oindA::IndexTuple{2}, cindA::IndexTuple{2}, + oindB::IndexTuple{N₂}, cindB::IndexTuple, + p1::IndexTuple, p2::IndexTuple, + syms::Union{Nothing, NTuple{3, Symbol}}) where {S, N₂} + codA = codomainind(A) + domA = domainind(A) + codB = codomainind(B) + domB = domainind(B) + oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) + + braidingtensor_levels = A.adjoint ? (2,1,1,2) : (1,2,2,1); + + if iszero(β) + fill!(C, β) + elseif β != 1 + rmul!(C, β) + end + + for (f1, f2) in fusiontrees(B) + if f1 == nothing && f2 == nothing + return TO.add!(α, B,:N,true,C, reverse(cindB),oindB); + end + + fmap = Dict{Tuple{typeof(f1),typeof(f2)},eltype(f1)}(); + + #transpose + for ((f1′,f2′),coeff) in transpose(f1,f2,cindB,oindB) + + #artin braid + braid_above = braidingtensor_levels[cindA[1]] > braidingtensor_levels[cindA[2]]; + for (f1′′,coeff′) in artin_braid(f1′,1,inv = braid_above) + nk = (f1′′,f2′); + nv = coeff′*coeff; + + fmap[nk] = get(fmap,nk,zero(nv)) + nv; + end + end + + for ((f1′,f2′),c) in fmap + TO._add!(c*α, B[f1, f2], true,C[f1′,f2′], (reverse(cindB)...,oindB...)); + end + end + C +end +function planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, + β, C::AbstractTensorMap{S}, + oindA::IndexTuple{N₁}, cindA::IndexTuple{N₁}, + oindB::IndexTuple{2}, cindB::IndexTuple{2}, + p1::IndexTuple, p2::IndexTuple, + syms::Union{Nothing, NTuple{3, Symbol}}) where {S, N₁} + codA = codomainind(A) + domA = domainind(A) + codB = codomainind(B) + domB = domainind(B) + oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) + + braidingtensor_levels = B.adjoint ? (2,1,1,2) : (1,2,2,1); + + if iszero(β) + fill!(C, β) + elseif β != 1 + rmul!(C, β) + end + + for (f1, f2) in fusiontrees(A) + if f1 == nothing && f2 == nothing + return TO.add!(α, A,:N,true,C, oindA,reverse(cindA)); + end + + fmap = Dict{Tuple{typeof(f1),typeof(f2)},eltype(f1)}(); + + #transpose + for ((f1′,f2′),coeff) in transpose(f1,f2,oindA,cindA) + + #artin braid + braid_above = braidingtensor_levels[cindB[1]] > braidingtensor_levels[cindB[2]]; + for (f2′′,coeff′) in artin_braid(f2′,1,inv = braid_above) + nk = (f1′,f2′′); + nv = coeff′*coeff; + + fmap[nk] = get(fmap,nk,zero(nv)) + nv; + end + end + + for ((f1′,f2′),c) in fmap + TO._add!(c*α, A[f1, f2], true,C[f1′,f2′], (oindA...,reverse(cindA)...)); + end + end + C +end From 227ff970ca61d4b9a40a43434b0d1b86a0556d72 Mon Sep 17 00:00:00 2001 From: maarten Date: Fri, 13 Aug 2021 11:45:54 +0200 Subject: [PATCH 07/15] I don't know how to rebase --- src/fusiontrees/manipulations.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/fusiontrees/manipulations.jl b/src/fusiontrees/manipulations.jl index a73e3287f..e664a5106 100644 --- a/src/fusiontrees/manipulations.jl +++ b/src/fusiontrees/manipulations.jl @@ -522,9 +522,11 @@ const TransposeKey{I<:Sector, N₁, N₂} = Tuple{<:FusionTree{I}, <:FusionTree{ function _transpose((f1, f2, p1, p2)::TransposeKey{I,N₁,N₂}) where {I<:Sector, N₁, N₂} N = N₁ + N₂ p = linearizepermutation(p1, p2, length(f1), length(f2)) + newtrees = repartition(f1, f2, N₁) + length(p) == 0 && return newtrees i1 = findfirst(==(1), p) @assert i1 !== nothing - newtrees = repartition(f1, f2, N₁) + i1 == 1 && return newtrees Nhalf = N >> 1 while 1 < i1 <= Nhalf local newtrees′ @@ -541,7 +543,7 @@ function _transpose((f1, f2, p1, p2)::TransposeKey{I,N₁,N₂}) where {I<:Secto newtrees = newtrees′ i1 -= 1 end - while Nhalf < i1 && Nhalf > 0 + while Nhalf < i1 local newtrees′ for ((f1a, f2a), coeffa) in newtrees for ((f1b, f2b), coeffb) in cycleclockwise(f1a, f2a) From 813a6d7531ef19fe8b92234808193555460887a3 Mon Sep 17 00:00:00 2001 From: maarten Date: Fri, 13 Aug 2021 12:01:25 +0200 Subject: [PATCH 08/15] 2 more special cases --- src/tensors/braidingtensor.jl | 111 +++++++++++++++++++++++++++++++++- 1 file changed, 109 insertions(+), 2 deletions(-) diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index 18174560e..37b4611f4 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -150,7 +150,7 @@ function planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S}, domA = domainind(A) codB = codomainind(B) domB = domainind(B) - oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) + oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) braidingtensor_levels = A.adjoint ? (2,1,1,2) : (1,2,2,1); @@ -196,7 +196,7 @@ function planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, domA = domainind(A) codB = codomainind(B) domB = domainind(B) - oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) + oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) braidingtensor_levels = B.adjoint ? (2,1,1,2) : (1,2,2,1); @@ -232,3 +232,110 @@ function planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, end C end + +function TensorKit.planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S}, + β, C::AbstractTensorMap{S}, + oindA::IndexTuple{0}, cindA::IndexTuple{4}, + oindB::IndexTuple{N₂}, cindB::IndexTuple, + p1::IndexTuple, p2::IndexTuple, + syms::Union{Nothing, NTuple{3, Symbol}}) where {S, N₂} + + braidingtensor_levels = A.adjoint ? (2,1,1,2) : (1,2,2,1); + codA = codomainind(A) + domA = domainind(A) + codB = codomainind(B) + domB = domainind(B) + + oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) + + if iszero(β) + fill!(C, β) + elseif β != 1 + rmul!(C, β) + end + + for (f1, f2) in fusiontrees(B) + if f1 == nothing && f2 == nothing + TO.trace!(α, B, :N,true,C, oindB, (),(cindB[1],cindB[2]), (cindB[3],cindB[4])) + break; + end + + local fmap; + + braid_above = braidingtensor_levels[cindB[2]] > braidingtensor_levels[cindB[3]]; + + for ((f1′,f2′),coeff) in transpose(f1,f2,cindB,oindB), #transpose + (f1′′,coeff′) in artin_braid(f1′,2,inv = braid_above), #artin braid + (f1_tr1,c_tr1) in elementary_trace(f1′′, 1), #trace + (f1_tr2,c_tr2) in elementary_trace(f1_tr1,1) #trace + nk = (f1_tr2,f2′); + nv = coeff*coeff′*c_tr1*c_tr2*α; + if @isdefined fmap + fmap[nk] = get(fmap,nk,zero(nv)) + nv; + else + fmap = Dict(nk=>nv); + end + + end + + for ((f1′,f2′),c) in fmap + TO._trace!(c, B[f1, f2], true,C[f1′,f2′], oindB, (cindB[1],cindB[2]), (cindB[3],cindB[4])) + end + end + + C +end + +function TensorKit.planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, + β, C::AbstractTensorMap{S}, + oindA::IndexTuple{N₁}, cindA::IndexTuple, + oindB::IndexTuple{0}, cindB::IndexTuple{4}, + p1::IndexTuple, p2::IndexTuple, + syms::Union{Nothing, NTuple{3, Symbol}}) where {S, N₁} + + braidingtensor_levels = B.adjoint ? (2,1,1,2) : (1,2,2,1); + + codA = codomainind(A) + domA = domainind(A) + codB = codomainind(B) + domB = domainind(B) + + oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) + + + if iszero(β) + fill!(C, β) + elseif β != 1 + rmul!(C, β) + end + + braid_above = braidingtensor_levels[cindA[2]] > braidingtensor_levels[cindA[3]]; + + for (f1, f2) in fusiontrees(A) + if f1 == nothing && f2 == nothing + TO.trace!(α, A, :N,true,C, oindA , (),(cindA[1],cindA[2]), (cindA[3],cindA[4])) + break; + end + + local fmap; + for ((f1′,f2′),coeff) in transpose(f1,f2,oindA,cindA), #transpose + (f2′′,coeff′) in artin_braid(f2′,2,inv = braid_above), #artin braid + (f2_tr1,c_tr1) in elementary_trace(f2′′, 1), #trace + (f2_tr2,c_tr2) in elementary_trace(f2_tr1,1) #trace + + nk = (f1′,f2_tr2); + nv = coeff*coeff′*c_tr1*c_tr2*α; + if @isdefined fmap + fmap[nk] = get(fmap,nk,zero(nv)) + nv; + else + fmap = Dict(nk=>nv); + end + end + + for ((f1′,f2′),c) in fmap + TO._trace!(c, A[f1, f2], true,C[f1′,f2′], oindA , (cindA[1],cindA[2]), (cindA[3],cindA[4])) + end + + end + C +end From 642c43748c5aeab7ccca9e04b91d18c92190a3ba Mon Sep 17 00:00:00 2001 From: maarten Date: Fri, 13 Aug 2021 12:25:15 +0200 Subject: [PATCH 09/15] last special cases --- src/tensors/braidingtensor.jl | 101 ++++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index 37b4611f4..5f2e4d87a 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -339,3 +339,104 @@ function TensorKit.planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTens end C end + +function TensorKit.planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S}, + β, C::AbstractTensorMap{S}, + oindA::IndexTuple{1}, cindA::IndexTuple{3}, + oindB::IndexTuple{N₂}, cindB::IndexTuple, + p1::IndexTuple, p2::IndexTuple, + syms::Union{Nothing, NTuple{3, Symbol}}) where {S, N₁, N₂} + + braidingtensor_levels = A.adjoint ? (2,1,1,2) : (1,2,2,1); + + codA = codomainind(A) + domA = domainind(A) + codB = codomainind(B) + domB = domainind(B) + + oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) + + for (f1, f2) in fusiontrees(B) + if f1 == nothing && f2 == nothing + TO.trace!(α, B,:N, true,C, (cindB[2],), oindB, (cindB[1],), (cindB[3],)) + break; + end + + local fmap + braid_above = braidingtensor_levels[cindB[1]] > braidingtensor_levels[cindB[2]]; + + for ((f1′,f2′),coeff) in transpose(f1,f2,cindB,oindB), + (f1′′,coeff′) in artin_braid(f1′,1,inv = braid_above), + (f1′′′,coeff′′) in elementary_trace(f1′′, 2) + nk = (f1′′′,f2′); + nv = coeff*coeff′*coeff′′*α + if @isdefined fmap + fmap[nk] = get(fmap,nk,zero(nv)) + nv + else + fmap = Dict(nk => nv) + end + end + + for ((f1′,f2′),c) in fmap + TO._trace!(c, B[f1, f2], true,C[f1′,f2′], (cindB[2],oindB...) , (cindB[1],), (cindB[3],)) + end + end + + + C +end + +function TensorKit.planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, + β, C::AbstractTensorMap{S}, + oindA::IndexTuple{N₁}, cindA::IndexTuple, + oindB::IndexTuple{1}, cindB::IndexTuple{3}, + p1::IndexTuple, p2::IndexTuple, + syms::Union{Nothing, NTuple{3, Symbol}}) where {S, N₁} + braidingtensor_levels = B.adjoint ? (2,1,1,2) : (1,2,2,1); + + codA = codomainind(A) + domA = domainind(A) + codB = codomainind(B) + domB = domainind(B) + + oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) + + if iszero(β) + fill!(C, β) + elseif β != 1 + rmul!(C, β) + end + + for (f1, f2) in fusiontrees(A) + if f1 == nothing && f2 == nothing + + TO.trace!(α, A,:N, true,C, oindA, (cindA[2],) , (cindA[1],), (cindA[3],)) + break; + end + + braid_above = braidingtensor_levels[cindA[1]] > braidingtensor_levels[cindA[2]]; + + local fmap + + for ((f1′,f2′),coeff) in transpose(f1,f2,cindA,oindA), + (f2′′,coeff′) in artin_braid(f2′,1,inv = braid_above), + (f2′′′,coeff′′) in elementary_trace(f2′′, 2) + + nk = (f1′,f2′′′); + nv = coeff*coeff′*coeff′′*α; + + if @isdefined fmap + fmap[nk] = get(fmap,nk,zero(nv)) + nv + else + famp = Dict(nk => nv); + end + + end + + for ((f1′,f2′),c) in fmap + TO._trace!(c, A[f1, f2], true,C[f1′,f2′], (oindA...,cindA[2]) , (cindA[1],), (cindA[3],)) + end + end + + C +end From 766f1038e177e9b68f0e1b54e8ea1d6719b639a6 Mon Sep 17 00:00:00 2001 From: maarten Date: Fri, 13 Aug 2021 13:40:44 +0200 Subject: [PATCH 10/15] bugfix --- src/tensors/braidingtensor.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index 5f2e4d87a..c719f4e18 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -262,7 +262,7 @@ function TensorKit.planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{ local fmap; - braid_above = braidingtensor_levels[cindB[2]] > braidingtensor_levels[cindB[3]]; + braid_above = braidingtensor_levels[cindA[2]] > braidingtensor_levels[cindA[3]]; for ((f1′,f2′),coeff) in transpose(f1,f2,cindB,oindB), #transpose (f1′′,coeff′) in artin_braid(f1′,2,inv = braid_above), #artin braid @@ -309,7 +309,7 @@ function TensorKit.planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTens rmul!(C, β) end - braid_above = braidingtensor_levels[cindA[2]] > braidingtensor_levels[cindA[3]]; + braid_above = braidingtensor_levels[cindB[2]] > braidingtensor_levels[cindB[3]]; for (f1, f2) in fusiontrees(A) if f1 == nothing && f2 == nothing @@ -363,7 +363,7 @@ function TensorKit.planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{ end local fmap - braid_above = braidingtensor_levels[cindB[1]] > braidingtensor_levels[cindB[2]]; + braid_above = braidingtensor_levels[cindA[1]] > braidingtensor_levels[cindA[2]]; for ((f1′,f2′),coeff) in transpose(f1,f2,cindB,oindB), (f1′′,coeff′) in artin_braid(f1′,1,inv = braid_above), @@ -414,11 +414,11 @@ function TensorKit.planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTens break; end - braid_above = braidingtensor_levels[cindA[1]] > braidingtensor_levels[cindA[2]]; + braid_above = braidingtensor_levels[cindB[1]] > braidingtensor_levels[cindB[2]]; local fmap - for ((f1′,f2′),coeff) in transpose(f1,f2,cindA,oindA), + for ((f1′,f2′),coeff) in transpose(f1,f2,oindA,cindA), (f2′′,coeff′) in artin_braid(f2′,1,inv = braid_above), (f2′′′,coeff′′) in elementary_trace(f2′′, 2) @@ -428,7 +428,7 @@ function TensorKit.planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTens if @isdefined fmap fmap[nk] = get(fmap,nk,zero(nv)) + nv else - famp = Dict(nk => nv); + fmap = Dict(nk => nv); end end From 2c3e12b64235f2dc919cd8aa70bb673f84ef39c9 Mon Sep 17 00:00:00 2001 From: maarten Date: Fri, 13 Aug 2021 13:59:49 +0200 Subject: [PATCH 11/15] touchup --- src/tensors/braidingtensor.jl | 68 +++++++++++++++++------------------ 1 file changed, 32 insertions(+), 36 deletions(-) diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index c719f4e18..91883d18a 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -152,7 +152,7 @@ function planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S}, domB = domainind(B) oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) - braidingtensor_levels = A.adjoint ? (2,1,1,2) : (1,2,2,1); + braidingtensor_levels = A.adjoint ? (1,2,2,1) : (2,1,1,2) if iszero(β) fill!(C, β) @@ -167,17 +167,15 @@ function planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S}, fmap = Dict{Tuple{typeof(f1),typeof(f2)},eltype(f1)}(); - #transpose - for ((f1′,f2′),coeff) in transpose(f1,f2,cindB,oindB) - #artin braid - braid_above = braidingtensor_levels[cindA[1]] > braidingtensor_levels[cindA[2]]; - for (f1′′,coeff′) in artin_braid(f1′,1,inv = braid_above) - nk = (f1′′,f2′); - nv = coeff′*coeff; + braid_above = braidingtensor_levels[cindA[1]] > braidingtensor_levels[cindA[2]]; + for ((f1′,f2′),coeff) in transpose(f1,f2,cindB,oindB), + (f1′′,coeff′) in artin_braid(f1′,1,inv = braid_above) - fmap[nk] = get(fmap,nk,zero(nv)) + nv; - end + nk = (f1′′,f2′); + nv = coeff′*coeff; + + fmap[nk] = get(fmap,nk,zero(nv)) + nv; end for ((f1′,f2′),c) in fmap @@ -198,7 +196,7 @@ function planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, domB = domainind(B) oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2) - braidingtensor_levels = B.adjoint ? (2,1,1,2) : (1,2,2,1); + braidingtensor_levels = B.adjoint ? (1,2,2,1) : (2,1,1,2); if iszero(β) fill!(C, β) @@ -213,17 +211,15 @@ function planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, fmap = Dict{Tuple{typeof(f1),typeof(f2)},eltype(f1)}(); - #transpose - for ((f1′,f2′),coeff) in transpose(f1,f2,oindA,cindA) + braid_above = braidingtensor_levels[cindB[1]] > braidingtensor_levels[cindB[2]]; - #artin braid - braid_above = braidingtensor_levels[cindB[1]] > braidingtensor_levels[cindB[2]]; - for (f2′′,coeff′) in artin_braid(f2′,1,inv = braid_above) - nk = (f1′,f2′′); - nv = coeff′*coeff; + for ((f1′,f2′),coeff) in transpose(f1,f2,oindA,cindA), + (f2′′,coeff′) in artin_braid(f2′,1,inv = braid_above) - fmap[nk] = get(fmap,nk,zero(nv)) + nv; - end + nk = (f1′,f2′′); + nv = coeff′*coeff; + + fmap[nk] = get(fmap,nk,zero(nv)) + nv; end for ((f1′,f2′),c) in fmap @@ -233,14 +229,14 @@ function planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, C end -function TensorKit.planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S}, +function planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S}, β, C::AbstractTensorMap{S}, oindA::IndexTuple{0}, cindA::IndexTuple{4}, oindB::IndexTuple{N₂}, cindB::IndexTuple, p1::IndexTuple, p2::IndexTuple, syms::Union{Nothing, NTuple{3, Symbol}}) where {S, N₂} - braidingtensor_levels = A.adjoint ? (2,1,1,2) : (1,2,2,1); + braidingtensor_levels = A.adjoint ? (1,2,2,1) : (2,1,1,2); codA = codomainind(A) domA = domainind(A) codB = codomainind(B) @@ -264,10 +260,10 @@ function TensorKit.planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{ braid_above = braidingtensor_levels[cindA[2]] > braidingtensor_levels[cindA[3]]; - for ((f1′,f2′),coeff) in transpose(f1,f2,cindB,oindB), #transpose - (f1′′,coeff′) in artin_braid(f1′,2,inv = braid_above), #artin braid - (f1_tr1,c_tr1) in elementary_trace(f1′′, 1), #trace - (f1_tr2,c_tr2) in elementary_trace(f1_tr1,1) #trace + for ((f1′,f2′),coeff) in transpose(f1,f2,cindB,oindB), + (f1′′,coeff′) in artin_braid(f1′,2,inv = braid_above), + (f1_tr1,c_tr1) in elementary_trace(f1′′, 1), + (f1_tr2,c_tr2) in elementary_trace(f1_tr1,1) nk = (f1_tr2,f2′); nv = coeff*coeff′*c_tr1*c_tr2*α; if @isdefined fmap @@ -286,14 +282,14 @@ function TensorKit.planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{ C end -function TensorKit.planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, +function planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, β, C::AbstractTensorMap{S}, oindA::IndexTuple{N₁}, cindA::IndexTuple, oindB::IndexTuple{0}, cindB::IndexTuple{4}, p1::IndexTuple, p2::IndexTuple, syms::Union{Nothing, NTuple{3, Symbol}}) where {S, N₁} - braidingtensor_levels = B.adjoint ? (2,1,1,2) : (1,2,2,1); + braidingtensor_levels = B.adjoint ? (1,2,2,1) : (2,1,1,2); codA = codomainind(A) domA = domainind(A) @@ -318,10 +314,10 @@ function TensorKit.planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTens end local fmap; - for ((f1′,f2′),coeff) in transpose(f1,f2,oindA,cindA), #transpose - (f2′′,coeff′) in artin_braid(f2′,2,inv = braid_above), #artin braid - (f2_tr1,c_tr1) in elementary_trace(f2′′, 1), #trace - (f2_tr2,c_tr2) in elementary_trace(f2_tr1,1) #trace + for ((f1′,f2′),coeff) in transpose(f1,f2,oindA,cindA), + (f2′′,coeff′) in artin_braid(f2′,2,inv = braid_above), + (f2_tr1,c_tr1) in elementary_trace(f2′′, 1), + (f2_tr2,c_tr2) in elementary_trace(f2_tr1,1) nk = (f1′,f2_tr2); nv = coeff*coeff′*c_tr1*c_tr2*α; @@ -340,14 +336,14 @@ function TensorKit.planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTens C end -function TensorKit.planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S}, +function planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S}, β, C::AbstractTensorMap{S}, oindA::IndexTuple{1}, cindA::IndexTuple{3}, oindB::IndexTuple{N₂}, cindB::IndexTuple, p1::IndexTuple, p2::IndexTuple, syms::Union{Nothing, NTuple{3, Symbol}}) where {S, N₁, N₂} - braidingtensor_levels = A.adjoint ? (2,1,1,2) : (1,2,2,1); + braidingtensor_levels = A.adjoint ? (1,2,2,1) : (2,1,1,2); codA = codomainind(A) domA = domainind(A) @@ -386,13 +382,13 @@ function TensorKit.planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{ C end -function TensorKit.planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, +function planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor, β, C::AbstractTensorMap{S}, oindA::IndexTuple{N₁}, cindA::IndexTuple, oindB::IndexTuple{1}, cindB::IndexTuple{3}, p1::IndexTuple, p2::IndexTuple, syms::Union{Nothing, NTuple{3, Symbol}}) where {S, N₁} - braidingtensor_levels = B.adjoint ? (2,1,1,2) : (1,2,2,1); + braidingtensor_levels = B.adjoint ? (1,2,2,1) : (2,1,1,2); codA = codomainind(A) domA = domainind(A) From 030b7766d13269efd70c9f25b16932024b1bafe8 Mon Sep 17 00:00:00 2001 From: maarten Date: Fri, 13 Aug 2021 14:37:43 +0200 Subject: [PATCH 12/15] artin_braid inference failure --- src/fusiontrees/manipulations.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/fusiontrees/manipulations.jl b/src/fusiontrees/manipulations.jl index e664a5106..1b0e26ab0 100644 --- a/src/fusiontrees/manipulations.jl +++ b/src/fusiontrees/manipulations.jl @@ -769,6 +769,12 @@ function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I<:Sector vertices = f.vertices u = one(I) + if BraidingStyle(I) isa NoBraiding + oneT = one(eltype(Fsymbol(u,u,u,u,u,u))) + else + oneT = one(eltype(Rsymbol(u,u,u))) * one(eltype(Fsymbol(u,u,u,u,u,u))) + end + if u in (uncoupled[i],uncoupled[i+1]) # the braid simplifies drastically, we are braiding a trivial sector a, b = uncoupled[i], uncoupled[i+1] uncoupled′ = TupleTools.setindex(uncoupled,b,i); @@ -789,10 +795,10 @@ function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I<:Sector end f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices′) - return fusiontreedict(I)(f′ => true) + return fusiontreedict(I)(f′ => oneT) end - oneT = one(eltype(Rsymbol(u,u,u))) * one(eltype(Fsymbol(u,u,u,u,u,u))) + BraidingStyle(I) isa NoBraiding && throw(SectorMismatch("cannot braid sector "*type_repr(I))) if i == 1 a, b = uncoupled[1], uncoupled[2] From 053160e8da22368e3c7349e0df8adb0158d39444 Mon Sep 17 00:00:00 2001 From: maarten Date: Fri, 13 Aug 2021 14:38:34 +0200 Subject: [PATCH 13/15] _one! --- src/tensors/braidingtensor.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index 91883d18a..4422f4672 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -108,7 +108,7 @@ function Base.copy!(t::TensorMap, b::BraidingTensor) fill!(t, zero(eltype(t))) for (f1, f2) in fusiontrees(t) if f1 == nothing || f2 == nothing - copy!(t.data,one(t.data)) + _one!(t.data) return t end From ccd0482ae09b57ac5678dadafcadbc68d771c3de Mon Sep 17 00:00:00 2001 From: maarten Date: Fri, 13 Aug 2021 15:01:54 +0200 Subject: [PATCH 14/15] possible workaround for planars? --- src/tensors/planar.jl | 30 +++++++++++++++++++++++++----- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/src/tensors/planar.jl b/src/tensors/planar.jl index df11a4626..6d48f03e4 100644 --- a/src/tensors/planar.jl +++ b/src/tensors/planar.jl @@ -67,8 +67,10 @@ end _conj_to_adjoint(ex) = ex function get_possible_planar_indices(ex::Expr) - @assert TO.istensorexpr(ex) - if TO.isgeneraltensor(ex) + #@assert TO.istensorexpr(ex) + if !TO.istensorexpr(ex) + return [[]] + elseif TO.isgeneraltensor(ex) _,leftind,rightind = TO.decomposegeneraltensor(ex) ind = planar_unique2(vcat(leftind, reverse(rightind))) return length(ind) == length(unique(ind)) ? Any[ind] : Any[] @@ -270,6 +272,17 @@ function _extract_contraction_pairs(rhs, lhs, pre, temporaries) a1 = _extract_contraction_pairs(rhs.args[2], (oind1, reverse(cind1)), pre, temporaries) a2 = _extract_contraction_pairs(rhs.args[3], (cind2, reverse(oind2)), pre, temporaries) end + + if TO.isscalarexpr(a1) || TO.isscalarexpr(a2) + rhs = Expr(:call, :*, a1, a2) + s = gensym() + newlhs = Expr(:typed_vcat, s, Expr(:tuple, oind1...), + Expr(:tuple, reverse(oind2)...)) + push!(temporaries, s) + push!(pre, Expr(:(:=), newlhs, rhs)) + return newlhs + end + # note that index order in _extract... is only a suggestion, now we have actual index order _, l1, r1, = TO.decomposegeneraltensor(a1) _, l2, r2, = TO.decomposegeneraltensor(a2) @@ -307,6 +320,9 @@ function _extract_contraction_pairs(rhs, lhs, pre, temporaries) args = [_extract_contraction_pairs(a, lhs, pre, temporaries) for a in rhs.args[2:end]] return Expr(rhs.head, rhs.args[1], args...) + elseif TO.isscalarexpr(rhs) + #do nothing? + return rhs else throw(ArgumentError("unknown tensor expression")) end @@ -500,11 +516,15 @@ function _update_temporaries(ex, temporaries) i = findfirst(==(lhs), temporaries) if i !== nothing rhs = ex.args[2] - if !(rhs isa Expr && rhs.head == :call && rhs.args[1] == :contract!) + if rhs isa Expr && rhs.head == :call && rhs.args[1] == :add! + newname = rhs.args[6] + temporaries[i] = newname + elseif rhs isa Expr && rhs.head == :call && rhs.args[1] == :contract! + newname = rhs.args[8] + temporaries[i] = newname + else @error "lhs = $lhs , rhs = $rhs" end - newname = rhs.args[8] - temporaries[i] = newname end elseif ex isa Expr for a in ex.args From b5e7ef5bb77baf009f5dfbc6e9b400c51a970980 Mon Sep 17 00:00:00 2001 From: maarten Date: Sat, 14 Aug 2021 13:40:17 +0200 Subject: [PATCH 15/15] proposed solution to #58 --- src/tensors/planar.jl | 80 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 64 insertions(+), 16 deletions(-) diff --git a/src/tensors/planar.jl b/src/tensors/planar.jl index 6d48f03e4..7ffbaa91b 100644 --- a/src/tensors/planar.jl +++ b/src/tensors/planar.jl @@ -450,8 +450,10 @@ function _construct_braidingtensors(ex::Expr) i += 1 end end - pre = Expr(:block) - for (t, construct_expr) in translatebraidings + + unresolved = Any[]; # list of indices that we couldn't yet figure out + indexmaps = Dict{Any,Any}(); # contains the expression to resolve the space at index indexmaps[i] + for (t,construct_expr) in translatebraidings obj, leftind, rightind = TO.decomposetensor(t) length(leftind) == length(rightind) == 2 || throw(ArgumentError("The name τ is reserved for the braiding, and should have two input and two output indices.")) @@ -462,25 +464,71 @@ function _construct_braidingtensors(ex::Expr) i2b, i1b, = leftind i1a, i2a, = rightind end - obj_and_pos = _findindex(i1a, list) - if !isnothing(obj_and_pos) - push!(construct_expr.args, Expr(:call, :space, obj_and_pos...)) + + obj_and_pos1a = _findindex(i1a, list) + obj_and_pos2a = _findindex(i2a, list) + obj_and_pos1b = _findindex(i1b, list) + obj_and_pos2b = _findindex(i2b, list) + + if !isnothing(obj_and_pos1a) + indexmaps[i1a] = Expr(:call, :space, obj_and_pos1a...); + indexmaps[i1b] = Expr(TO.prime, Expr(:call, :space, obj_and_pos1a...)); + elseif !isnothing(obj_and_pos1b) + indexmaps[i1a] = Expr(TO.prime, Expr(:call, :space, obj_and_pos1b...)); + indexmaps[i1b] = Expr(:call, :space, obj_and_pos1b...); else - obj_and_pos = _findindex(i1b, list) - isnothing(obj_and_pos) && - throw(ArgumentError("cannot determine space of index $i1a and $i1b of braiding tensor")) - push!(construct_expr.args, Expr(TO.prime, Expr(:call, :space, obj_and_pos...))) + push!(unresolved,(i1a,i1b)); end - obj_and_pos = _findindex(i2a, list) - if !isnothing(obj_and_pos) - push!(construct_expr.args, Expr(:call, :space, obj_and_pos...)) + if !isnothing(obj_and_pos2a) + indexmaps[i2a] = Expr(:call, :space, obj_and_pos2a...); + indexmaps[i2b] = Expr(TO.prime, Expr(:call, :space, obj_and_pos2a...)) + elseif !isnothing(obj_and_pos2b) + indexmaps[i2a] = Expr(TO.prime, Expr(:call, :space, obj_and_pos2b...)) + indexmaps[i2b] = Expr(:call, :space, obj_and_pos2b...); + else + push!(unresolved,(i2a,i2b)); + end + end + + # very simple loop that tries to resolve as many indices as possible + changed = true; + while changed == true + changed = false; + + i = 1; + while i<=length(unresolved) + (i1,i2) = unresolved[i]; + if i1 in keys(indexmaps) + changed = true; + indexmaps[i2] = Expr(TO.prime,indexmaps[i1]); + deleteat!(unresolved,i); + elseif i2 in keys(indexmaps) + changed = true; + indexmaps[i1] = Expr(TO.prime,indexmaps[i2]); + deleteat!(unresolved,i); + else + i+=1 + end + end + end + + !isempty(unresolved) && throw(ArgumentError("cannot determine the spaces of indices $(tuple(unresolved...)) for the braiding tensors in $(ex)")); + + pre = Expr(:block) + for (t, construct_expr) in translatebraidings + obj, leftind, rightind = TO.decomposetensor(t) + if _is_adjoint(obj) + i1b, i2b, = leftind + i2a, i1a, = rightind else - obj_and_pos = _findindex(i2b, list) - isnothing(obj_and_pos) && - throw(ArgumentError("cannot determine space of index $i2a and $i2b of braiding tensor")) - push!(construct_expr.args, Expr(TO.prime, Expr(:call, :space, obj_and_pos...))) + i2b, i1b, = leftind + i1a, i2a, = rightind end + + push!(construct_expr.args, indexmaps[i1a]); + push!(construct_expr.args, indexmaps[i2a]); + s = gensym() push!(pre.args, Expr(:(=), s, construct_expr)) ex = TO.replacetensorobjects(ex) do o, l, r