First draft of FusionTensor

Project.toml

@@ -1,7 +1,37 @@
 name = "FusionTensors"
 uuid = "e16ca583-1f51-4df0-8e12-57d32947d33e"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.1.0"
+version = "0.2.0"
+
+[deps]
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
+BroadcastMapConversion = "4a4adec5-520f-4750-bb37-d5e66b4ddeb2"
+GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5"
+HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721"
+LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
+LabelledNumbers = "f856a3a6-4152-4ec4-b2a7-02c1a55d7993"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
+Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
+SymmetrySectors = "f8a8ad64-adbc-4fce-92f7-ffe2bb36a86e"
+TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
+TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
+WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
 
 [compat]
+BlockArrays = "1.2.0"
+BlockSparseArrays = "0.2.0"
+BroadcastMapConversion = "0.1.0"
+GradedUnitRanges = "0.1.1"
+HalfIntegers = "1.6.0"
+LRUCache = "1.6.1"
+LabelledNumbers = "0.1.0"
+LinearAlgebra = "1.11.0"
+SparseArraysBase = "0.2.0"
+Strided = "2.2.0"
+SymmetrySectors = "0.1.1"
+TensorAlgebra = "0.1.0"
+TypeParameterAccessors = "0.1.0"
+WignerSymbols = "2.0.0"
 julia = "1.10"

src/FusionTensors.jl

@@ -1,5 +1,14 @@
 module FusionTensors
 
-# Write your package code here.
+include("fusion_trees/fusiontree.jl")
+include("fusion_trees/clebsch_gordan_tensors.jl")
+include("fusiontensor/fusedaxes.jl")
+include("fusiontensor/fusiontensor.jl")
+include("fusiontensor/base_interface.jl")
+include("fusiontensor/array_cast.jl")
+include("fusiontensor/linear_algebra_interface.jl")
+include("fusiontensor/tensor_algebra_interface.jl")
+include("permutedims/unitaries.jl")
+include("permutedims/permutedims.jl")
 
 end

src/fusion_trees/clebsch_gordan_tensors.jl

