explicit flip of domain tree

Author: ogauthe

src/fusion_trees/fusiontree.jl

@@ -4,7 +4,7 @@
 # compatibility with TensorKit conventions?
 
 using GradedUnitRanges:
-  AbstractGradedUnitRange, GradedUnitRanges, fusion_product, isdual, sector_type
+  AbstractGradedUnitRange, GradedUnitRanges, flip, fusion_product, isdual, sector_type
 using SymmetrySectors: ×, AbstractSector, SectorProduct, SymmetrySectors, arguments, trivial
 using TensorAlgebra: flatten_tuples
 
@@ -80,6 +80,15 @@ Base.length(::SectorFusionTree{<:Any,N}) where {N} = N
 
 # GradedUnitRanges interface
 GradedUnitRanges.sector_type(::Type{<:SectorFusionTree{S}}) where {S} = S
+function GradedUnitRanges.flip(f::SectorFusionTree)
+  return SectorFusionTree(
+    dual.(leaves(f)),
+    .!arrows(f),
+    dual(root_sector(f)),
+    dual.(branch_sectors(f)),
+    outer_multiplicity_indices(f),
+  )
+end
 
 # SymmetrySectors interface
 function SymmetrySectors.:×(f1::SectorFusionTree, f2::SectorFusionTree)

src/fusiontensor/array_cast.jl

@@ -165,7 +165,7 @@ end
 
 function contract_singlet_projector(f1::SectorFusionTree, f2::SectorFusionTree)
   f1_array = convert(Array, f1)
-  f2_array = convert(Array, f2)
+  f2_array = convert(Array, flip(f2))
   N_CO = length(f1)
   N_DO = length(f2)
   return contract(

src/fusiontensor/base_interface.jl

@@ -38,7 +38,9 @@ Base.Array(ft::FusionTensor) = Array(to_array(ft))
 
 # adjoint is costless: dual axes, swap codomain and domain, take data_matrix adjoint.
 # data_matrix coeff are not modified (beyond complex conjugation)
-transpose_mapping(d::Dict) = Dict([reverse(k) => transpose_mapping(v) for (k, v) in d])
+function transpose_mapping(d::Dict)
+  return Dict([reverse(flip.(k)) => transpose_mapping(v) for (k, v) in d])
+end
 function transpose_mapping(b::BlockIndexRange{2})
   new_block = Block(reverse(Tuple(Block(b))))
   return new_block[reverse(b.indices)...]

src/fusiontensor/fusiontensor.jl

@@ -9,6 +9,7 @@ using GradedUnitRanges:
   blocklabels,
   blockmergesort,
   dual,
+  flip,
   gradedrange,
   isdual,
   map_blocklabels,
@@ -71,7 +72,7 @@ function BlockArrays.findblock(ft::FusionTensor, f1::SectorFusionTree, f2::Secto
   # find outer block corresponding to fusion trees
   @assert typeof((f1, f2)) === keytype(trees_block_mapping(ft))
   b1 = find_sector_block.(leaves(f1), codomain_axes(ft))
-  b2 = find_sector_block.(leaves(f2), dual.(domain_axes(ft)))
+  b2 = find_sector_block.(leaves(f2), domain_axes(ft))
   return Block(b1..., b2...)
 end
 # TBD move to GradedUnitRanges? rename findfirst_sector?
@@ -137,16 +138,24 @@ function FusionTensor(mat::AbstractMatrix, legs::BlockedTuple{2})
   return ft
 end
 
+function flip_domain(nonflipped_col_axis, nonflipped_trees_to_ranges)
+  col_axis = dual(nonflipped_col_axis)
+  domain_trees_to_ranges_mapping = Dict(
+    map(((tree, v),) -> flip(tree) => v, collect(nonflipped_trees_to_ranges))
+  )
+  return col_axis, domain_trees_to_ranges_mapping
+end
+
 # empty matrix
 function FusionTensor(elt::Type, raw_legs::BlockedTuple{2})
   S, legs = sanitize_axes(raw_legs)
 
-  row_axis, codomain_trees_to_ranges_mapping = fuse_axes(S, first(blocks(legs)))
-  nondual_col_axis, domain_trees_to_ranges_mapping = fuse_axes(S, dual.(last(blocks(legs))))
+  row_axis, codomain_trees_to_ranges = fuse_axes(S, first(blocks(legs)))
+  col_axis, domain_trees_to_ranges = flip_domain(fuse_axes(S, dual.(last(blocks(legs))))...)
 
-  mat = initialize_data_matrix(elt, row_axis, nondual_col_axis)
+  mat = initialize_data_matrix(elt, row_axis, col_axis)
   tree_to_block_mapping = intersect_codomain_domain(
-    codomain_trees_to_ranges_mapping, domain_trees_to_ranges_mapping
+    codomain_trees_to_ranges, domain_trees_to_ranges
   )
   return FusionTensor(mat, legs, tree_to_block_mapping)
 end
@@ -170,9 +179,7 @@ function fusion_trees_external_multiplicities(
   )
 end
 
-function block_fusion_trees_external_multiplicities(
-  it::NTuple{N,AbstractUnitRange}
-) where {N}
+function block_fusion_trees_external_multiplicities(it::Tuple{Vararg{AbstractUnitRange}})
   block_sectors = only.(blocklabels.(it))
   block_mult = prod(length.(it))
   return build_trees(block_sectors, isdual.(it)) .=> block_mult
@@ -206,7 +213,7 @@ function intersect_codomain_domain(
   return Dict(
     map(
       Iterators.filter(
-        t -> root_sector(first(t[1])) == root_sector(first(t[2])),
+        t -> root_sector(first(first(t))) == dual(root_sector(first(t[2]))),
         Iterators.product(codomain_trees_to_ranges_mapping, domain_trees_to_ranges_mapping),
       ),
     ) do t
@@ -218,11 +225,11 @@ end
 function initialize_data_matrix(
   elt::Type{<:Number},
   mat_row_axis::AbstractGradedUnitRange,
-  nondual_col_axis::AbstractGradedUnitRange,
+  mat_col_axis::AbstractGradedUnitRange,
 )
   # non-abelian fusion trees have float eltype: need compatible type
   promoted = promote_type(elt, fusiontree_eltype(sector_type(mat_row_axis)))
-  mat = BlockSparseArray{promoted}(mat_row_axis, dual(nondual_col_axis))
+  mat = BlockSparseArray{promoted}(mat_row_axis, mat_col_axis)
   initialize_allowed_sectors!(mat)
   return mat
 end