diff --git a/.github/workflows/CI-Sectors.yml b/.github/workflows/CI-Sectors.yml deleted file mode 100644 index bb3ca81a..00000000 --- a/.github/workflows/CI-Sectors.yml +++ /dev/null @@ -1,53 +0,0 @@ -name: CI - TensorKitSectors -on: - push: - branches: - - 'master' - - 'main' - - 'release-' - paths: - - TensorKitSectors/ - tags: '*' - pull_request: - paths: - - TensorKitSectors/ - workflow_dispatch: - -concurrency: - group: ${{ github.workflow }}-${{ github.ref }} - cancel-in-progress: true - -jobs: - test: - name: TensorKitSectors - ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - runs-on: ${{ matrix.os }} - strategy: - fail-fast: false - matrix: - version: - - '1' # automatically expands to the latest stable 1.x release of Julia - os: - - ubuntu-latest - arch: - - x64 - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 - - uses: julia-actions/julia-buildpkg@latest - with: - project: TensorKitSectors - - uses: julia-actions/julia-runtest@latest - with: - project: TensorKitSectors - env: - JULIA_NUM_THREADS: 4 - - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v4 - env: - CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} - with: - file: lcov.info \ No newline at end of file diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 9959b2f9..8b2c5c73 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -21,8 +21,8 @@ jobs: fail-fast: false matrix: version: - - '1.8' # LTS version - - '1' # automatically expands to the latest stable 1.x release of Julia + - '1.10' # LTS version + # - '1' # TODO: add back when 1.10 is LTS os: - ubuntu-latest - macOS-latest diff --git a/Changelog.md b/Changelog.md index cd26928c..2c7979d8 100644 --- a/Changelog.md +++ b/Changelog.md @@ -42,3 +42,8 @@ uninitialized tensors is now `TensorMap{E}(undef, codomain ← domain)`, reminis This PR bumps the compatibility of `TensorOperations` to v5. This is a breaking change as there are some changes in the API. + +### TensorKitSectors + +This promotes TensorKitSectors to its own package, in order to make the dependencies +lighter and to separate the concerns of the two packages. diff --git a/Project.toml b/Project.toml index 4990f8fd..fc0ecf23 100644 --- a/Project.toml +++ b/Project.toml @@ -4,16 +4,15 @@ authors = ["Jutho Haegeman"] version = "1.0.0-DEV" [deps] -HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" +TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" -WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" [weakdeps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" @@ -29,19 +28,18 @@ ChainRulesCore = "1" ChainRulesTestUtils = "1" Combinatorics = "1" FiniteDifferences = "0.12" -HalfIntegers = "1" LRUCache = "1.0.2" LinearAlgebra = "1" PackageExtensionCompat = "1" Random = "1" Strided = "2" +TensorKitSectors = "0.1" TensorOperations = "5" Test = "1" TestExtras = "0.2" TupleTools = "1.1" VectorInterface = "0.4" -WignerSymbols = "1,2" -julia = "1.6" +julia = "1.10" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" @@ -49,12 +47,10 @@ ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" -HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a" -WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" [targets] -test = ["Aqua", "Combinatorics", "HalfIntegers", "LinearAlgebra", "TensorOperations", "Test", "TestExtras", "WignerSymbols", "ChainRulesCore", "ChainRulesTestUtils", "FiniteDifferences"] +test = ["Aqua", "Combinatorics", "LinearAlgebra", "TensorOperations", "Test", "TestExtras", "ChainRulesCore", "ChainRulesTestUtils", "FiniteDifferences"] diff --git a/TensorKitSectors/Project.toml b/TensorKitSectors/Project.toml deleted file mode 100644 index 7b915d17..00000000 --- a/TensorKitSectors/Project.toml +++ /dev/null @@ -1,21 +0,0 @@ -name = "TensorKitSectors" -uuid = "7ea117e4-0799-4282-ac6e-748b4bc2d0f5" -authors = ["Jutho Haegeman", "Lukas Devos"] -version = "0.1.0" - -[deps] -HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" -WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" - -[extras] -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" -TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a" -TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a" - -[targets] -test = ["Test", "TestExtras", "TestItemRunner", "TensorOperations", "Random", "Reexport"] diff --git a/TensorKitSectors/src/TensorKitSectors.jl b/TensorKitSectors/src/TensorKitSectors.jl deleted file mode 100644 index 3fcb93f2..00000000 --- a/TensorKitSectors/src/TensorKitSectors.jl +++ /dev/null @@ -1,69 +0,0 @@ -# Superselection sectors (quantum numbers): -# for defining graded vector spaces and invariant subspaces of tensor products -#==========================================================================================# - -module TensorKitSectors - -# exports -# ------- -export Sector, Group, AbstractIrrep -export Irrep - -export Nsymbol, Fsymbol, Rsymbol, Asymbol, Bsymbol -export sectorscalartype -export dim, sqrtdim, invsqrtdim, frobeniusschur, twist, fusiontensor, dual -export otimes, deligneproduct, times -export FusionStyle, UniqueFusion, MultipleFusion, SimpleFusion, GenericFusion, - MultiplicityFreeFusion -export BraidingStyle, NoBraiding, SymmetricBraiding, Bosonic, Fermionic, Anyonic -export SectorSet, SectorValues, findindex, vertex_ind2label, vertex_labeltype - -export pentagon_equation, hexagon_equation - -export Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep -export ProductSector -export FermionParity, FermionNumber, FermionSpin -export PlanarTrivial, FibonacciAnyon, IsingAnyon - -# unicode exports -# --------------- -export ⊠, ⊗, × -export ℤ, ℤ₂, ℤ₃, ℤ₄, U₁, SU, SU₂, CU₁ -export fℤ₂, fU₁, fSU₂ - -# imports -# ------- -using Base: SizeUnknown, HasLength, IsInfinite -using Base: HasEltype, EltypeUnknown -using Base.Iterators: product, filter -using Base: tuple_type_head, tuple_type_tail - -using LinearAlgebra: tr -using TensorOperations -using HalfIntegers -using WignerSymbols - -# includes -# -------- -include("auxiliary.jl") -include("sectors.jl") -include("trivial.jl") -include("groups.jl") -include("irreps.jl") # irreps of symmetry groups, with bosonic braiding -include("product.jl") # direct product of different sectors -include("fermions.jl") # irreps with defined fermionparity and fermionic braiding -include("anyons.jl") # non-group sectors - -# precompile -# ---------- -include("precompile.jl") - -function __precompile__() - for I in (Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep, - FermionParity, FermionNumber, FermionSpin, PlanarTrivial, FibonacciAnyon, - IsingAnyon) - precompile_sector(I) - end -end - -end # module TensorKitSectors diff --git a/TensorKitSectors/src/anyons.jl b/TensorKitSectors/src/anyons.jl deleted file mode 100644 index f32adf32..00000000 --- a/TensorKitSectors/src/anyons.jl +++ /dev/null @@ -1,294 +0,0 @@ -# PlanarTrivial -""" - struct PlanarTrivial <: Sector - PlanarTrivial() - -Represents a trivial anyon sector, i.e. a trivial sector without braiding. This is mostly -useful for testing. -""" -struct PlanarTrivial <: Sector end - -Base.IteratorSize(::Type{SectorValues{PlanarTrivial}}) = HasLength() -Base.length(::SectorValues{PlanarTrivial}) = 1 -Base.iterate(::SectorValues{PlanarTrivial}, i=0) = i == 0 ? (PlanarTrivial(), 1) : nothing -function Base.getindex(::SectorValues{PlanarTrivial}, i::Int) - return i == 1 ? PlanarTrivial() : - throw(BoundsError(values(PlanarTrivial), i)) -end -findindex(::SectorValues{PlanarTrivial}, c::PlanarTrivial) = 1 -Base.isless(::PlanarTrivial, ::PlanarTrivial) = false - -Base.one(::Type{PlanarTrivial}) = PlanarTrivial() -Base.conj(::PlanarTrivial) = PlanarTrivial() - -FusionStyle(::Type{PlanarTrivial}) = UniqueFusion() -BraidingStyle(::Type{PlanarTrivial}) = NoBraiding() -Base.isreal(::Type{PlanarTrivial}) = true - -Nsymbol(::Vararg{PlanarTrivial,3}) = 1 -Fsymbol(::Vararg{PlanarTrivial,6}) = 1 - -⊗(::PlanarTrivial, ::PlanarTrivial) = (PlanarTrivial(),) - -# FibonacciAnyons -""" - struct FibonacciAnyon <: Sector - FibonacciAnyon(s::Symbol) - -Represents the anyons of the Fibonacci modular fusion category. It can take two values, -corresponding to the trivial sector `FibonacciAnyon(:I)` and the non-trivial sector -`FibonacciAnyon(:τ)` with fusion rules ``τ ⊗ τ = 1 ⊕ τ``. - -## Fields -- `isone::Bool`: indicates whether the sector corresponds the to trivial anyon `:I` - (`true`), or the non-trivial anyon `:τ` (`false`). -""" -struct FibonacciAnyon <: Sector - isone::Bool - function FibonacciAnyon(s::Symbol) - s in (:I, :τ, :tau) || throw(ArgumentError("Unknown FibonacciAnyon $s.")) - return new(s === :I) - end -end - -Base.IteratorSize(::Type{SectorValues{FibonacciAnyon}}) = HasLength() -Base.length(::SectorValues{FibonacciAnyon}) = 2 -function Base.iterate(::SectorValues{FibonacciAnyon}, i=0) - return i == 0 ? (FibonacciAnyon(:I), 1) : (i == 1 ? (FibonacciAnyon(:τ), 2) : nothing) -end -function Base.getindex(S::SectorValues{FibonacciAnyon}, i) - if i == 1 - return FibonacciAnyon(:I) - elseif i == 2 - return FibonacciAnyon(:τ) - else - throw(BoundsError(S, i)) - end -end -findindex(::SectorValues{FibonacciAnyon}, s::FibonacciAnyon) = 2 - s.isone - -Base.convert(::Type{FibonacciAnyon}, s::Symbol) = FibonacciAnyon(s) -Base.one(::Type{FibonacciAnyon}) = FibonacciAnyon(:I) -Base.conj(s::FibonacciAnyon) = s - -const _goldenratio = Float64(MathConstants.golden) -dim(a::FibonacciAnyon) = isone(a) ? one(_goldenratio) : _goldenratio - -FusionStyle(::Type{FibonacciAnyon}) = SimpleFusion() -BraidingStyle(::Type{FibonacciAnyon}) = Anyonic() -Base.isreal(::Type{FibonacciAnyon}) = false - -⊗(a::FibonacciAnyon, b::FibonacciAnyon) = FibonacciIterator(a, b) - -struct FibonacciIterator - a::FibonacciAnyon - b::FibonacciAnyon -end -Base.IteratorSize(::Type{FibonacciIterator}) = Base.HasLength() -Base.IteratorEltype(::Type{FibonacciIterator}) = Base.HasEltype() -Base.length(iter::FibonacciIterator) = (isone(iter.a) || isone(iter.b)) ? 1 : 2 -Base.eltype(::Type{FibonacciIterator}) = FibonacciAnyon -function Base.iterate(iter::FibonacciIterator, state=1) - I = FibonacciAnyon(:I) - τ = FibonacciAnyon(:τ) - if state == 1 # first iteration - iter.a == I && return (iter.b, 2) - iter.b == I && return (iter.a, 2) - return (I, 2) - elseif state == 2 - (iter.a == iter.b == τ) && return (τ, 3) - return nothing - else - return nothing - end -end - -function Nsymbol(a::FibonacciAnyon, b::FibonacciAnyon, c::FibonacciAnyon) - return isone(a) + isone(b) + isone(c) != 2 -end # zero if one tau and two ones - -function Fsymbol(a::FibonacciAnyon, b::FibonacciAnyon, c::FibonacciAnyon, - d::FibonacciAnyon, e::FibonacciAnyon, f::FibonacciAnyon) - Nsymbol(a, b, e) || return zero(_goldenratio) - Nsymbol(e, c, d) || return zero(_goldenratio) - Nsymbol(b, c, f) || return zero(_goldenratio) - Nsymbol(a, f, d) || return zero(_goldenratio) - - I = FibonacciAnyon(:I) - τ = FibonacciAnyon(:τ) - if a == b == c == d == τ - if e == f == I - return +1 / _goldenratio - elseif e == f == τ - return -1 / _goldenratio - else - return +1 / sqrt(_goldenratio) - end - else - return one(_goldenratio) - end -end - -function Rsymbol(a::FibonacciAnyon, b::FibonacciAnyon, c::FibonacciAnyon) - Nsymbol(a, b, c) || return 0 * cispi(0 / 1) - if isone(a) || isone(b) - return cispi(0 / 1) - else - return isone(c) ? cispi(4 / 5) : cispi(-3 / 5) - end -end - -function Base.show(io::IO, a::FibonacciAnyon) - s = isone(a) ? ":I" : ":τ" - return get(io, :typeinfo, nothing) === FibonacciAnyon ? - print(io, s) : print(io, "FibonacciAnyon(", s, ")") -end - -Base.hash(a::FibonacciAnyon, h::UInt) = hash(a.isone, h) -Base.isless(a::FibonacciAnyon, b::FibonacciAnyon) = isless(!a.isone, !b.isone) - -# IsingAnyons -""" - struct IsingAnyon <: Sector - IsingAnyon(s::Symbol) - -Represents the anyons of the Ising modular fusion category. It can take three values, -corresponding to the trivial sector `IsingAnyon(:I)` and the non-trivial sectors -`IsingAnyon(:σ)` and `IsingAnyon(:ψ)`, with fusion rules ``ψ ⊗ ψ = 1``, ``σ ⊗ ψ = σ``, and -``σ ⊗ σ = 1 ⊕ ψ``. - -## Fields -- `s::Symbol`: the label of the represented anyon, which can be `:I`, `:σ`, or `:ψ`. -""" -struct IsingAnyon <: Sector - s::Symbol - function IsingAnyon(s::Symbol) - s == :sigma && (s = :σ) - s == :psi && (s = :ψ) - if !(s in (:I, :σ, :ψ)) - throw(ValueError("Unknown IsingAnyon $s.")) - end - return new(s) - end -end - -const all_isinganyons = (IsingAnyon(:I), IsingAnyon(:σ), IsingAnyon(:ψ)) - -Base.IteratorSize(::Type{SectorValues{IsingAnyon}}) = HasLength() -Base.length(::SectorValues{IsingAnyon}) = length(all_isinganyons) -Base.iterate(::SectorValues{IsingAnyon}, i=1) = iterate(all_isinganyons, i) -Base.getindex(S::SectorValues{IsingAnyon}, i) = getindex(all_isinganyons, i) - -function findindex(::SectorValues{IsingAnyon}, a::IsingAnyon) - a == all_isinganyons[1] && return 1 - a == all_isinganyons[2] && return 2 - return 3 -end - -Base.convert(::Type{IsingAnyon}, s::Symbol) = IsingAnyon(s) -Base.one(::Type{IsingAnyon}) = IsingAnyon(:I) -Base.conj(s::IsingAnyon) = s - -dim(a::IsingAnyon) = a.s == :σ ? sqrt(2) : 1.0 - -FusionStyle(::Type{IsingAnyon}) = SimpleFusion() -BraidingStyle(::Type{IsingAnyon}) = Anyonic() -Base.isreal(::Type{IsingAnyon}) = false - -⊗(a::IsingAnyon, b::IsingAnyon) = IsingIterator(a, b) - -struct IsingIterator - a::IsingAnyon - b::IsingAnyon -end - -Base.IteratorSize(::Type{IsingIterator}) = Base.HasLength() -Base.IteratorEltype(::Type{IsingIterator}) = Base.HasEltype() -Base.eltype(::Type{IsingIterator}) = IsingAnyon - -function Base.length(iter::IsingIterator) - σ = IsingAnyon(:σ) - return (iter.a == σ && iter.b == σ) ? 2 : 1 -end - -function Base.iterate(iter::IsingIterator, state=1) - I, σ, ψ = all_isinganyons - if state == 1 # first iteration - iter.a == I && return (iter.b, 2) - iter.b == I && return (iter.a, 2) - iter.a == σ && iter.b == ψ && return (σ, 2) - iter.a == ψ && iter.b == σ && return (σ, 2) - return (I, 2) - elseif state == 2 - (iter.a == iter.b == σ) && return (ψ, 3) - return nothing - else - return nothing - end -end - -function Nsymbol(a::IsingAnyon, b::IsingAnyon, c::IsingAnyon) - I, σ, ψ = all_isinganyons - return ((a == I && b == c) - || (b == I && a == c) - || (c == I && a == b) - || (a == σ && b == σ && c == ψ) - || (a == σ && b == ψ && c == σ) - || (a == ψ && b == σ && c == σ)) -end - -function Fsymbol(a::IsingAnyon, b::IsingAnyon, c::IsingAnyon, - d::IsingAnyon, e::IsingAnyon, f::IsingAnyon) - Nsymbol(a, b, e) || return 0.0 - Nsymbol(e, c, d) || return 0.0 - Nsymbol(b, c, f) || return 0.0 - Nsymbol(a, f, d) || return 0.0 - I, σ, ψ = all_isinganyons - if a == b == c == d == σ - if e == f == ψ - return -1.0 / sqrt(2.0) - else - return 1.0 / sqrt(2.0) - end - end - if e == f == σ - if a == c == σ && b == d == ψ - return -1.0 - elseif a == c == ψ && b == d == σ - return -1.0 - end - end - return 1.0 -end - -function Rsymbol(a::IsingAnyon, b::IsingAnyon, c::IsingAnyon) - Nsymbol(a, b, c) || return complex(0.0) - I, σ, ψ = all_isinganyons - if c == I - if b == a == σ - return cispi(-1 / 8) - elseif b == a == ψ - return complex(-1.0) - end - elseif c == σ && (a == σ && b == ψ || a == ψ && b == σ) - return -1.0im - elseif c == ψ && a == b == σ - return cispi(3 / 8) - end - return complex(1.0) -end - -function Base.show(io::IO, a::IsingAnyon) - if get(io, :typeinfo, nothing) === IsingAnyon - return print(io, ":$(a.s)") - else - return print(io, "IsingAnyon(:$(a.s))") - end -end - -Base.hash(s::IsingAnyon, h::UInt) = hash(s.s, h) - -function Base.isless(a::IsingAnyon, b::IsingAnyon) - vals = SectorValues{IsingAnyon}() - return isless(findindex(vals, a), findindex(vals, b)) -end diff --git a/TensorKitSectors/src/auxiliary.jl b/TensorKitSectors/src/auxiliary.jl deleted file mode 100644 index 3ea734cd..00000000 --- a/TensorKitSectors/src/auxiliary.jl +++ /dev/null @@ -1,90 +0,0 @@ -# kronecker product with arbitrary number of dimensions -_kron(A, B, C, D...) = _kron(_kron(A, B), C, D...) -function _kron(A, B) - sA = size(A) - sB = size(B) - s = map(*, sA, sB) - C = similar(A, promote_type(eltype(A), eltype(B)), s) - for IA in eachindex(IndexCartesian(), A) - for IB in eachindex(IndexCartesian(), B) - I = CartesianIndex(IB.I .+ (IA.I .- 1) .* sB) - C[I] = A[IA] * B[IB] - end - end - return C -end - -# Manhattan based distance enumeration: I is supposed to be one-based index -# TODO: is there any way to make this faster? - -# forward mapping from multidimensional to single Manhattan index -@inline function num_manhattan_points(d::Int, sz::Dims{1}) - return Int(sz[1] > d) -end -@inline function num_manhattan_points(d::Int, sz::Dims{N}) where {N} - d == 0 && return 1 - num = 0 - for i in 1:min(sz[1], d + 1) - num += num_manhattan_points(d - i + 1, Base.tail(sz)) - end - return num -end - -@inline localoffset(d, I::Tuple{Int}, sz::Tuple{Int}) = 0 -@inline function localoffset(d, I, sz) - offset = 0 - for i in 1:(I[1] - 1) - offset += num_manhattan_points(d - i + 1, Base.tail(sz)) - end - offset += localoffset(d - I[1] + 1, Base.tail(I), Base.tail(sz)) - return offset -end - -to_manhattan_index(::Tuple{}, ::Tuple{}) = 1 -to_manhattan_index(I::Tuple{Int}, ::Tuple{Int}) = I[1] -function to_manhattan_index(I::Dims{N}, sz::Dims{N}) where {N} - d = sum(I) - N # Manhattan distance to origin for one-based indices - d == 0 && return 1 - index = 1 - # count all the points with smaller Manhatten distance - for k in 0:(d - 1) - index += num_manhattan_points(k, sz) - end - return index + localoffset(d, I, sz) -end - -# inverse mapping -@inline invertlocaloffset(d, offset, sz::Tuple{Int}) = (d + 1,) -@inline function invertlocaloffset(d, offset, sz) - i₁ = 1 - while i₁ < sz[1] - jump = num_manhattan_points(d - i₁ + 1, Base.tail(sz)) - if offset < jump - break - end - offset -= jump - i₁ += 1 - end - return (i₁, invertlocaloffset(d - i₁ + 1, offset, Base.tail(sz))...) -end - -manhattan_to_multidimensional_index(index::Int, ::Dims{0}) = () -manhattan_to_multidimensional_index(index::Int, ::Dims{1}) = (index,) -function manhattan_to_multidimensional_index(index::Int, sz::Dims{N}) where {N} - # find the layer of the Manhattan distance - index == 1 && return ntuple(one, Val(N)) - offset = 1 - d = 1 - while true - currentlayer = num_manhattan_points(d, sz) - if index <= offset + currentlayer - break - end - d += 1 - offset += currentlayer - end - - # find the position within the layer - index -= offset - return invertlocaloffset(d, index - 1, sz) -end diff --git a/TensorKitSectors/src/fermions.jl b/TensorKitSectors/src/fermions.jl deleted file mode 100644 index a1d4b522..00000000 --- a/TensorKitSectors/src/fermions.jl +++ /dev/null @@ -1,105 +0,0 @@ -""" - FermionParity <: Sector - -Represents sectors with fermion parity. The fermion parity is a ℤ₂ quantum number that -yields an additional sign when two odd fermions are exchanged. - -## Fields -- `isodd::Bool`: indicates whether the fermion parity is odd (`true`) or even (`false`). - -See also: [`FermionNumber`](@ref), [`FermionSpin`](@ref) -""" -struct FermionParity <: Sector - isodd::Bool -end -const fℤ₂ = FermionParity -fermionparity(f::FermionParity) = f.isodd - -Base.convert(::Type{FermionParity}, a::FermionParity) = a -Base.convert(::Type{FermionParity}, a) = FermionParity(a) - -Base.IteratorSize(::Type{SectorValues{FermionParity}}) = HasLength() -Base.length(::SectorValues{FermionParity}) = 2 -function Base.iterate(::SectorValues{FermionParity}, i=0) - return i == 2 ? nothing : (FermionParity(i), i + 1) -end -function Base.getindex(::SectorValues{FermionParity}, i) - return 1 <= i <= 2 ? FermionParity(i - 1) : throw(BoundsError(values(FermionParity), i)) -end -findindex(::SectorValues{FermionParity}, f::FermionParity) = f.isodd ? 2 : 1 - -Base.one(::Type{FermionParity}) = FermionParity(false) -Base.conj(f::FermionParity) = f -dim(f::FermionParity) = 1 - -FusionStyle(::Type{FermionParity}) = UniqueFusion() -BraidingStyle(::Type{FermionParity}) = Fermionic() -Base.isreal(::Type{FermionParity}) = true - -⊗(a::FermionParity, b::FermionParity) = (FermionParity(a.isodd ⊻ b.isodd),) - -function Nsymbol(a::FermionParity, b::FermionParity, c::FermionParity) - return (a.isodd ⊻ b.isodd) == c.isodd -end -function Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I<:FermionParity} - return Int(Nsymbol(a, b, e) * Nsymbol(e, c, d) * Nsymbol(b, c, f) * Nsymbol(a, f, d)) -end -function Rsymbol(a::I, b::I, c::I) where {I<:FermionParity} - return a.isodd && b.isodd ? -Int(Nsymbol(a, b, c)) : Int(Nsymbol(a, b, c)) -end -twist(a::FermionParity) = a.isodd ? -1 : +1 - -function fusiontensor(a::I, b::I, c::I) where {I<:FermionParity} - @warn "FermionParity Arrays do not preserve categorical properties." maxlog = 1 - return fill(Int(Nsymbol(a, b, c)), (1, 1, 1, 1)) -end - -function Base.show(io::IO, a::FermionParity) - if get(io, :typeinfo, nothing) === typeof(a) - print(io, Int(a.isodd)) - else - print(io, type_repr(typeof(a)), "(", Int(a.isodd), ")") - end -end -type_repr(::Type{FermionParity}) = "FermionParity" - -Base.hash(f::FermionParity, h::UInt) = hash(f.isodd, h) -Base.isless(a::FermionParity, b::FermionParity) = isless(a.isodd, b.isodd) - -# Common fermionic combinations -# ----------------------------- - -""" - const FermionNumber = U1Irrep ⊠ FermionParity - FermionNumber(a::Int) - -Represents the fermion number as the direct product of a ``U₁`` irrep `a` and a fermion -parity, with the restriction that the fermion parity is odd if and only if `a` is odd. - -See also: [`U1Irrep`](@ref), [`FermionParity`](@ref) -""" -const FermionNumber = U1Irrep ⊠ FermionParity -const fU₁ = FermionNumber -FermionNumber(a::Int) = U1Irrep(a) ⊠ FermionParity(isodd(a)) -type_repr(::Type{FermionNumber}) = "FermionNumber" - -# convenience default converter -> allows Vect[FermionNumber](1 => 1) -Base.convert(::Type{FermionNumber}, a::Int) = FermionNumber(a) - -""" - const FermionSpin = SU2Irrep ⊠ FermionParity - FermionSpin(j::Real) - -Represents the fermion spin as the direct product of a ``SU₂`` irrep `j` and a fermion -parity, with the restriction that the fermion parity is odd if `2 * j` is odd. - -See also: [`SU2Irrep`](@ref), [`FermionParity`](@ref) -""" -const FermionSpin = SU2Irrep ⊠ FermionParity -const fSU₂ = FermionSpin -FermionSpin(j::Real) = (s = SU2Irrep(j); - s ⊠ FermionParity(isodd(twice(s.j)))) -type_repr(::Type{FermionSpin}) = "FermionSpin" - -# convenience default converter -> allows Vect[FermionSpin](1 => 1) -Base.convert(::Type{FermionSpin}, a::Real) = FermionSpin(a) diff --git a/TensorKitSectors/src/groups.jl b/TensorKitSectors/src/groups.jl deleted file mode 100644 index 2d34e069..00000000 --- a/TensorKitSectors/src/groups.jl +++ /dev/null @@ -1,68 +0,0 @@ -# Groups -#------------------------------------------------------------------------------# -abstract type Group end -abstract type AbelianGroup <: Group end - -abstract type ℤ{N} <: AbelianGroup end -abstract type U₁ <: AbelianGroup end -abstract type SU{N} <: Group end -abstract type CU₁ <: Group end - -const ℤ₂ = ℤ{2} -const ℤ₃ = ℤ{3} -const ℤ₄ = ℤ{4} -const SU₂ = SU{2} - -type_repr(::Type{ℤ₂}) = "ℤ₂" -type_repr(::Type{ℤ₃}) = "ℤ₃" -type_repr(::Type{ℤ₄}) = "ℤ₄" -type_repr(::Type{SU₂}) = "SU₂" -type_repr(T::Type) = repr(T) - -const GroupTuple = Tuple{Vararg{Group}} - -abstract type ProductGroup{T<:GroupTuple} <: Group end - -""" - ×(G::Vararg{Type{<:Group}}) -> ProductGroup{Tuple{G...}} - times(G::Vararg{Type{<:Group}}) -> ProductGroup{Tuple{G...}} - -Construct the direct product of a (list of) groups. -""" -function ×(::Vararg{Type{<:Group}}) end -const times = × - -×(a::Type{<:Group}, b::Type{<:Group}, c::Type{<:Group}...) = ×(×(a, b), c...) -×(G::Type{<:Group}) = ProductGroup{Tuple{G}} -×(G1::Type{ProductGroup{Tuple{}}}, -G2::Type{ProductGroup{T}}) where {T<:GroupTuple} = G2 -function ×(G1::Type{ProductGroup{T1}}, - G2::Type{ProductGroup{T2}}) where {T1<:GroupTuple,T2<:GroupTuple} - return tuple_type_head(T1) × (ProductGroup{tuple_type_tail(T1)} × G2) -end -×(G1::Type{ProductGroup{Tuple{}}}, G2::Type{<:Group}) = ProductGroup{Tuple{G2}} -function ×(G1::Type{ProductGroup{T}}, G2::Type{<:Group}) where {T<:GroupTuple} - return Base.tuple_type_head(T) × (ProductGroup{Base.tuple_type_tail(T)} × G2) -end -function ×(G1::Type{<:Group}, G2::Type{ProductGroup{T}}) where {T<:GroupTuple} - return ProductGroup{Base.tuple_type_cons(G1, T)} -end -×(G1::Type{<:Group}, G2::Type{<:Group}) = ProductGroup{Tuple{G1,G2}} - -function type_repr(G::Type{<:ProductGroup}) - T = G.parameters[1] - groups = T.parameters - if length(groups) == 1 - s = "ProductGroup{Tuple{" * type_repr(groups[1]) * "}}" - else - s = "(" - for i in 1:length(groups) - if i != 1 - s *= " × " - end - s *= type_repr(groups[i]) - end - s *= ")" - end - return s -end diff --git a/TensorKitSectors/src/irreps.jl b/TensorKitSectors/src/irreps.jl deleted file mode 100644 index e0096a81..00000000 --- a/TensorKitSectors/src/irreps.jl +++ /dev/null @@ -1,483 +0,0 @@ -# Sectors corresponding to irreducible representations of compact groups -#==============================================================================# -# Irreps of groups -#------------------------------------------------------------------------------# -""" - abstract type AbstractIrrep{G<:Group} <: Sector end - -Abstract supertype for sectors which corresponds to irreps (irreducible representations) of -a group `G`. As we assume unitary representations, these would be finite groups or compact -Lie groups. Note that this could also include projective rather than linear representations. - -Actual concrete implementations of those irreps can be obtained as `Irrep[G]`, or via their -actual name, which generically takes the form `(asciiG)Irrep`, i.e. the ASCII spelling of -the group name followed by `Irrep`. - -All irreps have [`BraidingStyle`](@ref) equal to `Bosonic()` and thus trivial twists. -""" -abstract type AbstractIrrep{G<:Group} <: Sector end # irreps have integer quantum dimensions -BraidingStyle(::Type{<:AbstractIrrep}) = Bosonic() - -struct IrrepTable end -""" - const Irrep - -A constant of a singleton type used as `Irrep[G]` with `G<:Group` a type of group, to -construct or obtain a concrete subtype of `AbstractIrrep{G}` that implements the data -structure used to represent irreducible representations of the group `G`. -""" -const Irrep = IrrepTable() - -type_repr(::Type{<:AbstractIrrep{G}}) where {G<:Group} = "Irrep[" * type_repr(G) * "]" -function Base.show(io::IO, c::AbstractIrrep) - I = typeof(c) - if get(io, :typeinfo, nothing) !== I - print(io, type_repr(I), "(") - for k in 1:fieldcount(I) - k > 1 && print(io, ", ") - print(io, getfield(c, k)) - end - print(io, ")") - else - fieldcount(I) > 1 && print(io, "(") - for k in 1:fieldcount(I) - k > 1 && print(io, ", ") - print(io, getfield(c, k)) - end - fieldcount(I) > 1 && print(io, ")") - end -end - -const AbelianIrrep{G} = AbstractIrrep{G} where {G<:AbelianGroup} -FusionStyle(::Type{<:AbelianIrrep}) = UniqueFusion() -Base.isreal(::Type{<:AbelianIrrep}) = true - -Nsymbol(a::I, b::I, c::I) where {I<:AbelianIrrep} = c == first(a ⊗ b) -function Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I<:AbelianIrrep} - return Int(Nsymbol(a, b, e) * Nsymbol(e, c, d) * Nsymbol(b, c, f) * Nsymbol(a, f, d)) -end -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)) - -function fusiontensor(a::I, b::I, c::I) where {I<:AbelianIrrep} - return fill(Int(Nsymbol(a, b, c)), (1, 1, 1, 1)) -end - -# ZNIrrep: irreps of Z_N are labelled by integers mod N; do we ever want N > 64? -""" - struct ZNIrrep{N} <: AbstractIrrep{ℤ{N}} - ZNIrrep{N}(n::Integer) - Irrep[ℤ{N}](n::Integer) - -Represents irreps of the group ``ℤ_N`` for some value of `N<64`. (We need 2*(N-1) <= 127 in -order for a ⊗ b to work correctly.) For `N` equals `2`, `3` or `4`, `ℤ{N}` can be replaced -by `ℤ₂`, `ℤ₃`, `ℤ₄`. An arbitrary `Integer` `n` can be provided to the constructor, but only -the value `mod(n, N)` is relevant. - -## Fields -- `n::Int8`: the integer label of the irrep, modulo `N`. -""" -struct ZNIrrep{N} <: AbstractIrrep{ℤ{N}} - n::Int8 - function ZNIrrep{N}(n::Integer) where {N} - @assert N < 64 - return new{N}(mod(n, N)) - end -end -Base.getindex(::IrrepTable, ::Type{ℤ{N}}) where {N} = ZNIrrep{N} -Base.convert(Z::Type{<:ZNIrrep}, n::Real) = Z(n) -const Z2Irrep = ZNIrrep{2} -const Z3Irrep = ZNIrrep{3} -const Z4Irrep = ZNIrrep{4} - -Base.one(::Type{ZNIrrep{N}}) where {N} = ZNIrrep{N}(0) -Base.conj(c::ZNIrrep{N}) where {N} = ZNIrrep{N}(-c.n) -⊗(c1::ZNIrrep{N}, c2::ZNIrrep{N}) where {N} = (ZNIrrep{N}(c1.n + c2.n),) - -Base.IteratorSize(::Type{SectorValues{ZNIrrep{N}}}) where {N} = HasLength() -Base.length(::SectorValues{ZNIrrep{N}}) where {N} = N -function Base.iterate(::SectorValues{ZNIrrep{N}}, i=0) where {N} - return i == N ? nothing : (ZNIrrep{N}(i), i + 1) -end -function Base.getindex(::SectorValues{ZNIrrep{N}}, i::Int) where {N} - return 1 <= i <= N ? ZNIrrep{N}(i - 1) : throw(BoundsError(values(ZNIrrep{N}), i)) -end -findindex(::SectorValues{ZNIrrep{N}}, c::ZNIrrep{N}) where {N} = c.n + 1 - -Base.hash(c::ZNIrrep{N}, h::UInt) where {N} = hash(c.n, h) -Base.isless(c1::ZNIrrep{N}, c2::ZNIrrep{N}) where {N} = isless(c1.n, c2.n) - -# U1Irrep: irreps of U1 are labelled by integers -""" - struct U1Irrep <: AbstractIrrep{U₁} - U1Irrep(charge::Real) - Irrep[U₁](charge::Real) - -Represents irreps of the group ``U₁``. The irrep is labelled by a charge, which should be -an integer for a linear representation. However, it is often useful to allow half integers -to represent irreps of ``U₁`` subgroups of ``SU₂``, such as the Sz of spin-1/2 system. -Hence, the charge is stored as a `HalfInt` from the package HalfIntegers.jl, but can be -entered as arbitrary `Real`. The sequence of the charges is: 0, 1/2, -1/2, 1, -1, ... - -## Fields -- `charge::HalfInt`: the label of the irrep, which can be any half integer. -""" -struct U1Irrep <: AbstractIrrep{U₁} - charge::HalfInt -end -Base.getindex(::IrrepTable, ::Type{U₁}) = U1Irrep -Base.convert(::Type{U1Irrep}, c::Real) = U1Irrep(c) - -Base.one(::Type{U1Irrep}) = U1Irrep(0) -Base.conj(c::U1Irrep) = U1Irrep(-c.charge) -⊗(c1::U1Irrep, c2::U1Irrep) = (U1Irrep(c1.charge + c2.charge),) - -Base.IteratorSize(::Type{SectorValues{U1Irrep}}) = IsInfinite() -function Base.iterate(::SectorValues{U1Irrep}, i=0) - return i <= 0 ? (U1Irrep(half(i)), (-i + 1)) : (U1Irrep(half(i)), -i) -end -function Base.getindex(::SectorValues{U1Irrep}, i::Int) - i < 1 && throw(BoundsError(values(U1Irrep), i)) - return U1Irrep(iseven(i) ? half(i >> 1) : -half(i >> 1)) -end -function findindex(::SectorValues{U1Irrep}, c::U1Irrep) - return (n = twice(c.charge); 2 * abs(n) + (n <= 0)) -end - -Base.hash(c::U1Irrep, h::UInt) = hash(c.charge, h) -@inline function Base.isless(c1::U1Irrep, c2::U1Irrep) - return isless(abs(c1.charge), abs(c2.charge)) || zero(HalfInt) < c1.charge == -c2.charge -end - -# Non-abelian groups -#------------------------------------------------------------------------------# -# SU2Irrep: irreps of SU2 are labelled by half integers j -struct SU2IrrepException <: Exception end -function Base.show(io::IO, ::SU2IrrepException) - return print(io, - "Irreps of (bosonic or fermionic) `SU₂` should be labelled by non-negative half integers, i.e. elements of `Rational{Int}` with denominator 1 or 2") -end - -""" - struct SU2Irrep <: AbstractIrrep{SU₂} - SU2Irrep(j::Real) - Irrep[SU₂](j::Real) - -Represents irreps of the group ``SU₂``. The irrep is labelled by a half integer `j` which -can be entered as an abitrary `Real`, but is stored as a `HalfInt` from the HalfIntegers.jl -package. - -## Fields -- `j::HalfInt`: the label of the irrep, which can be any non-negative half integer. -""" -struct SU2Irrep <: AbstractIrrep{SU₂} - j::HalfInt - function SU2Irrep(j) - j >= zero(j) || error("Not a valid SU₂ irrep") - return new(j) - end -end -Base.getindex(::IrrepTable, ::Type{SU₂}) = SU2Irrep -Base.convert(::Type{SU2Irrep}, j::Real) = SU2Irrep(j) - -const _su2one = SU2Irrep(zero(HalfInt)) -Base.one(::Type{SU2Irrep}) = _su2one -Base.conj(s::SU2Irrep) = s -⊗(s1::SU2Irrep, s2::SU2Irrep) = SectorSet{SU2Irrep}(abs(s1.j - s2.j):(s1.j + s2.j)) - -Base.IteratorSize(::Type{SectorValues{SU2Irrep}}) = IsInfinite() -Base.iterate(::SectorValues{SU2Irrep}, i=0) = (SU2Irrep(half(i)), i + 1) -function Base.getindex(::SectorValues{SU2Irrep}, i::Int) - return 1 <= i ? SU2Irrep(half(i - 1)) : throw(BoundsError(values(SU2Irrep), i)) -end -findindex(::SectorValues{SU2Irrep}, s::SU2Irrep) = twice(s.j) + 1 - -dim(s::SU2Irrep) = twice(s.j) + 1 - -FusionStyle(::Type{SU2Irrep}) = SimpleFusion() -sectorscalartype(::Type{SU2Irrep}) = Float64 -Base.isreal(::Type{SU2Irrep}) = true - -Nsymbol(sa::SU2Irrep, sb::SU2Irrep, sc::SU2Irrep) = WignerSymbols.δ(sa.j, sb.j, sc.j) - -function Fsymbol(s1::SU2Irrep, s2::SU2Irrep, s3::SU2Irrep, - s4::SU2Irrep, s5::SU2Irrep, s6::SU2Irrep) - if all(==(_su2one), (s1, s2, s3, s4, s5, s6)) - return 1.0 - else - return sqrtdim(s5) * sqrtdim(s6) * - WignerSymbols.racahW(sectorscalartype(SU2Irrep), s1.j, s2.j, - s4.j, s3.j, s5.j, s6.j) - end -end - -function Rsymbol(sa::SU2Irrep, sb::SU2Irrep, sc::SU2Irrep) - Nsymbol(sa, sb, sc) || return zero(sectorscalartype(SU2Irrep)) - return iseven(convert(Int, sa.j + sb.j - sc.j)) ? one(sectorscalartype(SU2Irrep)) : - -one(sectorscalartype(SU2Irrep)) -end - -function fusiontensor(a::SU2Irrep, b::SU2Irrep, c::SU2Irrep) - C = Array{Float64}(undef, dim(a), dim(b), dim(c), 1) - ja, jb, jc = a.j, b.j, c.j - - for kc in 1:dim(c), kb in 1:dim(b), ka in 1:dim(a) - C[ka, kb, kc, 1] = WignerSymbols.clebschgordan(ja, ja + 1 - ka, jb, jb + 1 - kb, jc, - jc + 1 - kc) - end - return C -end - -Base.hash(s::SU2Irrep, h::UInt) = hash(s.j, h) -Base.isless(s1::SU2Irrep, s2::SU2Irrep) = isless(s1.j, s2.j) - -# U₁ ⋊ C (U₁ and charge conjugation) -""" - struct CU1Irrep <: AbstractIrrep{CU₁} - CU1Irrep(j, s = ifelse(j>zero(j), 2, 0)) - Irrep[CU₁](j, s = ifelse(j>zero(j), 2, 0)) - -Represents irreps of the group ``U₁ ⋊ C`` (``U₁`` and charge conjugation or reflection), -which is also known as just `O₂`. - -## Fields -- `j::HalfInt`: the value of the ``U₁`` charge. -- `s::Int`: the representation of charge conjugation. - -They can take values: -* if `j == 0`, `s = 0` (trivial charge conjugation) or - `s = 1` (non-trivial charge conjugation) -* if `j > 0`, `s = 2` (two-dimensional representation) -""" -struct CU1Irrep <: AbstractIrrep{CU₁} - j::HalfInt # value of the U1 charge - s::Int # rep of charge conjugation: - # if j == 0, s = 0 (trivial) or s = 1 (non-trivial), - # else s = 2 (two-dimensional representation) - # Let constructor take the actual half integer value j - function CU1Irrep(j::Real, s::Integer=ifelse(j > zero(j), 2, 0)) - if ((j > zero(j) && s == 2) || (j == zero(j) && (s == 0 || s == 1))) - new(j, s) - else - error("Not a valid CU₁ irrep") - end - end -end -Base.getindex(::IrrepTable, ::Type{CU₁}) = CU1Irrep -Base.convert(::Type{CU1Irrep}, (j, s)::Tuple{Real,Integer}) = CU1Irrep(j, s) - -Base.IteratorSize(::Type{SectorValues{CU1Irrep}}) = IsInfinite() -function Base.iterate(::SectorValues{CU1Irrep}, state=(0, 0)) - j, s = state - if iszero(j) && s == 0 - return CU1Irrep(j, s), (j, 1) - elseif iszero(j) && s == 1 - return CU1Irrep(j, s), (j + 1, 2) - else - return CU1Irrep(half(j), s), (j + 1, 2) - end -end -function Base.getindex(::SectorValues{CU1Irrep}, i::Int) - i < 1 && throw(BoundsError(values(CU1Irrep), i)) - if i == 1 - return CU1Irrep(0, 0) - elseif i == 2 - return CU1Irrep(0, 1) - else - return CU1Irrep(half(i - 2), 2) - end -end -findindex(::SectorValues{CU1Irrep}, c::CU1Irrep) = twice(c.j) + iszero(c.j) + c.s - -Base.hash(c::CU1Irrep, h::UInt) = hash(c.s, hash(c.j, h)) -function Base.isless(c1::CU1Irrep, c2::CU1Irrep) - return isless(c1.j, c2.j) || (c1.j == c2.j == zero(HalfInt) && c1.s < c2.s) -end - -# CU1Irrep(j::Real, s::Int = ifelse(j>0, 2, 0)) = CU1Irrep(convert(HalfInteger, j), s) - -Base.convert(::Type{CU1Irrep}, j::Real) = CU1Irrep(j) -Base.convert(::Type{CU1Irrep}, js::Tuple{Real,Int}) = CU1Irrep(js...) - -Base.one(::Type{CU1Irrep}) = CU1Irrep(zero(HalfInt), 0) -Base.conj(c::CU1Irrep) = c - -struct CU1ProdIterator - a::CU1Irrep - b::CU1Irrep -end -function Base.iterate(p::CU1ProdIterator, s::Int=1) - if s == 1 - if p.a.j == p.b.j == zero(HalfInt) - return CU1Irrep(zero(HalfInt), xor(p.a.s, p.b.s)), 4 - elseif p.a.j == zero(HalfInt) - return p.b, 4 - elseif p.b.j == zero(HalfInt) - return p.a, 4 - elseif p.a == p.b # != zero - return one(CU1Irrep), 2 - else - return CU1Irrep(abs(p.a.j - p.b.j)), 3 - end - elseif s == 2 - return CU1Irrep(zero(HalfInt), 1), 3 - elseif s == 3 - CU1Irrep(p.a.j + p.b.j), 4 - else - return nothing - end -end -function Base.length(p::CU1ProdIterator) - if p.a.j == zero(HalfInt) || p.b.j == zero(HalfInt) - return 1 - elseif p.a == p.b - return 3 - else - return 2 - end -end -Base.eltype(::Type{CU1ProdIterator}) = CU1Irrep - -⊗(a::CU1Irrep, b::CU1Irrep) = CU1ProdIterator(a, b) - -dim(c::CU1Irrep) = ifelse(c.j == zero(HalfInt), 1, 2) - -FusionStyle(::Type{CU1Irrep}) = SimpleFusion() -sectorscalartype(::Type{CU1Irrep}) = Float64 -Base.isreal(::Type{CU1Irrep}) = true - -function Nsymbol(a::CU1Irrep, b::CU1Irrep, c::CU1Irrep) - return ifelse(c.s == 0, (a.j == b.j) & ((a.s == b.s == 2) | (a.s == b.s)), - ifelse(c.s == 1, (a.j == b.j) & ((a.s == b.s == 2) | (a.s != b.s)), - (c.j == a.j + b.j) | (c.j == abs(a.j - b.j)))) -end - -function Fsymbol(a::CU1Irrep, b::CU1Irrep, c::CU1Irrep, - d::CU1Irrep, e::CU1Irrep, f::CU1Irrep) - Nabe = convert(Int, Nsymbol(a, b, e)) - Necd = convert(Int, Nsymbol(e, c, d)) - Nbcf = convert(Int, Nsymbol(b, c, f)) - Nafd = convert(Int, Nsymbol(a, f, d)) - - T = sectorscalartype(CU1Irrep) - Nabe * Necd * Nbcf * Nafd == 0 && return zero(T) - - op = CU1Irrep(0, 0) - om = CU1Irrep(0, 1) - - if a == op || b == op || c == op - return one(T) - end - if (a == b == om) || (a == c == om) || (b == c == om) - return one(T) - end - if a == om - if d.j == zero(HalfInt) - return one(T) - else - return (d.j == c.j - b.j) ? -one(T) : one(T) - end - end - if b == om - return (d.j == abs(a.j - c.j)) ? -one(T) : one(T) - end - if c == om - return (d.j == a.j - b.j) ? -one(T) : one(T) - end - # from here on, a, b, c are neither 0+ or 0- - s = T(sqrt(2) / 2) - if a == b == c - if d == a - if e.j == 0 - if f.j == 0 - return f.s == 1 ? T(-1 // 2) : T(1 // 2) - else - return e.s == 1 ? -s : s - end - else - return f.j == 0 ? s : zero(T) - end - else - return one(T) - end - end - if a == b # != c - if d == c - if f.j == b.j + c.j - return e.s == 1 ? -s : s - else - return s - end - else - return one(T) - end - end - if b == c - if d == a - if e.j == a.j + b.j - return s - else - return f.s == 1 ? -s : s - end - else - return one(T) - end - end - if a == c - if d == b - if e.j == f.j - return zero(T) - else - return one(T) - end - else - return d.s == 1 ? -one(T) : one(T) - end - end - if d == om - return b.j == a.j + c.j ? -one(T) : one(T) - end - return one(T) -end - -function Rsymbol(a::CU1Irrep, b::CU1Irrep, c::CU1Irrep) - R = convert(sectorscalartype(CU1Irrep), Nsymbol(a, b, c)) - return c.s == 1 && a.j > 0 ? -R : R -end - -function fusiontensor(a::CU1Irrep, b::CU1Irrep, c::CU1Irrep) - C = fill(0.0, dim(a), dim(b), dim(c), 1) - !Nsymbol(a, b, c) && return C - if c.j == 0 - if a.j == b.j == 0 - C[1, 1, 1, 1] = 1.0 - else - if c.s == 0 - C[1, 2, 1, 1] = 1.0 / sqrt(2) - C[2, 1, 1, 1] = 1.0 / sqrt(2) - else - C[1, 2, 1, 1] = 1.0 / sqrt(2) - C[2, 1, 1, 1] = -1.0 / sqrt(2) - end - end - elseif a.j == 0 - C[1, 1, 1, 1] = 1.0 - C[1, 2, 2, 1] = a.s == 1 ? -1.0 : 1.0 - elseif b.j == 0 - C[1, 1, 1, 1] = 1.0 - C[2, 1, 2, 1] = b.s == 1 ? -1.0 : 1.0 - elseif c.j == a.j + b.j - C[1, 1, 1, 1] = 1.0 - C[2, 2, 2, 1] = 1.0 - elseif c.j == a.j - b.j - C[1, 2, 1, 1] = 1.0 - C[2, 1, 2, 1] = 1.0 - elseif c.j == b.j - a.j - C[2, 1, 1, 1] = 1.0 - C[1, 2, 2, 1] = 1.0 - end - return C -end -frobeniusschur(::CU1Irrep) = 1 diff --git a/TensorKitSectors/src/precompile.jl b/TensorKitSectors/src/precompile.jl deleted file mode 100644 index 9e216017..00000000 --- a/TensorKitSectors/src/precompile.jl +++ /dev/null @@ -1,32 +0,0 @@ -""" - precompile_sector(I::Type{<:Sector}) - -Precompile common methods for the given sector type. -""" -function precompile_sector(::Type{I}) where {I<:Sector} - precompile(Nsymbol, (I, I, I)) - precompile(Fsymbol, (I, I, I, I, I, I)) - precompile(Rsymbol, (I, I, I)) - precompile(Asymbol, (I, I, I)) - precompile(Bsymbol, (I, I, I)) - - precompile(⊗, (I,)) - precompile(⊗, (I, I)) - precompile(⊗, (I, I, I)) - precompile(⊗, (I, I, I, I)) - - precompile(FusionStyle, (I,)) - precompile(BraidingStyle, (I,)) - - precompile(dim, (I,)) - precompile(sqrtdim, (I,)) - precompile(invsqrtdim, (I,)) - precompile(dual, (I,)) - precompile(twist, (I,)) - precompile(frobeniusschur, (I,)) - - try - precompile(fusiontensor, (I, I, I)) - catch - end -end diff --git a/TensorKitSectors/src/product.jl b/TensorKitSectors/src/product.jl deleted file mode 100644 index a98d7721..00000000 --- a/TensorKitSectors/src/product.jl +++ /dev/null @@ -1,303 +0,0 @@ -# Deligne tensor product of different sectors: ⊠ -#==============================================================================# -const SectorTuple = Tuple{Vararg{Sector}} - -""" - ProductSector{T<:SectorTuple} - -Represents the Deligne tensor product of sectors. The type parameter `T` is a tuple of the -component sectors. The recommended way to construct a `ProductSector` is using the -[`deligneproduct`](@ref) (`⊠`) operator on the components. -""" -struct ProductSector{T<:SectorTuple} <: Sector - sectors::T -end - -Base.getindex(s::ProductSector, i::Int) = getindex(s.sectors, i) -Base.iterate(s::ProductSector, args...) = iterate(s.sectors, args...) -Base.indexed_iterate(s::ProductSector, args...) = Base.indexed_iterate(s.sectors, args...) - -_sectors(::Type{Tuple{}}) = () -Base.@pure function _sectors(::Type{T}) where {T<:SectorTuple} - return (Base.tuple_type_head(T), _sectors(Base.tuple_type_tail(T))...) -end - -function Base.IteratorSize(::Type{SectorValues{ProductSector{T}}}) where {T<:SectorTuple} - return Base.IteratorSize(Base.Iterators.product(map(values, _sectors(T))...)) -end -function Base.size(::SectorValues{ProductSector{T}}) where {T<:SectorTuple} - return map(s -> length(values(s)), _sectors(T)) -end -Base.length(P::SectorValues{<:ProductSector}) = *(size(P)...) - -function _length(iter::SectorValues{I}) where {I<:Sector} - return Base.IteratorSize(iter) === Base.IsInfinite() ? typemax(Int) : length(iter) -end -function _size(::SectorValues{ProductSector{T}}) where {T<:SectorTuple} - return map(s -> _length(values(s)), _sectors(T)) -end -function Base.getindex(P::SectorValues{ProductSector{T}}, i::Int) where {T<:SectorTuple} - I = manhattan_to_multidimensional_index(i, _size(P)) - return ProductSector{T}(getindex.(values.(_sectors(T)), I)) -end -function findindex(P::SectorValues{ProductSector{T}}, - c::ProductSector{T}) where {T<:SectorTuple} - return to_manhattan_index(findindex.(values.(_sectors(T)), c.sectors), _size(P)) -end - -function Base.iterate(P::SectorValues{ProductSector{T}}, i=1) where {T<:SectorTuple} - Base.IteratorSize(P) != Base.IsInfinite() && i > length(P) && return nothing - return getindex(P, i), i + 1 -end - -ProductSector{T}(args...) where {T<:SectorTuple} = ProductSector{T}(args) -function Base.convert(::Type{ProductSector{T}}, t::Tuple) where {T<:SectorTuple} - return ProductSector{T}(convert(T, t)) -end - -Base.one(::Type{ProductSector{T}}) where {I<:Sector,T<:Tuple{I}} = ProductSector((one(I),)) -function Base.one(::Type{ProductSector{T}}) where {I<:Sector,T<:Tuple{I,Vararg{Sector}}} - return one(I) ⊠ one(ProductSector{Base.tuple_type_tail(T)}) -end - -Base.conj(p::ProductSector) = ProductSector(map(conj, p.sectors)) -function ⊗(p1::P, p2::P) where {P<:ProductSector} - if FusionStyle(P) isa UniqueFusion - (P(first(product(map(⊗, p1.sectors, p2.sectors)...))),) - else - return SectorSet{P}(product(map(⊗, p1.sectors, p2.sectors)...)) - end -end - -function Nsymbol(a::P, b::P, c::P) where {P<:ProductSector} - return prod(map(Nsymbol, a.sectors, b.sectors, c.sectors)) -end - -_firstsector(x::ProductSector) = x.sectors[1] -_tailsector(x::ProductSector) = ProductSector(Base.tail(x.sectors)) - -function Fsymbol(a::P, b::P, c::P, d::P, e::P, f::P) where {P<:ProductSector} - heads = map(_firstsector, (a, b, c, d, e, f)) - tails = map(_tailsector, (a, b, c, d, e, f)) - f₁ = Fsymbol(heads...) - f₂ = Fsymbol(tails...) - if f₁ isa Number || f₂ isa Number - f₁ * f₂ - else - _kron(f₁, f₂) - end -end -function Fsymbol(a::P, b::P, c::P, d::P, e::P, - f::P) where {P<:ProductSector{<:Tuple{Sector}}} - return Fsymbol(map(_firstsector, (a, b, c, d, e, f))...) -end - -function Rsymbol(a::P, b::P, c::P) where {P<:ProductSector} - heads = map(_firstsector, (a, b, c)) - tails = map(_tailsector, (a, b, c)) - r1 = Rsymbol(heads...) - r2 = Rsymbol(tails...) - if r1 isa Number || r2 isa Number - r1 * r2 - else - _kron(r1, r2) - end -end -function Rsymbol(a::P, b::P, c::P) where {P<:ProductSector{<:Tuple{Sector}}} - return Rsymbol(map(_firstsector, (a, b, c))...) -end - -function Bsymbol(a::P, b::P, c::P) where {P<:ProductSector} - heads = map(_firstsector, (a, b, c)) - tails = map(_tailsector, (a, b, c)) - b1 = Bsymbol(heads...) - b2 = Bsymbol(tails...) - if b1 isa Number || b2 isa Number - b1 * b2 - else - _kron(b1, b2) - end -end -function Bsymbol(a::P, b::P, c::P) where {P<:ProductSector{<:Tuple{Sector}}} - return Bsymbol(map(_firstsector, (a, b, c))...) -end - -function Asymbol(a::P, b::P, c::P) where {P<:ProductSector} - heads = map(_firstsector, (a, b, c)) - tails = map(_tailsector, (a, b, c)) - a1 = Asymbol(heads...) - a2 = Asymbol(tails...) - if a1 isa Number || a2 isa Number - a1 * a2 - else - _kron(a1, a2) - end -end -function Asymbol(a::P, b::P, c::P) where {P<:ProductSector{<:Tuple{Sector}}} - return Asymbol(map(_firstsector, (a, b, c))...) -end - -frobeniusschur(p::ProductSector) = prod(map(frobeniusschur, p.sectors)) - -function fusiontensor(a::P, b::P, c::P) where {P<:ProductSector} - return _kron(fusiontensor(map(_firstsector, (a, b, c))...), - fusiontensor(map(_tailsector, (a, b, c))...)) -end - -function fusiontensor(a::P, b::P, c::P) where {P<:ProductSector{<:Tuple{Sector}}} - return fusiontensor(map(_firstsector, (a, b, c))...) -end - -function FusionStyle(::Type{<:ProductSector{T}}) where {T<:SectorTuple} - return Base.:&(map(FusionStyle, _sectors(T))...) -end -function BraidingStyle(::Type{<:ProductSector{T}}) where {T<:SectorTuple} - return Base.:&(map(BraidingStyle, _sectors(T))...) -end -Base.isreal(::Type{<:ProductSector{T}}) where {T<:SectorTuple} = _isreal(T) -_isreal(::Type{Tuple{}}) = true -function _isreal(T::Type{<:SectorTuple}) - return isreal(Base.tuple_type_head(T)) && _isreal(Base.tuple_type_tail(T)) -end - -fermionparity(P::ProductSector) = mapreduce(fermionparity, xor, P.sectors) - -dim(p::ProductSector) = *(dim.(p.sectors)...) - -Base.isequal(p1::ProductSector, p2::ProductSector) = isequal(p1.sectors, p2.sectors) -Base.hash(p::ProductSector, h::UInt) = hash(p.sectors, h) -function Base.isless(p1::P, p2::P) where {P<:ProductSector} - return isless(findindex(values(P), p1), findindex(values(P), p2)) -end - -# Default construction from tensor product of sectors -#----------------------------------------------------- -⊠(s1, s2, s3, s4...) = ⊠(⊠(s1, s2), s3, s4...) -const deligneproduct = ⊠ - -""" - ⊠(s₁::Sector, s₂::Sector) - deligneproduct(s₁::Sector, s₂::Sector) - -Given two sectors `s₁` and `s₂`, which label an isomorphism class of simple objects in a -fusion category ``C₁`` and ``C₂``, `s1 ⊠ s2` (obtained as `\\boxtimes+TAB`) labels the -isomorphism class of simple objects in the Deligne tensor product category ``C₁ ⊠ C₂``. - -The Deligne tensor product also works in the type domain and for spaces and tensors. For -group representations, we have `Irrep[G₁] ⊠ Irrep[G₂] == Irrep[G₁ × G₂]`. -""" -⊠(s1::Sector, s2::Sector) = ProductSector((s1, s2)) -⊠(s1::Trivial, s2::Trivial) = s1 -⊠(s1::Sector, s2::Trivial) = s1 -⊠(s1::Trivial, s2::Sector) = s2 -⊠(p1::ProductSector, s2::Trivial) = p1 -⊠(p1::ProductSector, s2::Sector) = ProductSector(tuple(p1.sectors..., s2)) -⊠(s1::Trivial, p2::ProductSector) = p2 -⊠(s1::Sector, p2::ProductSector) = ProductSector(tuple(s1, p2.sectors...)) -⊠(p1::ProductSector, p2::ProductSector) = ProductSector(tuple(p1.sectors..., p2.sectors...)) - -⊠(I1::Type{Trivial}, I2::Type{Trivial}) = Trivial -⊠(I1::Type{Trivial}, I2::Type{<:ProductSector}) = I2 -⊠(I1::Type{Trivial}, I2::Type{<:Sector}) = I2 - -⊠(I1::Type{<:ProductSector}, I2::Type{Trivial}) = I1 -@static if VERSION >= v"1.8" - Base.@assume_effects :foldable function ⊠(I1::Type{<:ProductSector}, - I2::Type{<:ProductSector}) - T1 = I1.parameters[1] - T2 = I2.parameters[1] - return ProductSector{Tuple{T1.parameters...,T2.parameters...}} - end -else - Base.@pure function ⊠(I1::Type{<:ProductSector}, I2::Type{<:ProductSector}) - T1 = I1.parameters[1] - T2 = I2.parameters[1] - return ProductSector{Tuple{T1.parameters...,T2.parameters...}} - end -end -⊠(I1::Type{<:ProductSector}, I2::Type{<:Sector}) = I1 ⊠ ProductSector{Tuple{I2}} - -⊠(I1::Type{<:Sector}, I2::Type{Trivial}) = I1 -⊠(I1::Type{<:Sector}, I2::Type{<:ProductSector}) = ProductSector{Tuple{I1}} ⊠ I2 -⊠(I1::Type{<:Sector}, I2::Type{<:Sector}) = ProductSector{Tuple{I1,I2}} - -function Base.show(io::IO, P::ProductSector) - sectors = P.sectors - compact = get(io, :typeinfo, nothing) === typeof(P) - sep = compact ? ", " : " ⊠ " - print(io, "(") - for i in 1:length(sectors) - i == 1 || print(io, sep) - io2 = compact ? IOContext(io, :typeinfo => typeof(sectors[i])) : io - print(io2, sectors[i]) - end - return print(io, ")") -end - -function type_repr(P::Type{<:ProductSector}) - sectors = P.parameters[1].parameters - if length(sectors) == 1 - s = "ProductSector{Tuple{" * type_repr(sectors[1]) * "}}" - else - s = "(" - for i in 1:length(sectors) - if i != 1 - s *= " ⊠ " - end - s *= type_repr(sectors[i]) - end - s *= ")" - end - return s -end - -#============================================================================== -TODO: the following would implement pretty-printing of product sectors, i.e. -`ProductSector{Tuple{Irrep[G]}}` would be printed as `Irrep[G]`, and -`ProductSector{Tuple{Irrep[G]}}(x)` would be printed as `Irrep[G](x)`. -However, defining show for a type is considered type piracy/treason, and can lead -to unexpected behavior. While we can avoid this by only defining show for the -instances, this would lead to the following behavior: - -```julia-repl -julia> [Irrep[ℤ₂ × U₁](0, 0)] -1-element Vector{TensorKit.ProductSector{Tuple{Z2Irrep, U1Irrep}}}: - (0, 0) -``` - -See Julia issues #29988, #29428, #22363, #28983. - -Base.show(io::IO, P::Type{<:ProductSector}) = print(io, type_repr(P)) -==============================================================================# -function Base.show(io::IO, P::ProductSector{T}) where {T<:Tuple{Vararg{AbstractIrrep}}} - sectors = P.sectors - get(io, :typeinfo, nothing) === typeof(P) || print(io, type_repr(typeof(P))) - print(io, "(") - for i in 1:length(sectors) - i == 1 || print(io, ", ") - print(IOContext(io, :typeinfo => typeof(sectors[i])), sectors[i]) - end - return print(io, ")") -end - -function type_repr(::Type{ProductSector{T}}) where {T<:Tuple{Vararg{AbstractIrrep}}} - sectors = T.parameters - s = "Irrep[" - for i in 1:length(sectors) - if i != 1 - s *= " × " - end - s *= type_repr(supertype(sectors[i]).parameters[1]) - end - s *= "]" - return s -end - -function Base.getindex(::IrrepTable, ::Type{ProductGroup{Gs}}) where {Gs<:GroupTuple} - G1 = tuple_type_head(Gs) - Grem = tuple_type_tail(Gs) - return ProductSector{Tuple{Irrep[G1]}} ⊠ Irrep[ProductGroup{tuple_type_tail(Gs)}] -end -function Base.getindex(::IrrepTable, ::Type{ProductGroup{Tuple{G}}}) where {G<:Group} - return ProductSector{Tuple{Irrep[G]}} -end diff --git a/TensorKitSectors/src/sectors.jl b/TensorKitSectors/src/sectors.jl deleted file mode 100644 index 47ca56ce..00000000 --- a/TensorKitSectors/src/sectors.jl +++ /dev/null @@ -1,436 +0,0 @@ -""" - abstract type Sector end - -Abstract type for representing the (isomorphism classes of) simple objects in (unitary -and pivotal) (pre-)fusion categories, e.g. the irreducible representations of a finite or -compact group. Subtypes `I<:Sector` as the set of labels of a `GradedSpace`. - -Every new `I<:Sector` should implement the following methods: -* `one(::Type{I})`: unit element of `I` -* `conj(a::I)`: ``a̅``, conjugate or dual label of ``a`` -* `⊗(a::I, b::I)`: iterable with unique fusion outputs of ``a ⊗ b`` - (i.e. don't repeat in case of multiplicities) -* `Nsymbol(a::I, b::I, c::I)`: number of times `c` appears in `a ⊗ b`, i.e. the - multiplicity -* `FusionStyle(::Type{I})`: `UniqueFusion()`, `SimpleFusion()` or - `GenericFusion()` -* `BraidingStyle(::Type{I})`: `Bosonic()`, `Fermionic()`, `Anyonic()`, ... -* `Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I)`: F-symbol: scalar (in case of - `UniqueFusion`/`SimpleFusion`) or matrix (in case of `GenericFusion`) -* `Rsymbol(a::I, b::I, c::I)`: R-symbol: scalar (in case of - `UniqueFusion`/`SimpleFusion`) or matrix (in case of `GenericFusion`) -and optionally -* `dim(a::I)`: quantum dimension of sector `a` -* `frobeniusschur(a::I)`: Frobenius-Schur indicator of `a` -* `Bsymbol(a::I, b::I, c::I)`: B-symbol: scalar (in case of - `UniqueFusion`/`SimpleFusion`) or matrix (in case of `GenericFusion`) -* `twist(a::I)` -> twist of sector `a` -and optionally, if `FusionStyle(I) isa GenericFusion` -* `vertex_ind2label(i::Int, a::I, b::I, c::I)` -> a custom label for the `i`th copy of - `c` appearing in `a ⊗ b` - -Furthermore, `iterate` and `Base.IteratorSize` should be made to work for the singleton type -[`SectorValues{I}`](@ref). -""" -abstract type Sector end - -# iterator over the values (i.e., elements of representative set of simple objects) -# in the sector -""" - struct SectorValues{I<:Sector} - -Singleton type to represent an iterator over the possible values of type `I`, whose -instance is obtained as `values(I)`. For a new `I::Sector`, the following should be defined -* `Base.iterate(::SectorValues{I}[, state])`: iterate over the values -* `Base.IteratorSize(::Type{SectorValues{I}})`: `HasLenght()`, `SizeUnkown()` - or `IsInfinite()` depending on whether the number of values of type `I` is finite - (and sufficiently small) or infinite; for a large number of values, `SizeUnknown()` is - recommend because this will trigger the use of `GenericGradedSpace`. -If `IteratorSize(I) == HasLength()`, also the following must be implemented: -* `Base.length(::SectorValues{I})`: the number of different values -* `Base.getindex(::SectorValues{I}, i::Int)`: a mapping between an index `i` and an - instance of `I` -* `findindex(::SectorValues{I}, c::I)`: reverse mapping between a value `c::I` and an - index `i::Integer ∈ 1:length(values(I))` -""" -struct SectorValues{I<:Sector} end -Base.IteratorEltype(::Type{<:SectorValues}) = HasEltype() -Base.eltype(::Type{SectorValues{I}}) where {I<:Sector} = I -Base.values(::Type{I}) where {I<:Sector} = SectorValues{I}() - -""" - one(::Sector) -> Sector - one(::Type{<:Sector}) -> Sector - -Return the unit element within this type of sector. -""" -Base.one(a::Sector) = one(typeof(a)) - -""" - dual(a::Sector) -> Sector - -Return the conjugate label `conj(a)`. -""" -dual(a::Sector) = conj(a) - -""" - sectorscalartype(I::Type{<:Sector}) -> Type - -Return the scalar type of the topological data (Fsymbol, Rsymbol) of the sector `I`. -""" -function sectorscalartype(::Type{I}) where {I<:Sector} - if BraidingStyle(I) isa NoBraiding - return eltype(Core.Compiler.return_type(Fsymbol, NTuple{6,I})) - else - Feltype = eltype(Core.Compiler.return_type(Fsymbol, NTuple{6,I})) - Reltype = eltype(Core.Compiler.return_type(Rsymbol, NTuple{3,I})) - return Base.promote_op(*, Feltype, Reltype) - end -end - -""" - isreal(::Type{<:Sector}) -> Bool - -Return whether the topological data (Fsymbol, Rsymbol) of the sector is real or not (in -which case it is complex). -""" -Base.isreal(I::Type{<:Sector}) = sectorscalartype(I) <: Real - -# FusionStyle: the most important aspect of Sector -#--------------------------------------------- -""" - ⊗(a::I, b::I...) where {I<:Sector} - otimes(a::I, b::I...) where {I<:Sector} - -Return an iterable of elements of `c::I` that appear in the fusion product `a ⊗ b`. - -Note that every element `c` should appear at most once, fusion degeneracies (if -`FusionStyle(I) == GenericFusion()`) should be accessed via `Nsymbol(a, b, c)`. -""" -function ⊗ end -const otimes = ⊗ - -⊗(I::Sector) = (I,) - -# NOTE: the following inline is extremely important for performance, especially -# in the case of UniqueFusion, because ⊗(...) is computed very often -@inline function ⊗(a::I, b::I, c::I, rest::Vararg{I}) where {I<:Sector} - if FusionStyle(I) isa UniqueFusion - return a ⊗ first(⊗(b, c, rest...)) - else - s = Set{I}() - for d in ⊗(b, c, rest...) - for e in a ⊗ d - push!(s, e) - end - end - return s - end -end - -""" - Nsymbol(a::I, b::I, c::I) where {I<:Sector} -> Integer - -Return an `Integer` representing the number of times `c` appears in the fusion product -`a ⊗ b`. Could be a `Bool` if `FusionStyle(I) == UniqueFusion()` or `SimpleFusion()`. -""" -function Nsymbol end - -# trait to describe the fusion of superselection sectors -""" - FusionStyle(::Sector) - FusionStyle(I::Type{<:Sector}) - -Trait to describe the fusion behavior of sectors of type `I`, which can be either -* `UniqueFusion()`: single fusion output when fusing two sectors; -* `SimpleFusion()`: multiple outputs, but every output occurs at most one, - also known as multiplicity-free (e.g. irreps of ``SU(2)``); -* `GenericFusion()`: multiple outputs that can occur more than once (e.g. irreps - of ``SU(3)``). - -There is an abstract supertype `MultipleFusion` of which both `SimpleFusion` and -`GenericFusion` are subtypes. Furthermore, there is a type alias `MultiplicityFreeFusion` -for those fusion types which do not require muliplicity labels, i.e. -`MultiplicityFreeFusion = Union{UniqueFusion,SimpleFusion}`. -""" -abstract type FusionStyle end -FusionStyle(a::Sector) = FusionStyle(typeof(a)) - -struct UniqueFusion <: FusionStyle end # unique fusion output when fusing two sectors -abstract type MultipleFusion <: FusionStyle end -struct SimpleFusion <: MultipleFusion end # multiple fusion but multiplicity free -struct GenericFusion <: MultipleFusion end # multiple fusion with multiplicities -const MultiplicityFreeFusion = Union{UniqueFusion,SimpleFusion} - -# combine fusion properties of tensor products of sectors -Base.:&(f::F, ::F) where {F<:FusionStyle} = f -Base.:&(f₁::FusionStyle, f₂::FusionStyle) = f₂ & f₁ - -Base.:&(::SimpleFusion, ::UniqueFusion) = SimpleFusion() -Base.:&(::GenericFusion, ::UniqueFusion) = GenericFusion() -Base.:&(::GenericFusion, ::SimpleFusion) = GenericFusion() - -""" - Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I<:Sector} - -Return the F-symbol ``F^{abc}_d`` that associates the two different fusion orders of sectors -`a`, `b` and `c` into an ouput sector `d`, using either an intermediate sector ``a ⊗ b → e`` -or ``b ⊗ c → f``: -``` -a-<-μ-<-e-<-ν-<-d a-<-λ-<-d - ∨ ∨ -> Fsymbol(a,b,c,d,e,f)[μ,ν,κ,λ] ∨ - b c f - v - b-<-κ - ∨ - c -``` -If `FusionStyle(I)` is `UniqueFusion` or `SimpleFusion`, the F-symbol is a number. Otherwise -it is a rank 4 array of size -`(Nsymbol(a, b, e), Nsymbol(e, c, d), Nsymbol(b, c, f), Nsymbol(a, f, d))`. -""" -function Fsymbol end - -# If a I::Sector with `fusion(I) == GenericFusion` fusion wants to have custom vertex -# labels, a specialized method for `vertindex2label` should be added -""" - vertex_ind2label(k::Int, a::I, b::I, c::I) where {I<:Sector} - -Convert the index `k` of the fusion vertex (a,b)->c into a label. For -`FusionStyle(I) == UniqueFusion()` or `FusionStyle(I) == MultipleFusion()`, where every -fusion output occurs only once and `k == 1`, the default is to suppress vertex labels by -setting them equal to `nothing`. For `FusionStyle(I) == GenericFusion()`, the default is to -just use `k`, unless a specialized method is provided. -""" -function vertex_ind2label(k::Int, a::I, b::I, c::I) where {I<:Sector} - return _ind2label(FusionStyle(I), k::Int, a::I, b::I, c::I) -end -_ind2label(::UniqueFusion, k, a, b, c) = nothing -_ind2label(::SimpleFusion, k, a, b, c) = nothing -_ind2label(::GenericFusion, k, a, b, c) = k - -""" - vertex_labeltype(I::Type{<:Sector}) -> Type - -Return the type of labels for the fusion vertices of sectors of type `I`. -""" -vertex_labeltype(I::Type{<:Sector}) = typeof(vertex_ind2label(1, one(I), one(I), one(I))) - -# properties that can be determined in terms of the F symbol -# TODO: find mechanism for returning these numbers with custom type T<:AbstractFloat -""" - dim(a::Sector) - -Return the (quantum) dimension of the sector `a`. -""" -function dim(a::Sector) - if FusionStyle(a) isa UniqueFusion - 1 - elseif FusionStyle(a) isa SimpleFusion - abs(1 / Fsymbol(a, conj(a), a, a, one(a), one(a))) - else - abs(1 / Fsymbol(a, conj(a), a, a, one(a), one(a))[1]) - end -end -sqrtdim(a::Sector) = (FusionStyle(a) isa UniqueFusion) ? 1 : sqrt(dim(a)) -invsqrtdim(a::Sector) = (FusionStyle(a) isa UniqueFusion) ? 1 : inv(sqrt(dim(a))) - -""" - frobeniusschur(a::Sector) - -Return the Frobenius-Schur indicator of a sector `a`. -""" -function frobeniusschur(a::Sector) - if FusionStyle(a) isa UniqueFusion || FusionStyle(a) isa SimpleFusion - sign(Fsymbol(a, conj(a), a, a, one(a), one(a))) - else - sign(Fsymbol(a, conj(a), a, a, one(a), one(a))[1]) - end -end - -# Not necessary -function Asymbol(a::I, b::I, c::I) where {I<:Sector} - if FusionStyle(I) isa UniqueFusion || FusionStyle(I) isa SimpleFusion - (sqrtdim(a) * sqrtdim(b) * invsqrtdim(c)) * - conj(frobeniusschur(a) * Fsymbol(dual(a), a, b, b, one(a), c)) - else - reshape((sqrtdim(a) * sqrtdim(b) * invsqrtdim(c)) * - conj(frobeniusschur(a) * Fsymbol(dual(a), a, b, b, one(a), c)), - (Nsymbol(a, b, c), Nsymbol(dual(a), c, b))) - end -end - -""" - Bsymbol(a::I, b::I, c::I) where {I<:Sector} - -Return the value of ``B^{ab}_c`` which appears in transforming a splitting vertex -into a fusion vertex using the transformation -``` -a -<-μ-<- c a -<-ν-<- c - ∨ -> √(dim(c)/dim(a)) * Bsymbol(a,b,c)[μ,ν] ∧ - b dual(b) -``` -If `FusionStyle(I)` is `UniqueFusion()` or `SimpleFusion()`, the B-symbol is a -number. Otherwise it is a square matrix with row and column size -`Nsymbol(a, b, c) == Nsymbol(c, dual(b), a)`. -""" -function Bsymbol(a::I, b::I, c::I) where {I<:Sector} - if FusionStyle(I) isa UniqueFusion || FusionStyle(I) isa SimpleFusion - (sqrtdim(a) * sqrtdim(b) * invsqrtdim(c)) * Fsymbol(a, b, dual(b), a, c, one(a)) - else - reshape((sqrtdim(a) * sqrtdim(b) * invsqrtdim(c)) * - Fsymbol(a, b, dual(b), a, c, one(a)), - (Nsymbol(a, b, c), Nsymbol(c, dual(b), a))) - end -end - -# Braiding: -#------------------------------------------------- -# trait to describe type to denote how the elementary spaces in a tensor product space -# interact under permutations or actions of the braid group -""" - BraidingStyle(::Sector) -> ::BraidingStyle - BraidingStyle(I::Type{<:Sector}) -> ::BraidingStyle - -Return the type of braiding and twist behavior of sectors of type `I`, which can be either -* `Bosonic()`: symmetric braiding with trivial twist (i.e. identity) -* `Fermionic()`: symmetric braiding with non-trivial twist (squares to identity) -* `Anyonic()`: general ``R_(a,b)^c`` phase or matrix (depending on `SimpleFusion` or - `GenericFusion` fusion) and arbitrary twists - -Note that `Bosonic` and `Fermionic` are subtypes of `SymmetricBraiding`, which means that -braids are in fact equivalent to crossings (i.e. braiding twice is an identity: -`isone(Rsymbol(b,a,c)*Rsymbol(a,b,c)) == true`) and permutations are uniquely defined. -""" -abstract type BraidingStyle end -BraidingStyle(a::Sector) = BraidingStyle(typeof(a)) - -abstract type HasBraiding <: BraidingStyle end -struct NoBraiding <: BraidingStyle end -abstract type SymmetricBraiding <: HasBraiding end # symmetric braiding => actions of permutation group are well defined -struct Bosonic <: SymmetricBraiding end # all twists are one -struct Fermionic <: SymmetricBraiding end # twists one and minus one -struct Anyonic <: HasBraiding end - -Base.:&(b::B, ::B) where {B<:BraidingStyle} = b -Base.:&(B1::BraidingStyle, B2::BraidingStyle) = B2 & B1 -Base.:&(::Bosonic, ::Fermionic) = Fermionic() -Base.:&(::Bosonic, ::Anyonic) = Anyonic() -Base.:&(::Fermionic, ::Anyonic) = Anyonic() -Base.:&(::Bosonic, ::NoBraiding) = NoBraiding() -Base.:&(::Fermionic, ::NoBraiding) = NoBraiding() -Base.:&(::Anyonic, ::NoBraiding) = NoBraiding() - -""" - Rsymbol(a::I, b::I, c::I) where {I<:Sector} - -Returns the R-symbol ``R^{ab}_c`` that maps between ``c → a ⊗ b`` and ``c → b ⊗ a`` as in -``` -a -<-μ-<- c b -<-ν-<- c - ∨ -> Rsymbol(a,b,c)[μ,ν] v - b a -``` -If `FusionStyle(I)` is `UniqueFusion()` or `SimpleFusion()`, the R-symbol is a -number. Otherwise it is a square matrix with row and column size -`Nsymbol(a,b,c) == Nsymbol(b,a,c)`. -""" -function Rsymbol end - -# properties that can be determined in terms of the R symbol - -""" - twist(a::Sector) - -Return the twist of a sector `a` -""" -twist(a::Sector) = sum(dim(b) / dim(a) * tr(Rsymbol(a, a, b)) for b in a ⊗ a) - -# Pentagon and Hexagon equations -#------------------------------------------------------------------------------- -# Consistency equations for F- and R-symbols -function pentagon_equation(a::I, b::I, c::I, d::I; kwargs...) where {I<:Sector} - for f in ⊗(a, b), h in ⊗(c, d) - for g in ⊗(f, c), i in ⊗(b, h) - for e in intersect(⊗(g, d), ⊗(a, i)) - if FusionStyle(I) isa MultiplicityFreeFusion - p1 = Fsymbol(f, c, d, e, g, h) * Fsymbol(a, b, h, e, f, i) - p2 = zero(p1) - for j in ⊗(b, c) - p2 += Fsymbol(a, b, c, g, f, j) * - Fsymbol(a, j, d, e, g, i) * - Fsymbol(b, c, d, i, j, h) - end - else - @tensor p1[λ, μ, ν, κ, ρ, σ] := Fsymbol(f, c, d, e, g, h)[λ, μ, ν, τ] * - Fsymbol(a, b, h, e, f, i)[κ, τ, ρ, σ] - p2 = zero(p1) - for j in ⊗(b, c) - @tensor p2[λ, μ, ν, κ, ρ, σ] += Fsymbol(a, b, c, g, f, j)[κ, λ, α, - β] * - Fsymbol(a, j, d, e, g, i)[β, μ, τ, - σ] * - Fsymbol(b, c, d, i, j, h)[α, τ, ν, - ρ] - end - end - isapprox(p1, p2; kwargs...) || return false - end - end - end - return true -end - -function hexagon_equation(a::I, b::I, c::I; kwargs...) where {I<:Sector} - BraidingStyle(I) isa NoBraiding && - throw(ArgumentError("Hexagon equation only defined for sectors with braiding")) - for e in ⊗(c, a), f in ⊗(c, b) - for d in intersect(⊗(e, b), ⊗(a, f)) - if FusionStyle(I) isa MultiplicityFreeFusion - p1 = Rsymbol(a, c, e) * Fsymbol(a, c, b, d, e, f) * Rsymbol(b, c, f) - p2 = zero(p1) - for g in ⊗(a, b) - p2 += Fsymbol(c, a, b, d, e, g) * Rsymbol(g, c, d) * - Fsymbol(a, b, c, d, g, f) - end - else - @tensor p1[α, β, μ, ν] := Rsymbol(a, c, e)[α, λ] * - Fsymbol(a, c, b, d, e, f)[λ, β, γ, ν] * - Rsymbol(b, c, f)[γ, μ] - p2 = zero(p1) - for g in ⊗(a, b) - @tensor p2[α, β, μ, ν] += Fsymbol(c, a, b, d, e, g)[α, β, δ, - σ] * - Rsymbol(g, c, d)[σ, ψ] * - Fsymbol(a, b, c, d, g, f)[δ, ψ, μ, ν] - end - end - isapprox(p1, p2; kwargs...) || return false - end - end - return true -end - -# SectorSet: -#------------------------------------------------------------------------------- -# Custum generator to represent sets of sectors with type inference -struct SectorSet{I<:Sector,F,S} - f::F - set::S -end -SectorSet{I}(::Type{F}, set::S) where {I<:Sector,F,S} = SectorSet{I,Type{F},S}(F, set) -SectorSet{I}(f::F, set::S) where {I<:Sector,F,S} = SectorSet{I,F,S}(f, set) -SectorSet{I}(set) where {I<:Sector} = SectorSet{I}(identity, set) - -Base.IteratorEltype(::Type{<:SectorSet}) = HasEltype() -Base.IteratorSize(::Type{SectorSet{I,F,S}}) where {I<:Sector,F,S} = Base.IteratorSize(S) - -Base.eltype(::SectorSet{I}) where {I<:Sector} = I -Base.length(s::SectorSet) = length(s.set) -Base.size(s::SectorSet) = size(s.set) - -function Base.iterate(s::SectorSet{I}, args...) where {I<:Sector} - next = iterate(s.set, args...) - next === nothing && return nothing - val, state = next - return convert(I, s.f(val)), state -end diff --git a/TensorKitSectors/src/trivial.jl b/TensorKitSectors/src/trivial.jl deleted file mode 100644 index be60db31..00000000 --- a/TensorKitSectors/src/trivial.jl +++ /dev/null @@ -1,36 +0,0 @@ -""" - Trivial - -Singleton type to represent the trivial sector, i.e. the trivial representation of the -trivial group. This is equivalent to `Rep[ℤ₁]`, or the unit object of the category `Vect` of -ordinary vector spaces. -""" -struct Trivial <: Sector end - -Base.show(io::IO, ::Trivial) = print(io, "Trivial()") - -# iterators -Base.IteratorSize(::Type{SectorValues{Trivial}}) = HasLength() -Base.length(::SectorValues{Trivial}) = 1 -Base.iterate(::SectorValues{Trivial}, i=false) = return i ? nothing : (Trivial(), true) -function Base.getindex(::SectorValues{Trivial}, i::Int) - return i == 1 ? Trivial() : throw(BoundsError(values(Trivial), i)) -end -findindex(::SectorValues{Trivial}, c::Trivial) = 1 - -# basic properties -Base.one(::Type{Trivial}) = Trivial() -Base.conj(::Trivial) = Trivial() - -Base.isreal(::Type{Trivial}) = true -Base.isless(::Trivial, ::Trivial) = false - -# fusion rules -⊗(::Trivial, ::Trivial) = (Trivial(),) -Nsymbol(::Trivial, ::Trivial, ::Trivial) = true -FusionStyle(::Type{Trivial}) = UniqueFusion() -Fsymbol(::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial) = 1 - -# braiding rules -Rsymbol(::Trivial, ::Trivial, ::Trivial) = 1 -BraidingStyle(::Type{Trivial}) = Bosonic() diff --git a/TensorKitSectors/test/newsectors.jl b/TensorKitSectors/test/newsectors.jl deleted file mode 100644 index f84b4c94..00000000 --- a/TensorKitSectors/test/newsectors.jl +++ /dev/null @@ -1,83 +0,0 @@ -""" - module NewSectors - -Bare-bones implementation of a new sector type, `NewSU2Irrep`, which is set to be of -`GenericFusion` type and `Bosonic` braiding style. -""" -module NewSectors - -export NewSU2Irrep - -using HalfIntegers -using WignerSymbols -using TensorKitSectors - -import TensorKitSectors: FusionStyle, BraidingStyle, Nsymbol, Fsymbol, Rsymbol, dim, - fusiontensor, ⊗ - -struct NewSU2Irrep <: Sector - j::HalfInt - function NewSU2Irrep(j) - j >= zero(j) || error("Not a valid SU₂ irrep") - return new(j) - end -end -Base.convert(::Type{NewSU2Irrep}, j::Real) = NewSU2Irrep(j) - -const _su2one = NewSU2Irrep(zero(HalfInt)) -Base.one(::Type{NewSU2Irrep}) = _su2one -Base.conj(s::NewSU2Irrep) = s -function ⊗(s1::NewSU2Irrep, s2::NewSU2Irrep) - return SectorSet{NewSU2Irrep}(abs(s1.j - s2.j):(s1.j + s2.j)) -end - -Base.IteratorSize(::Type{SectorValues{NewSU2Irrep}}) = Base.IsInfinite() -Base.iterate(::SectorValues{NewSU2Irrep}, i=0) = (NewSU2Irrep(half(i)), i + 1) - -FusionStyle(::Type{NewSU2Irrep}) = GenericFusion() -BraidingStyle(::Type{NewSU2Irrep}) = Bosonic() -Base.isreal(::Type{NewSU2Irrep}) = true - -function Nsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep) - return convert(Int, WignerSymbols.δ(sa.j, sb.j, sc.j)) -end - -function Fsymbol(s1::NewSU2Irrep, s2::NewSU2Irrep, s3::NewSU2Irrep, - s4::NewSU2Irrep, s5::NewSU2Irrep, s6::NewSU2Irrep) - n1 = Nsymbol(s1, s2, s5) - n2 = Nsymbol(s5, s3, s4) - n3 = Nsymbol(s2, s3, s6) - n4 = Nsymbol(s1, s6, s4) - f = all(==(_su2one), (s1, s2, s3, s4, s5, s6)) ? 1.0 : - sqrt(dim(s5) * dim(s6)) * WignerSymbols.racahW(Float64, s1.j, s2.j, - s4.j, s3.j, s5.j, s6.j) - return fill(f, (n1, n2, n3, n4)) -end - -function Rsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep) - Nsymbol(sa, sb, sc) > 0 || return fill(0.0, (0, 0)) - return fill(iseven(convert(Int, sa.j + sb.j - sc.j)) ? 1.0 : -1.0, (1, 1)) -end - -dim(s::NewSU2Irrep) = twice(s.j) + 1 - -function fusiontensor(a::NewSU2Irrep, b::NewSU2Irrep, c::NewSU2Irrep) - C = Array{Float64}(undef, dim(a), dim(b), dim(c), 1) - ja, jb, jc = a.j, b.j, c.j - - for kc in 1:dim(c), kb in 1:dim(b), ka in 1:dim(a) - C[ka, kb, kc, 1] = WignerSymbols.clebschgordan(ja, ja + 1 - ka, jb, jb + 1 - kb, jc, - jc + 1 - kc) - end - return C -end - -Base.hash(s::NewSU2Irrep, h::UInt) = hash(s.j, h) -Base.isless(s1::NewSU2Irrep, s2::NewSU2Irrep) = isless(s1.j, s2.j) - -function Base.getindex(::SectorValues{NewSU2Irrep}, i::Int) - return 1 <= i ? NewSU2Irrep(half(i - 1)) : throw(BoundsError(values(NewSU2Irrep), i)) -end -TensorKitSectors.findindex(::SectorValues{NewSU2Irrep}, s::NewSU2Irrep) = twice(s.j) + 1 - -end # module NewSectors diff --git a/TensorKitSectors/test/runtests.jl b/TensorKitSectors/test/runtests.jl deleted file mode 100644 index 9763d559..00000000 --- a/TensorKitSectors/test/runtests.jl +++ /dev/null @@ -1,41 +0,0 @@ -using Test -using TestExtras -using Random -# using TensorKit: TensorKitSectors -using TensorKitSectors -using TensorOperations -using Base.Iterators: take, product -using LinearAlgebra: LinearAlgebra - -const TKS = TensorKitSectors - -include("testsetup.jl") -using .TestSetup -include("newsectors.jl") -using .NewSectors - -const sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewSU2Irrep, - FibonacciAnyon, IsingAnyon, FermionParity, - FermionParity ⊠ FermionParity, - Z3Irrep ⊠ Z4Irrep, FermionParity ⊠ U1Irrep ⊠ SU2Irrep, - FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, NewSU2Irrep ⊠ NewSU2Irrep, - NewSU2Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ NewSU2Irrep, - Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon) - -@testset "$(TensorKitSectors.type_repr(I))" for I in sectorlist - @include("sectors.jl") -end - -@testset "Deligne product" begin - sectorlist′ = (Trivial, sectorlist...) - for I1 in sectorlist′, I2 in sectorlist′ - a = first(smallset(I1)) - b = first(smallset(I2)) - - @constinferred a ⊠ b - @constinferred a ⊠ b ⊠ a - @constinferred a ⊠ b ⊠ a ⊠ b - @constinferred I1 ⊠ I2 - @test typeof(a ⊠ b) == I1 ⊠ I2 - end -end diff --git a/TensorKitSectors/test/sectors.jl b/TensorKitSectors/test/sectors.jl deleted file mode 100644 index 289e18b5..00000000 --- a/TensorKitSectors/test/sectors.jl +++ /dev/null @@ -1,102 +0,0 @@ -Istr = TKS.type_repr(I) -@testset "Sector $Istr: Basic properties" begin - s = (randsector(I), randsector(I), randsector(I)) - @test eval(Meta.parse(sprint(show, I))) == I - @test eval(Meta.parse(TKS.type_repr(I))) == I - @test eval(Meta.parse(sprint(show, s[1]))) == s[1] - @test @constinferred(hash(s[1])) == hash(deepcopy(s[1])) - @test @constinferred(one(s[1])) == @constinferred(one(I)) - @constinferred dual(s[1]) - @constinferred dim(s[1]) - @constinferred frobeniusschur(s[1]) - @constinferred Nsymbol(s...) - @constinferred Rsymbol(s...) - @constinferred Bsymbol(s...) - @constinferred Fsymbol(s..., s...) - it = @constinferred s[1] ⊗ s[2] - @constinferred ⊗(s..., s...) -end -@testset "Sector $Istr: Value iterator" begin - @test eltype(values(I)) == I - sprev = one(I) - for (i, s) in enumerate(values(I)) - @test !isless(s, sprev) # confirm compatibility with sort order - @test s == @constinferred (values(I)[i]) - @test findindex(values(I), s) == i - sprev = s - i >= 10 && break - end - @test one(I) == first(values(I)) - @test (@constinferred findindex(values(I), one(I))) == 1 - for s in smallset(I) - @test (@constinferred values(I)[findindex(values(I), s)]) == s - end -end -if BraidingStyle(I) isa Bosonic && hasfusiontensor(I) - @testset "Sector $Istr: fusion tensor and F-move and R-move" begin - for a in smallset(I), b in smallset(I) - for c in ⊗(a, b) - X1 = permutedims(fusiontensor(a, b, c), (2, 1, 3, 4)) - X2 = fusiontensor(b, a, c) - l = dim(a) * dim(b) * dim(c) - R = LinearAlgebra.transpose(Rsymbol(a, b, c)) - sz = (l, convert(Int, Nsymbol(a, b, c))) - @test reshape(X1, sz) ≈ reshape(X2, sz) * R - end - end - for a in smallset(I), b in smallset(I), c in smallset(I) - for e in ⊗(a, b), f in ⊗(b, c) - for d in intersect(⊗(e, c), ⊗(a, f)) - X1 = fusiontensor(a, b, e) - X2 = fusiontensor(e, c, d) - Y1 = fusiontensor(b, c, f) - Y2 = fusiontensor(a, f, d) - @tensor f1[-1, -2, -3, -4] := conj(Y2[a, f, d, -4]) * - conj(Y1[b, c, f, -3]) * - X1[a, b, e, -1] * X2[e, c, d, -2] - if FusionStyle(I) isa MultiplicityFreeFusion - f2 = fill(Fsymbol(a, b, c, d, e, f) * dim(d), (1, 1, 1, 1)) - else - f2 = Fsymbol(a, b, c, d, e, f) * dim(d) - end - @test isapprox(f1, f2; atol=1e-12, rtol=1e-12) - end - end - end - end -end -@testset "Sector $Istr: Unitarity of F-move" begin - for a in smallset(I), b in smallset(I), c in smallset(I) - for d in ⊗(a, b, c) - es = collect(intersect(⊗(a, b), map(dual, ⊗(c, dual(d))))) - fs = collect(intersect(⊗(b, c), map(dual, ⊗(dual(d), a)))) - if FusionStyle(I) isa MultiplicityFreeFusion - @test length(es) == length(fs) - F = [Fsymbol(a, b, c, d, e, f) for e in es, f in fs] - else - Fblocks = Vector{Any}() - for e in es - for f in fs - Fs = Fsymbol(a, b, c, d, e, f) - push!(Fblocks, - reshape(Fs, - (size(Fs, 1) * size(Fs, 2), - size(Fs, 3) * size(Fs, 4)))) - end - end - F = hvcat(length(fs), Fblocks...) - end - @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12) - end - end -end -@testset "Sector $Istr: Pentagon equation" begin - for a in smallset(I), b in smallset(I), c in smallset(I), d in smallset(I) - @test pentagon_equation(a, b, c, d; atol=1e-12, rtol=1e-12) - end -end -@testset "Sector $Istr: Hexagon equation" begin - for a in smallset(I), b in smallset(I), c in smallset(I) - @test hexagon_equation(a, b, c; atol=1e-12, rtol=1e-12) - end -end diff --git a/TensorKitSectors/test/testsetup.jl b/TensorKitSectors/test/testsetup.jl deleted file mode 100644 index ea76781c..00000000 --- a/TensorKitSectors/test/testsetup.jl +++ /dev/null @@ -1,55 +0,0 @@ -""" - module TestSetup - -Basic utility for testing sectors. -""" -module TestSetup - -export smallset, randsector, hasfusiontensor -export @include - -using TensorKitSectors, Test, Random -using TestExtras -using Base.Iterators: take, product - -smallset(::Type{I}) where {I<:Sector} = take(values(I), 5) -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) - return length(s) > 6 ? rand(s, 6) : s -end -function smallset(::Type{ProductSector{Tuple{I1,I2,I3}}}) where {I1,I2,I3} - iter = product(smallset(I1), smallset(I2), smallset(I3)) - s = collect(i ⊠ j ⊠ k for (i, j, k) in iter if dim(i) * dim(j) * dim(k) <= 6) - return length(s) > 6 ? rand(s, 6) : s -end -function randsector(::Type{I}) where {I<:Sector} - s = collect(smallset(I)) - a = rand(s) - while a == one(a) # don't use trivial label - a = rand(s) - end - return a -end -function hasfusiontensor(I::Type{<:Sector}) - try - fusiontensor(one(I), one(I), one(I)) - return true - catch e - if e isa MethodError - return false - else - rethrow(e) - end - end -end - -# include file as-is in local scope -macro include(filename::String) - dir = dirname(string(__source__.file)) - filepath = joinpath(dir, filename) - source = "quote; " * read(filepath, String) * "; end" - return esc(Meta.parse(source).args[1]) -end - -end # module TestSetup diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 22f5f235..30f423d6 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -101,8 +101,10 @@ const TO = TensorOperations using LRUCache -using HalfIntegers -using WignerSymbols +using TensorKitSectors +import TensorKitSectors: dim, BraidingStyle, FusionStyle, ⊠, ⊗ +import TensorKitSectors: dual, type_repr +import TensorKitSectors: twist using Base: @boundscheck, @propagate_inbounds, OneTo, tail, front, tuple_type_head, tuple_type_tail, tuple_type_cons, @@ -165,18 +167,6 @@ IndexError() = IndexError{Nothing}(nothing) Base.show(io::IO, ::IndexError{Nothing}) = print(io, "IndexError()") Base.show(io::IO, e::IndexError) = print(io, "IndexError(", e.message, ")") -# typerepr -# type_repr(T::Type) = repr(T) - -# Definitions and methods for superselection sectors (quantum numbers) -#---------------------------------------------------------------------- - -include(joinpath(@__DIR__, "..", "TensorKitSectors", "src", "TensorKitSectors.jl")) -using .TensorKitSectors -import .TensorKitSectors: dim, BraidingStyle, FusionStyle, ⊠, ⊗ -import .TensorKitSectors: dual, type_repr -import .TensorKitSectors: twist - # Constructing and manipulating fusion trees and iterators thereof #------------------------------------------------------------------ include("fusiontrees/fusiontrees.jl")