@@ -0,0 +1,136 @@
+# This file defines Clebsch-Gordan tensors
+# one tensor is defined from 3 simple objects s1, s2 and s3
+# and contains the coefficients fusing s1 ⊗ s2 -> s3
+
+using HalfIntegers: half
+using WignerSymbols: clebschgordan
+
+using GradedUnitRanges: dual, fusion_product
+using SymmetrySectors:
+  ⊗,
+  AbelianStyle,
+  AbstractSector,
+  NotAbelianStyle,
+  SymmetryStyle,
+  O2,
+  SU,
+  istrivial,
+  quantum_dimension,
+  sector_label,
+  trivial,
+  zero_odd
+using TensorAlgebra: contract
+
+function symbol_1j(s::AbstractSector)
+  cgt = clebsch_gordan_tensor(s, dual(s), trivial(s), 1)
+  return sqrt(quantum_dimension(s)) * cgt[:, :, 1]
+end
+
+function clebsch_gordan_tensor(
+  s1::AbstractSector,
+  s2::AbstractSector,
+  s3::AbstractSector,
+  arrow1::Bool,
+  arrow2::Bool,
+  inner_mult_index::Int,
+)
+  cgt = clebsch_gordan_tensor(s1, s2, s3, inner_mult_index)
+  if arrow1
+    flip1 = symbol_1j(s1)
+    cgt = contract((1, 2, 3), flip1, (4, 1), cgt, (4, 2, 3))
+  end
+  if arrow2
+    flip2 = symbol_1j(s2)
+    cgt = contract((1, 2, 3), flip2, (4, 2), cgt, (1, 4, 3))
+  end
+  return cgt
+end
+
+function clebsch_gordan_tensor(s1::S, s2::S, s3::S, outer_mult_index::Int) where {S}
+  return clebsch_gordan_tensor(SymmetryStyle(S), s1, s2, s3, outer_mult_index)
+end
+
+function clebsch_gordan_tensor(
+  ::AbelianStyle, s1::S, s2::S, s3::S, outer_mult_index::Int
+) where {S}
+  @assert outer_mult_index == 1
+  return s1 ⊗ s2 == s3 ? ones((1, 1, 1)) : zeros((1, 1, 1))
+end
+
+function clebsch_gordan_tensor(::NotAbelianStyle, s1::O2, s2::O2, s3::O2, ::Int)
+  return clebsch_gordan_tensor(s1, s2, s3)  # no outer multiplicity
+end
+
+function clebsch_gordan_tensor(s1::O2, s2::O2, s3::O2)
+  d1 = quantum_dimension(s1)
+  d2 = quantum_dimension(s2)
+  d3 = quantum_dimension(s3)
+  cgt = zeros((d1, d2, d3))
+  s3 ∉ blocklabels(fusion_product(s1, s2)) && return cgt
+
+  # adapted from TensorKit
+  l1 = sector_label(s1)
+  l2 = sector_label(s2)
+  l3 = sector_label(s3)
+  if l3 <= 0  # 0even or 0odd
+    if l1 <= 0 && l2 <= 0
+      cgt[1, 1, 1, 1] = 1.0
+    else
+      if istrivial(s3)
+        cgt[1, 2, 1, 1] = 1.0 / sqrt(2)
+        cgt[2, 1, 1, 1] = 1.0 / sqrt(2)
+      else
+        cgt[1, 2, 1, 1] = 1.0 / sqrt(2)
+        cgt[2, 1, 1, 1] = -1.0 / sqrt(2)
+      end
+    end
+  elseif l1 <= 0  # 0even or 0odd
+    cgt[1, 1, 1, 1] = 1.0
+    cgt[1, 2, 2, 1] = s1 == zero_odd(O2) ? -1.0 : 1.0
+  elseif l2 == 0
+    cgt[1, 1, 1, 1] = 1.0
+    cgt[2, 1, 2, 1] = s2 == zero_odd(O2) ? -1.0 : 1.0
+  elseif l3 == l1 + l2
+    cgt[1, 1, 1, 1] = 1.0
+    cgt[2, 2, 2, 1] = 1.0
+  elseif l3 == l1 - l2
+    cgt[1, 2, 1, 1] = 1.0
+    cgt[2, 1, 2, 1] = 1.0
+  elseif l3 == l2 - l1
+    cgt[2, 1, 1, 1] = 1.0
+    cgt[1, 2, 2, 1] = 1.0
+  end
+  return cgt
+end
+
+function clebsch_gordan_tensor(::NotAbelianStyle, s1::SU{2}, s2::SU{2}, s3::SU{2}, ::Int)
+  return clebsch_gordan_tensor(s1, s2, s3)  # no outer multiplicity
+end
+
+function clebsch_gordan_tensor(s1::SU{2}, s2::SU{2}, s3::SU{2})
+  d1 = quantum_dimension(s1)
+  d2 = quantum_dimension(s2)
+  d3 = quantum_dimension(s3)
+  j1 = half(d1 - 1)
+  j2 = half(d2 - 1)
+  j3 = half(d3 - 1)
+  cgtensor = Array{Float64,3}(undef, (d1, d2, d3))
+  for (i, j, k) in Iterators.product(1:d1, 1:d2, 1:d3)
+    m1 = j1 - i + 1
+    m2 = j2 - j + 1
+    m3 = j3 - k + 1
+    cgtensor[i, j, k] = clebschgordan(j1, m1, j2, m2, j3, m3)
+  end
+  return cgtensor
+end
+
+function clebsch_gordan_tensor(
+  ::NotAbelianStyle, s1::SU{3}, s2::SU{3}, s3::SU{3}, outer_mult_index::Int
+)
+  d1 = quantum_dimension(s1)
+  d2 = quantum_dimension(s2)
+  d3 = quantum_dimension(s3)
+  cgtensor = zeros(d1, d2, d3)
+  # dummy
+  return cgtensor
+end