diff --git a/.github/workflows/CI-Sectors.yml b/.github/workflows/CI-Sectors.yml new file mode 100644 index 00000000..b7f8be6c --- /dev/null +++ b/.github/workflows/CI-Sectors.yml @@ -0,0 +1,49 @@ +name: CI-sectors +on: + push: + branches: + - 'master' + - 'main' + - 'release-' + tags: '*' + pull_request: + 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@v1 + - 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/.gitignore b/.gitignore index ae7f128c..5df92310 100644 --- a/.gitignore +++ b/.gitignore @@ -5,4 +5,5 @@ *-old __pycache__ .ipynb* -Manifest.toml \ No newline at end of file +Manifest.toml +.vscode \ No newline at end of file diff --git a/Changelog.md b/Changelog.md new file mode 100644 index 00000000..be534eed --- /dev/null +++ b/Changelog.md @@ -0,0 +1,31 @@ +# Planned changes for v1.0 + +Features that are planned to be implemented before the release of v1.0.0, in no particular order. + +- [ ] Separate `Sectors` module +- [ ] Make `TrivialTensorMap` and `TensorMap` be the same +- [ ] Simplify `TensorMap` type to hide `rowr` and `colr` +- [ ] Change block order in `rowr` / `colr` to speed up particular contractions +- [ ] Make `AdjointTensorMap` generic +- [ ] Rewrite planar operations in order to be AD-compatible +- [ ] Fix rrules for fermionic tensors +- [ ] Fix GPU support +- [ ] Proper threading support +- [ ] Rewrite documentation + +# Changelog + +### `AbstractTensorMap{E,S,N₁,N₂}` + +This adds the scalar type as a parameter to the `AbstractTensorMap` type. This is useful in +the contexts where different types of tensors are used (`AdjointTensor`, `BraidingTensor`, +...), which still have the same scalartype. Additionally, this removes many specializations +for methods in order to reduce ambiguity errors. + +### `copy(BraidingTensor` + +This PR changes the behaviour of `copy` as an instantiator for creating a `TensorMap` from a +`BraidingTensor`. The rationale is that while this does sometimes happen in Julia `Base`, +this is always in the context of lazy wrapper types, for which it makes sense to not copy +the parent and then create a new wrapper. `BraidingTensor` does not wrap anything, so this +definition makes less sense. diff --git a/Project.toml b/Project.toml index 5299ed3f..7edb6cd5 100644 --- a/Project.toml +++ b/Project.toml @@ -27,8 +27,8 @@ ChainRulesTestUtils = "1" Combinatorics = "1" FiniteDifferences = "0.12" HalfIntegers = "1" -LinearAlgebra = "1" LRUCache = "1.0.2" +LinearAlgebra = "1" PackageExtensionCompat = "1" Random = "1" Strided = "2" diff --git a/TensorKitSectors/Project.toml b/TensorKitSectors/Project.toml new file mode 100644 index 00000000..7b915d17 --- /dev/null +++ b/TensorKitSectors/Project.toml @@ -0,0 +1,21 @@ +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 new file mode 100644 index 00000000..2a5b456d --- /dev/null +++ b/TensorKitSectors/src/TensorKitSectors.jl @@ -0,0 +1,68 @@ +# 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 dim, sqrtdim, isqrtdim, 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/src/sectors/anyons.jl b/TensorKitSectors/src/anyons.jl similarity index 100% rename from src/sectors/anyons.jl rename to TensorKitSectors/src/anyons.jl diff --git a/TensorKitSectors/src/auxiliary.jl b/TensorKitSectors/src/auxiliary.jl new file mode 100644 index 00000000..46486da8 --- /dev/null +++ b/TensorKitSectors/src/auxiliary.jl @@ -0,0 +1,15 @@ +# 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 diff --git a/src/sectors/fermions.jl b/TensorKitSectors/src/fermions.jl similarity index 100% rename from src/sectors/fermions.jl rename to TensorKitSectors/src/fermions.jl diff --git a/src/sectors/groups.jl b/TensorKitSectors/src/groups.jl similarity index 87% rename from src/sectors/groups.jl rename to TensorKitSectors/src/groups.jl index eab121fa..2d34e069 100644 --- a/src/sectors/groups.jl +++ b/TensorKitSectors/src/groups.jl @@ -17,11 +17,21 @@ 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{}}}, diff --git a/src/sectors/irreps.jl b/TensorKitSectors/src/irreps.jl similarity index 100% rename from src/sectors/irreps.jl rename to TensorKitSectors/src/irreps.jl diff --git a/TensorKitSectors/src/precompile.jl b/TensorKitSectors/src/precompile.jl new file mode 100644 index 00000000..491595b9 --- /dev/null +++ b/TensorKitSectors/src/precompile.jl @@ -0,0 +1,32 @@ +""" + 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(isqrtdim, (I,)) + precompile(dual, (I,)) + precompile(twist, (I,)) + precompile(frobeniusschur, (I,)) + + try + precompile(fusiontensor, (I, I, I)) + catch + end +end diff --git a/src/sectors/product.jl b/TensorKitSectors/src/product.jl similarity index 98% rename from src/sectors/product.jl rename to TensorKitSectors/src/product.jl index 3bc45526..4331504c 100644 --- a/src/sectors/product.jl +++ b/TensorKitSectors/src/product.jl @@ -74,7 +74,7 @@ function Nsymbol(a::P, b::P, c::P) where {P<:ProductSector} end _firstsector(x::ProductSector) = x.sectors[1] -_tailsector(x::ProductSector) = ProductSector(tail(x.sectors)) +_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)) @@ -196,7 +196,6 @@ group representations, we have `Irrep[G₁] ⊠ Irrep[G₂] == Irrep[G₁ × G ⊠(s1::Sector, p2::ProductSector) = ProductSector(tuple(s1, p2.sectors...)) ⊠(p1::ProductSector, p2::ProductSector) = ProductSector(tuple(p1.sectors..., p2.sectors...)) -# grow types from the left using Base.tuple_type_cons ⊠(I1::Type{Trivial}, I2::Type{Trivial}) = Trivial ⊠(I1::Type{Trivial}, I2::Type{<:ProductSector}) = I2 ⊠(I1::Type{Trivial}, I2::Type{<:Sector}) = I2 diff --git a/src/sectors/sectors.jl b/TensorKitSectors/src/sectors.jl similarity index 90% rename from src/sectors/sectors.jl rename to TensorKitSectors/src/sectors.jl index bfb83312..527c1856 100644 --- a/src/sectors/sectors.jl +++ b/TensorKitSectors/src/sectors.jl @@ -1,6 +1,3 @@ -# Superselection sectors (quantum numbers): -# for defining graded vector spaces and invariant subspaces of tensor products -#==============================================================================# """ abstract type Sector end @@ -61,25 +58,6 @@ Base.IteratorEltype(::Type{<:SectorValues}) = HasEltype() Base.eltype(::Type{SectorValues{I}}) where {I<:Sector} = I Base.values(::Type{I}) where {I<:Sector} = SectorValues{I}() -# Define a sector for ungraded vector spaces -""" - 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()") - -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 - """ one(::Sector) -> Sector one(::Type{<:Sector}) -> Sector @@ -87,7 +65,6 @@ findindex(::SectorValues{Trivial}, c::Trivial) = 1 Return the unit element within this type of sector. """ Base.one(a::Sector) = one(typeof(a)) -Base.one(::Type{Trivial}) = Trivial() """ dual(a::Sector) -> Sector @@ -95,7 +72,6 @@ Base.one(::Type{Trivial}) = Trivial() Return the conjugate label `conj(a)`. """ dual(a::Sector) = conj(a) -Base.conj(::Trivial) = Trivial() """ isreal(::Type{<:Sector}) -> Bool @@ -112,9 +88,6 @@ function Base.isreal(I::Type{<:Sector}) return (eltype(Fsymbol(u, u, u, u, u, u)) <: Real) end end -Base.isreal(::Type{Trivial}) = true - -Base.isless(::Trivial, ::Trivial) = false # FusionStyle: the most important aspect of Sector #--------------------------------------------- @@ -127,10 +100,27 @@ Return an iterable of elements of `c::I` that appear in the fusion product `a Note that every element `c` should appear at most once, fusion degeneracies (if `FusionStyle(I) == GenericFusion()`) should be accessed via `Nsymbol(a, b, c)`. """ -⊗(::Trivial, ::Trivial) = (Trivial(),) -⊗(I::Sector) = (I,) +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 @@ -138,17 +128,13 @@ Return an `Integer` representing the number of times `c` appears in the fusion p `a ⊗ b`. Could be a `Bool` if `FusionStyle(I) == UniqueFusion()` or `SimpleFusion()`. """ function Nsymbol end -Nsymbol(::Trivial, ::Trivial, ::Trivial) = true # trait to describe the fusion of superselection sectors abstract type FusionStyle end -struct UniqueFusion <: FusionStyle # unique fusion output when fusing two sectors -end +struct UniqueFusion <: FusionStyle end # unique fusion output when fusing two sectors abstract type MultipleFusion <: FusionStyle end -struct SimpleFusion <: MultipleFusion # multiple fusion but multiplicity free -end -struct GenericFusion <: MultipleFusion # multiple fusion with multiplicities -end +struct SimpleFusion <: MultipleFusion end # multiple fusion but multiplicity free +struct GenericFusion <: MultipleFusion end # multiple fusion with multiplicities const MultiplicityFreeFusion = Union{UniqueFusion,SimpleFusion} """ @@ -166,23 +152,6 @@ There is an abstract supertype `MultipleFusion` of which both `SimpleFusion` and `MultiplicityFreeFusion = Union{UniqueFusion,SimpleFusion}`. """ FusionStyle(a::Sector) = FusionStyle(typeof(a)) -FusionStyle(::Type{Trivial}) = UniqueFusion() - -# 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 """ Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I<:Sector} @@ -204,7 +173,6 @@ 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 -Fsymbol(::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial) = 1 """ Rsymbol(a::I, b::I, c::I) where {I<:Sector} @@ -220,7 +188,6 @@ number. Otherwise it is a square matrix with row and column size `Nsymbol(a,b,c) == Nsymbol(b,a,c)`. """ function Rsymbol end -Rsymbol(::Trivial, ::Trivial, ::Trivial) = 1 # If a I::Sector with `fusion(I) == GenericFusion` fusion wants to have custom vertex # labels, a specialized method for `vertindex2label` should be added @@ -366,7 +333,6 @@ 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. """ BraidingStyle(a::Sector) = BraidingStyle(typeof(a)) -BraidingStyle(::Type{Trivial}) = Bosonic() # Pentagon and Hexagon equations #------------------------------------------------------------------------------- @@ -457,10 +423,3 @@ function Base.iterate(s::SectorSet{I}, args...) where {I<:Sector} val, state = next return convert(I, s.f(val)), state end - -# possible sectors -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 diff --git a/TensorKitSectors/src/trivial.jl b/TensorKitSectors/src/trivial.jl new file mode 100644 index 00000000..be60db31 --- /dev/null +++ b/TensorKitSectors/src/trivial.jl @@ -0,0 +1,36 @@ +""" + 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/test/newsectors.jl b/TensorKitSectors/test/newsectors.jl similarity index 51% rename from test/newsectors.jl rename to TensorKitSectors/test/newsectors.jl index 4aa5e606..c1bbe803 100644 --- a/test/newsectors.jl +++ b/TensorKitSectors/test/newsectors.jl @@ -1,13 +1,21 @@ -# NewSU2Irrep -# Test of a bare-bones sector implementation, which is set to be `DegenerateNonAbelian` -# (even though it is not) +""" + 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, WignerSymbols, TensorKit +using HalfIntegers +using WignerSymbols +using TensorKitSectors + +import TensorKitSectors: FusionStyle, BraidingStyle, Nsymbol, Fsymbol, Rsymbol, dim, + fusiontensor, ⊗ -struct NewSU2Irrep <: TensorKit.Sector +struct NewSU2Irrep <: Sector j::HalfInt function NewSU2Irrep(j) j >= zero(j) || error("Not a valid SU₂ irrep") @@ -19,27 +27,23 @@ Base.convert(::Type{NewSU2Irrep}, j::Real) = NewSU2Irrep(j) const _su2one = NewSU2Irrep(zero(HalfInt)) Base.one(::Type{NewSU2Irrep}) = _su2one Base.conj(s::NewSU2Irrep) = s -function TensorKit.:⊗(s1::NewSU2Irrep, s2::NewSU2Irrep) - return TensorKit.SectorSet{NewSU2Irrep}(abs(s1.j - s2.j):(s1.j + s2.j)) +function ⊗(s1::NewSU2Irrep, s2::NewSU2Irrep) + return SectorSet{NewSU2Irrep}(abs(s1.j - s2.j):(s1.j + s2.j)) end -Base.IteratorSize(::Type{TensorKit.SectorValues{NewSU2Irrep}}) = Base.IsInfinite() -Base.iterate(::TensorKit.SectorValues{NewSU2Irrep}, i=0) = (NewSU2Irrep(half(i)), i + 1) -# Base.getindex(::SectorValues{NewSU2Irrep}, i::Int) = -# 1 <= i ? NewSU2Irrep(half(i-1)) : throw(BoundsError(values(NewSU2Irrep), i)) -# findindex(::SectorValues{NewSU2Irrep}, s::NewSU2Irrep) = twice(s.j)+1 +Base.IteratorSize(::Type{SectorValues{NewSU2Irrep}}) = Base.IsInfinite() +Base.iterate(::SectorValues{NewSU2Irrep}, i=0) = (NewSU2Irrep(half(i)), i + 1) -# TensorKit.dim(s::NewSU2Irrep) = twice(s.j)+1 -# -TensorKit.FusionStyle(::Type{NewSU2Irrep}) = GenericFusion() -TensorKit.BraidingStyle(::Type{NewSU2Irrep}) = Bosonic() +FusionStyle(::Type{NewSU2Irrep}) = GenericFusion() +BraidingStyle(::Type{NewSU2Irrep}) = Bosonic() Base.isreal(::Type{NewSU2Irrep}) = true -function TensorKit.Nsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep) +function Nsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep) return convert(Int, WignerSymbols.δ(sa.j, sb.j, sc.j)) end -function TensorKit.Fsymbol(s1::NewSU2Irrep, s2::NewSU2Irrep, s3::NewSU2Irrep, - s4::NewSU2Irrep, s5::NewSU2Irrep, s6::NewSU2Irrep) + +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) @@ -49,13 +53,15 @@ function TensorKit.Fsymbol(s1::NewSU2Irrep, s2::NewSU2Irrep, s3::NewSU2Irrep, s4.j, s3.j, s5.j, s6.j) return fill(f, (n1, n2, n3, n4)) end -function TensorKit.Rsymbol(sa::NewSU2Irrep, sb::NewSU2Irrep, sc::NewSU2Irrep) + +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 -TensorKit.dim(s::NewSU2Irrep) = twice(s.j) + 1 -function TensorKit.fusiontensor(a::NewSU2Irrep, b::NewSU2Irrep, c::NewSU2Irrep) +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 @@ -69,4 +75,4 @@ end Base.hash(s::NewSU2Irrep, h::UInt) = hash(s.j, h) Base.isless(s1::NewSU2Irrep, s2::NewSU2Irrep) = isless(s1.j, s2.j) -end +end # module NewSectors diff --git a/TensorKitSectors/test/runtests.jl b/TensorKitSectors/test/runtests.jl new file mode 100644 index 00000000..3e1ddc7a --- /dev/null +++ b/TensorKitSectors/test/runtests.jl @@ -0,0 +1,40 @@ +using Test +using TestExtras +using Random +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 new file mode 100644 index 00000000..53d20ee7 --- /dev/null +++ b/TensorKitSectors/test/sectors.jl @@ -0,0 +1,111 @@ +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 + if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector + @test_throws ArgumentError values(I)[i] + @test_throws ArgumentError findindex(values(I), s) + elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int}) + @test s == @constinferred (values(I)[i]) + @test findindex(values(I), s) == i + end + sprev = s + i >= 10 && break + end + @test one(I) == first(values(I)) + if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector + @test_throws ArgumentError findindex(values(I), one(I)) + elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int}) + @test (@constinferred findindex(values(I), one(I))) == 1 + for s in smallset(I) + @test (@constinferred values(I)[findindex(values(I), s)]) == s + end + 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 new file mode 100644 index 00000000..ea76781c --- /dev/null +++ b/TensorKitSectors/test/testsetup.jl @@ -0,0 +1,55 @@ +""" + 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/ext/TensorKitChainRulesCoreExt.jl b/ext/TensorKitChainRulesCoreExt.jl index b36976fd..cb38d009 100644 --- a/ext/TensorKitChainRulesCoreExt.jl +++ b/ext/TensorKitChainRulesCoreExt.jl @@ -22,8 +22,8 @@ function _repartition(p::Union{IndexTuple,Index2Tuple}, ::Index2Tuple{N₁}) whe return _repartition(p, N₁) end function _repartition(p::Union{IndexTuple,Index2Tuple}, - ::AbstractTensorMap{<:Any,N₁}) where {N₁} - return _repartition(p, N₁) + t::AbstractTensorMap) + return _repartition(p, numout(t)) end TensorKit.block(t::ZeroTangent, c::Sector) = t @@ -59,8 +59,8 @@ function ChainRulesCore.rrule(::typeof(Base.copy), t::AbstractTensorMap) end ChainRulesCore.ProjectTo(::T) where {T<:AbstractTensorMap} = ProjectTo{T}() -function (::ProjectTo{T1})(x::T2) where {S,N1,N2,T1<:AbstractTensorMap{S,N1,N2}, - T2<:AbstractTensorMap{S,N1,N2}} +function (::ProjectTo{T1})(x::T2) where {S,N1,N2,T1<:AbstractTensorMap{<:Any,S,N1,N2}, + T2<:AbstractTensorMap{<:Any,S,N1,N2}} T1 === T2 && return x y = similar(x, scalartype(T1)) for (c, b) in blocks(y) diff --git a/src/TensorKit.jl b/src/TensorKit.jl index a796c9e5..e640f065 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -163,11 +163,16 @@ 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) +# type_repr(T::Type) = repr(T) # Definitions and methods for superselection sectors (quantum numbers) #---------------------------------------------------------------------- -include("sectors/sectors.jl") + +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 #------------------------------------------------------------------ diff --git a/src/planar/planaroperations.jl b/src/planar/planaroperations.jl index ee23bcd7..1d763631 100644 --- a/src/planar/planaroperations.jl +++ b/src/planar/planaroperations.jl @@ -1,34 +1,40 @@ # planar versions of tensor operations add!, trace! and contract! -function planaradd!(C::AbstractTensorMap{S,N₁,N₂}, - A::AbstractTensorMap{S}, - p::Index2Tuple{N₁,N₂}, - α::Number, - β::Number, - backend::Backend...) where {S,N₁,N₂} +function planaradd!(C::AbstractTensorMap, + A::AbstractTensorMap, + p::Index2Tuple, + α::Number, β::Number, + backend::Backend...) return add_transpose!(C, A, p, α, β, backend...) end -function planartrace!(C::AbstractTensorMap{S,N₁,N₂}, - A::AbstractTensorMap{S}, - p::Index2Tuple{N₁,N₂}, - q::Index2Tuple{N₃,N₃}, +function planartrace!(C::AbstractTensorMap, + A::AbstractTensorMap, + (p₁, p₂)::Index2Tuple, + (q₁, q₂)::Index2Tuple, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂,N₃} + backend::Backend...) + (S = spacetype(C)) == spacetype(A) || + throw(SpaceMismatch("incompatible spacetypes")) if BraidingStyle(sectortype(S)) == Bosonic() return trace_permute!(C, A, p, q, α, β, backend...) end + (N₃ = length(q₁)) == length(q₂) || + throw(IndexError("number of trace indices does not match")) + N₁, N₂ = length(p₁), length(p₂) @boundscheck begin - all(i -> space(A, p[1][i]) == space(C, i), 1:N₁) || - throw(SpaceMismatch("trace: A = $(codomain(A))←$(domain(A)), - C = $(codomain(C))←$(domain(C)), p1 = $(p1), p2 = $(p2)")) - all(i -> space(A, p[2][i]) == space(C, N₁ + i), 1:N₂) || - throw(SpaceMismatch("trace: A = $(codomain(A))←$(domain(A)), - C = $(codomain(C))←$(domain(C)), p1 = $(p1), p2 = $(p2)")) - all(i -> space(A, q[1][i]) == dual(space(A, q[2][i])), 1:N₃) || - throw(SpaceMismatch("trace: A = $(codomain(A))←$(domain(A)), - q1 = $(q1), q2 = $(q2)")) + numout(C) == N₁ || throw(IndexError("number of output indices does not match")) + numin(C) == N₂ || throw(IndexError("number of input indices does not match")) + all(i -> space(A, p₁[i]) == space(C, i), 1:N₁) || + throw(SpaceMismatch("trace: A = $(space(A)), + C = $(space(C)), p₁ = $(p₁), p₂ = $(p₂)")) + all(i -> space(A, p₂[i]) == space(C, N₁ + i), 1:N₂) || + throw(SpaceMismatch("trace: A = $(space(A)), + C = $(space(C)), p₁ = $(p₁), p₂ = $(p₂)")) + all(i -> space(A, q₁[i]) == dual(space(A, q₂[i])), 1:N₃) || + throw(SpaceMismatch("trace: A = $(space(A)), + q1 = $(q₁), q2 = $(q₂)")) end if iszero(β) @@ -37,23 +43,24 @@ function planartrace!(C::AbstractTensorMap{S,N₁,N₂}, rmul!(C, β) end for (f₁, f₂) in fusiontrees(A) - for ((f₁′, f₂′), coeff) in planar_trace(f₁, f₂, p..., q...) - TO.tensortrace!(C[f₁′, f₂′], p, A[f₁, f₂], q, :N, α * coeff, true, backend...) + for ((f₁′, f₂′), coeff) in planar_trace(f₁, f₂, p₁, p₂, q₁, q₂) + TO.tensortrace!(C[f₁′, f₂′], (p₁, p₂), A[f₁, f₂], (q₁, q₂), :N, α * coeff, true, + backend...) end end return C end -function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, - A::AbstractTensorMap{S}, +function planarcontract!(C::AbstractTensorMap, + A::AbstractTensorMap, pA::Index2Tuple, - B::AbstractTensorMap{S}, + B::AbstractTensorMap, pB::Index2Tuple, - pAB::Index2Tuple{N₁,N₂}, + pAB::Index2Tuple, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} - if BraidingStyle(sectortype(S)) == Bosonic() + backend::Backend...) + if BraidingStyle(sectortype(C)) == Bosonic() return contract!(C, A, pA, B, pB, pAB, α, β, backend...) end diff --git a/src/spaces/productspace.jl b/src/spaces/productspace.jl index 07b08076..e27be39b 100644 --- a/src/spaces/productspace.jl +++ b/src/spaces/productspace.jl @@ -208,7 +208,7 @@ Base.hash(P::ProductSpace, h::UInt) = hash(P.spaces, h) # Default construction from product of spaces #--------------------------------------------- -⊗(V::ElementarySpace...) = ProductSpace(V...) +⊗(V::ElementarySpace, Vrest::ElementarySpace...) = ProductSpace(V, Vrest...) ⊗(P::ProductSpace) = P function ⊗(P1::ProductSpace{S}, P2::ProductSpace{S}) where {S<:ElementarySpace} return ProductSpace{S}(tuple(P1.spaces..., P2.spaces...)) diff --git a/src/tensors/abstracttensor.jl b/src/tensors/abstracttensor.jl index 55195e4f..447538b4 100644 --- a/src/tensors/abstracttensor.jl +++ b/src/tensors/abstracttensor.jl @@ -3,45 +3,45 @@ # Abstract Tensor type #---------------------- """ - abstract type AbstractTensorMap{S<:IndexSpace, N₁, N₂} end + abstract type AbstractTensorMap{E<:Number, S<:IndexSpace, N₁, N₂} end -Abstract supertype of all tensor maps, i.e. linear maps between tensor products -of vector spaces of type `S<:IndexSpace`. An `AbstractTensorMap` maps from -an input space of type `ProductSpace{S, N₂}` to an output space of type -`ProductSpace{S, N₁}`. +Abstract supertype of all tensor maps, i.e. linear maps between tensor products of vector +spaces of type `S<:IndexSpace`, with element type `E`. An `AbstractTensorMap` maps from an +input space of type `ProductSpace{S, N₂}` to an output space of type `ProductSpace{S, N₁}`. """ -abstract type AbstractTensorMap{S<:IndexSpace,N₁,N₂} end +abstract type AbstractTensorMap{E<:Number,S<:IndexSpace,N₁,N₂} end + """ - AbstractTensor{S<:IndexSpace, N} = AbstractTensorMap{S, N, 0} + AbstractTensor{E,S,N} = AbstractTensorMap{E,S,N,0} -Abstract supertype of all tensors, i.e. elements in the tensor product space -of type `ProductSpace{S, N}`, built from elementary spaces of type `S<:IndexSpace`. +Abstract supertype of all tensors, i.e. elements in the tensor product space of type +`ProductSpace{S, N}`, with element type `E`. -An `AbstractTensor{S, N}` is actually a special case `AbstractTensorMap{S, N, 0}`, -i.e. a tensor map with only a non-trivial output space. +An `AbstractTensor{E, S, N}` is actually a special case `AbstractTensorMap{E, S, N, 0}`, +i.e. a tensor map with only non-trivial output spaces. """ -const AbstractTensor{S<:IndexSpace,N} = AbstractTensorMap{S,N,0} +const AbstractTensor{E,S,N} = AbstractTensorMap{E,S,N,0} # tensor characteristics #------------------------ -Base.eltype(::Union{T,Type{T}}) where {T<:AbstractTensorMap} = scalartype(T) +Base.eltype(::Type{<:AbstractTensorMap{E}}) where {E} = E """ spacetype(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Type{S<:IndexSpace} Return the type of the elementary space `S` of a tensor. """ -spacetype(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = S +spacetype(::Type{<:AbstractTensorMap{E,S}}) where {E,S<:IndexSpace} = S """ sectortype(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Type{I<:Sector} Return the type of sector `I` of a tensor. """ -sectortype(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = sectortype(S) +sectortype(::Type{T}) where {T<:AbstractTensorMap} = sectortype(spacetype(T)) -function InnerProductStyle(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} - return InnerProductStyle(S) +function InnerProductStyle(::Type{T}) where {T<:AbstractTensorMap} + return InnerProductStyle(spacetype(T)) end """ @@ -49,7 +49,7 @@ end Return the type of field `𝕂` of a tensor. """ -field(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = field(S) +field(::Type{T}) where {T<:AbstractTensorMap} = field(spacetype(T)) """ numout(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Int @@ -58,7 +58,7 @@ Return the number of output spaces of a tensor. This is equivalent to the number See also [`numin`](@ref) and [`numind`](@ref). """ -numout(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} = N₁ +numout(::Type{<:AbstractTensorMap{E,S,N₁}}) where {E,S,N₁} = N₁ """ numin(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Int @@ -67,7 +67,7 @@ Return the number of input spaces of a tensor. This is equivalent to the number See also [`numout`](@ref) and [`numind`](@ref). """ -numin(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} = N₂ +numin(::Type{<:AbstractTensorMap{E,S,N₁,N₂}}) where {E,S,N₁,N₂} = N₂ """ numind(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Int @@ -77,7 +77,7 @@ total number of spaces in the domain and codomain of that tensor. See also [`numout`](@ref) and [`numin`](@ref). """ -numind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} = N₁ + N₂ +numind(::Type{T}) where {T<:AbstractTensorMap} = numin(T) + numout(T) function similarstoragetype(TT::Type{<:AbstractTensorMap}, ::Type{T}) where {T} return Core.Compiler.return_type(similar, Tuple{storagetype(TT),Type{T}}) @@ -97,7 +97,7 @@ similarstoragetype(t::AbstractTensorMap, T) = similarstoragetype(typeof(t), T) const order = numind @doc """ - codomain(t::AbstractTensorMap{S,N₁,N₂}, [i::Int]) -> ProductSpace{S,N₁} + codomain(t::AbstractTensorMap{E,S,N₁,N₂}, [i::Int]) -> ProductSpace{S,N₁} Return the codomain of a tensor, i.e. the product space of the output spaces. If `i` is specified, return the `i`-th output space. Implementations should provide `codomain(t)`. @@ -109,7 +109,7 @@ codomain(t::AbstractTensorMap, i) = codomain(t)[i] target(t::AbstractTensorMap) = codomain(t) # categorical terminology @doc """ - domain(t::AbstractTensorMap{S,N₁,N₂}, [i::Int]) -> ProductSpace{S,N₂} + domain(t::AbstractTensorMap{E,S,N₁,N₂}, [i::Int]) -> ProductSpace{S,N₂} Return the domain of a tensor, i.e. the product space of the input spaces. If `i` is specified, return the `i`-th input space. Implementations should provide `domain(t)`. @@ -121,7 +121,7 @@ domain(t::AbstractTensorMap, i) = domain(t)[i] source(t::AbstractTensorMap) = domain(t) # categorical terminology """ - space(t::AbstractTensorMap{S,N₁,N₂}, [i::Int]) -> HomSpace{S,N₁,N₂} + space(t::AbstractTensorMap{E,S,N₁,N₂}, [i::Int]) -> HomSpace{S,N₁,N₂} The index information of a tensor, i.e. the `HomSpace` of its domain and codomain. If `i` is specified, return the `i`-th index space. """ @@ -143,8 +143,8 @@ Return all indices of the codomain of a tensor. See also [`domainind`](@ref) and [`allind`](@ref). """ -function codomainind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} - return ntuple(n -> n, N₁) +function codomainind(::Type{T}) where {T<:AbstractTensorMap} + return ntuple(identity, numout(T)) end codomainind(t::AbstractTensorMap) = codomainind(typeof(t)) @@ -155,8 +155,8 @@ Return all indices of the domain of a tensor. See also [`codomainind`](@ref) and [`allind`](@ref). """ -function domainind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} - return ntuple(n -> N₁ + n, N₂) +function domainind(::Type{T}) where {T<:AbstractTensorMap} + return ntuple(n -> numout(T) + n, numin(T)) end domainind(t::AbstractTensorMap) = domainind(typeof(t)) @@ -167,13 +167,13 @@ Return all indices of a tensor, i.e. the indices of its domain and codomain. See also [`codomainind`](@ref) and [`domainind`](@ref). """ -function allind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} - return ntuple(n -> n, N₁ + N₂) +function allind(::Type{T}) where {T<:AbstractTensorMap} + return ntuple(n -> n, numind(T)) end allind(t::AbstractTensorMap) = allind(typeof(t)) -function adjointtensorindex(::AbstractTensorMap{<:IndexSpace,N₁,N₂}, i) where {N₁,N₂} - return ifelse(i <= N₁, N₂ + i, i - N₁) +function adjointtensorindex(t::AbstractTensorMap, i) + return ifelse(i <= numout(t), numin(t) + i, i - numout(t)) end function adjointtensorindices(t::AbstractTensorMap, indices::IndexTuple) @@ -258,7 +258,7 @@ end # Conversion to Array: #---------------------- # probably not optimized for speed, only for checking purposes -function Base.convert(::Type{Array}, t::AbstractTensorMap{S,N₁,N₂}) where {S,N₁,N₂} +function Base.convert(::Type{Array}, t::AbstractTensorMap) I = sectortype(t) if I === Trivial convert(Array, t[]) diff --git a/src/tensors/adjoint.jl b/src/tensors/adjoint.jl index 43ad110d..073361bf 100644 --- a/src/tensors/adjoint.jl +++ b/src/tensors/adjoint.jl @@ -1,19 +1,19 @@ # AdjointTensorMap: lazy adjoint #==========================================================# """ - struct AdjointTensorMap{S<:IndexSpace, N₁, N₂, ...} <: AbstractTensorMap{S, N₁, N₂} + struct AdjointTensorMap{E, S, N₁, N₂, ...} <: AbstractTensorMap{E, S, N₁, N₂} Specific subtype of [`AbstractTensorMap`](@ref) that is a lazy wrapper for representing the adjoint of an instance of [`TensorMap`](@ref). """ -struct AdjointTensorMap{S<:IndexSpace,N₁,N₂,I<:Sector,A,F₁,F₂} <: - AbstractTensorMap{S,N₁,N₂} - parent::TensorMap{S,N₂,N₁,I,A,F₂,F₁} +struct AdjointTensorMap{E,S,N₁,N₂,I,A,F₁,F₂} <: + AbstractTensorMap{E,S,N₁,N₂} + parent::TensorMap{E,S,N₂,N₁,I,A,F₂,F₁} end #! format: off -const AdjointTrivialTensorMap{S<:IndexSpace,N₁,N₂,A<:DenseMatrix} = - AdjointTensorMap{S,N₁,N₂,Trivial,A,Nothing,Nothing} +const AdjointTrivialTensorMap{E,S,N₁,N₂,A<:DenseMatrix} = + AdjointTensorMap{E,S,N₁,N₂,Trivial,A,Nothing,Nothing} #! format: on # Constructor: construct from taking adjoint of a tensor @@ -29,8 +29,8 @@ domain(t::AdjointTensorMap) = codomain(t.parent) blocksectors(t::AdjointTensorMap) = blocksectors(t.parent) #! format: off -storagetype(::Type{<:AdjointTensorMap{<:IndexSpace,N₁,N₂,Trivial,A}}) where {N₁,N₂,A<:DenseMatrix} = A -storagetype(::Type{<:AdjointTensorMap{<:IndexSpace,N₁,N₂,I,<:SectorDict{I,A}}}) where {N₁, N₂,I<:Sector,A<:DenseMatrix} = A +storagetype(::Type{<:AdjointTrivialTensorMap{E,S,N₁,N₂,A}}) where {E,S,N₁,N₂,A<:DenseMatrix} = A +storagetype(::Type{<:AdjointTensorMap{E,S,N₁,N₂,I,<:SectorDict{I,A}}}) where {E,S,N₁,N₂,I<:Sector,A<:DenseMatrix} = A #! format: on dim(t::AdjointTensorMap) = dim(t.parent) @@ -44,8 +44,8 @@ blocks(t::AdjointTensorMap) = (c => b' for (c, b) in blocks(t.parent)) fusiontrees(::AdjointTrivialTensorMap) = ((nothing, nothing),) fusiontrees(t::AdjointTensorMap) = TensorKeyIterator(t.parent.colr, t.parent.rowr) -function Base.getindex(t::AdjointTensorMap{S,N₁,N₂,I}, - f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {S,N₁,N₂,I} +function Base.getindex(t::AdjointTensorMap{E,S,N₁,N₂,I}, + f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {E,S,N₁,N₂,I} c = f₁.coupled @boundscheck begin c == f₂.coupled || throw(SectorMismatch()) @@ -55,9 +55,9 @@ function Base.getindex(t::AdjointTensorMap{S,N₁,N₂,I}, t.parent.colr[c][f₁]])', (dims(codomain(t), f₁.uncoupled)..., dims(domain(t), f₂.uncoupled)...)) end -@propagate_inbounds function Base.setindex!(t::AdjointTensorMap{S,N₁,N₂}, v, +@propagate_inbounds function Base.setindex!(t::AdjointTensorMap{E,S,N₁,N₂,I}, v, f₁::FusionTree{I,N₁}, - f₂::FusionTree{I,N₂}) where {S,N₁,N₂,I} + f₂::FusionTree{I,N₂}) where {E,S,N₁,N₂,I} return copy!(getindex(t, f₁, f₂), v) end @@ -90,16 +90,16 @@ end function Base.summary(t::AdjointTensorMap) return print("AdjointTensorMap(", codomain(t), " ← ", domain(t), ")") end -function Base.show(io::IO, t::AdjointTensorMap{S}) where {S<:IndexSpace} +function Base.show(io::IO, t::AdjointTensorMap) if get(io, :compact, false) print(io, "AdjointTensorMap(", codomain(t), " ← ", domain(t), ")") return end println(io, "AdjointTensorMap(", codomain(t), " ← ", domain(t), "):") - if sectortype(S) == Trivial + if sectortype(t) === Trivial Base.print_array(io, t[]) println(io) - elseif FusionStyle(sectortype(S)) isa UniqueFusion + elseif FusionStyle(sectortype(t)) isa UniqueFusion for (f₁, f₂) in fusiontrees(t) println(io, "* Data for sector ", f₁.uncoupled, " ← ", f₂.uncoupled, ":") Base.print_array(io, t[f₁, f₂]) diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index 4b8489d6..98a2e08e 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -2,7 +2,7 @@ # special (2,2) tensor that implements a standard braiding operation #====================================================================# """ - struct BraidingTensor{S<:IndexSpace} <: AbstractTensorMap{S, 2, 2} + struct BraidingTensor{E,S<:IndexSpace} <: AbstractTensorMap{E, S, 2, 2} BraidingTensor(V1::S, V2::S, adjoint::Bool=false) where {S<:IndexSpace} Specific subtype of [`AbstractTensorMap`](@ref) for representing the braiding tensor that @@ -11,12 +11,11 @@ braids the first input over the second input; its inverse can be obtained as the It holds that `domain(BraidingTensor(V1, V2)) == V1 ⊗ V2` and `codomain(BraidingTensor(V1, V2)) == V2 ⊗ V1`. """ -struct BraidingTensor{S<:IndexSpace,A} <: AbstractTensorMap{S,2,2} +struct BraidingTensor{E,S} <: AbstractTensorMap{E,S,2,2} V1::S V2::S adjoint::Bool - function BraidingTensor{S,A}(V1::S, V2::S, - adjoint::Bool=false) where {S<:IndexSpace,A<:DenseMatrix} + function BraidingTensor{E,S}(V1::S, V2::S, adjoint::Bool=false) where {E,S<:IndexSpace} for a in sectors(V1) for b in sectors(V2) for c in (a ⊗ b) @@ -25,30 +24,31 @@ struct BraidingTensor{S<:IndexSpace,A} <: AbstractTensorMap{S,2,2} end end end - return new{S,A}(V1, V2, adjoint) + return new{E,S}(V1, V2, adjoint) # partial construction: only construct rowr and colr when needed end end function BraidingTensor(V1::S, V2::S, adjoint::Bool=false) where {S<:IndexSpace} if BraidingStyle(sectortype(S)) isa SymmetricBraiding - return BraidingTensor{S,Matrix{Float64}}(V1, V2, adjoint) + return BraidingTensor{Float64,S}(V1, V2, adjoint) else - return BraidingTensor{S,Matrix{ComplexF64}}(V1, V2, adjoint) + return BraidingTensor{ComplexF64,S}(V1, V2, adjoint) end end -function BraidingTensor(V::HomSpace{S}, adjoint::Bool=false) where {S<:IndexSpace} +function BraidingTensor(V::HomSpace, adjoint::Bool=false) domain(V) == reverse(codomain(V)) || throw(SpaceMismatch("Cannot define a braiding on $V")) return BraidingTensor(V[1], V[2], adjoint) end -function Base.adjoint(b::BraidingTensor{S,A}) where {S<:IndexSpace,A<:DenseMatrix} - return BraidingTensor{S,A}(b.V1, b.V2, !b.adjoint) +function Base.adjoint(b::BraidingTensor{E,S}) where {E,S} + return BraidingTensor{E,S}(b.V1, b.V2, !b.adjoint) end domain(b::BraidingTensor) = b.adjoint ? b.V2 ⊗ b.V1 : b.V1 ⊗ b.V2 codomain(b::BraidingTensor) = b.adjoint ? b.V1 ⊗ b.V2 : b.V2 ⊗ b.V1 -storagetype(::Type{BraidingTensor{S,A}}) where {S<:IndexSpace,A<:DenseMatrix} = A +# TODO: check if this can be removed/is correct +storagetype(::Type{BraidingTensor{E,S}}) where {E,S} = Matrix{E} blocksectors(b::BraidingTensor) = blocksectors(b.V1 ⊗ b.V2) hasblock(b::BraidingTensor, s::Sector) = s ∈ blocksectors(b) @@ -87,8 +87,8 @@ function fusiontrees(b::BraidingTensor) return TensorKeyIterator(rowr, colr) end -function Base.getindex(b::BraidingTensor{S}) where {S} - sectortype(S) == Trivial || throw(SectorMismatch()) +function Base.getindex(b::BraidingTensor) + sectortype(b) === Trivial || throw(SectorMismatch()) (V1, V2) = domain(b) d = (dim(V2), dim(V1), dim(V1), dim(V2)) return sreshape(StridedView(block(b, Trivial())), d) @@ -127,7 +127,9 @@ end Base.similar(::BraidingTensor, T::Type, P::TensorMapSpace) = TensorMap(undef, T, P) -Base.copy(b::BraidingTensor) = copy!(similar(b), b) +# efficient copy constructor +Base.copy(b::BraidingTensor) = BraidingTensor(copy(b.V1), copy(b.V2), b.adjoint) + function Base.copy!(t::TensorMap, b::BraidingTensor) space(t) == space(b) || throw(SectorMismatch()) fill!(t, zero(scalartype(t))) @@ -148,8 +150,8 @@ function Base.copy!(t::TensorMap, b::BraidingTensor) end return t end -TensorMap(b::BraidingTensor) = copy(b) -Base.convert(::Type{TensorMap}, b::BraidingTensor) = copy(b) +TensorMap(b::BraidingTensor) = copy!(similar(b), b) +Base.convert(::Type{TensorMap}, b::BraidingTensor) = TensorMap(b) function block(b::BraidingTensor, s::Sector) sectortype(b) == typeof(s) || throw(SectorMismatch()) @@ -192,14 +194,15 @@ blocks(b::BraidingTensor) = blocks(TensorMap(b)) # Index manipulations # ------------------- has_shared_permute(t::BraidingTensor, args...) = false -function add_transform!(tdst::AbstractTensorMap{S,N₁,N₂}, - tsrc::BraidingTensor{S}, - (p₁, p₂)::Index2Tuple{N₁,N₂}, +function add_transform!(tdst::AbstractTensorMap, + tsrc::BraidingTensor, + (p₁, p₂)::Index2Tuple, fusiontreetransform, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} - return add_transform!(tdst, copy(tsrc), (p₁, p₂), fusiontreetransform, α, β, backend...) + backend::Backend...) + return add_transform!(tdst, TensorMap(tsrc), (p₁, p₂), fusiontreetransform, α, β, + backend...) end # VectorInterface @@ -209,30 +212,37 @@ end # TensorOperations # ---------------- # TODO: implement specialized methods -function TO.tensoradd!(C::AbstractTensorMap{S,N₁,N₂}, pC::Index2Tuple{N₁,N₂}, - A::BraidingTensor{S}, conjA::Symbol, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} - return TO.tensoradd!(C, pC, copy(A), conjA, α, β, backend...) + +function TO.tensoradd!(C::AbstractTensorMap, pC::Index2Tuple, + A::BraidingTensor, conjA::Symbol, α::Number, β::Number, + backend::Backend...) + return TO.tensoradd!(C, pC, TensorMap(A), conjA, α, β, backend...) end # Planar operations # ----------------- -function planaradd!(C::AbstractTensorMap{S,N₁,N₂}, - A::BraidingTensor{S}, - p::Index2Tuple{N₁,N₂}, +# TODO: implement specialized methods + +function planaradd!(C::AbstractTensorMap, + A::BraidingTensor, p::Index2Tuple, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} - return planaradd!(C, copy(A), p, α, β, backend...) + backend::Backend...) + return planaradd!(C, TensorMap(A), p, α, β, backend...) end -function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, - A::BraidingTensor{S}, - (oindA, cindA)::Index2Tuple{2,2}, - B::AbstractTensorMap{S}, - (cindB, oindB)::Index2Tuple{2,N₃}, - (p1, p2)::Index2Tuple{N₁,N₂}, +function planarcontract!(C::AbstractTensorMap, + A::BraidingTensor, + (oindA, cindA)::Index2Tuple, + B::AbstractTensorMap, + (cindB, oindB)::Index2Tuple, + (p1, p2)::Index2Tuple, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂,N₃} + backend::Backend...) + # special case only defined for contracting 2 indices + length(oindA) == length(cindA) == 2 || + return planarcontract!(C, TensorMap(A), (oindA, cindA), B, (cindB, oindB), (p1, p2), + α, β, backend...) + codA, domA = codomainind(A), domainind(A) codB, domB = codomainind(B), domainind(B) oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, @@ -271,14 +281,19 @@ function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, end return C end -function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, - A::AbstractTensorMap{S}, - (oindA, cindA)::Index2Tuple{N₃,2}, - B::BraidingTensor{S}, - (cindB, oindB)::Index2Tuple{2,2}, - (p1, p2)::Index2Tuple{N₁,N₂}, +function planarcontract!(C::AbstractTensorMap, + A::AbstractTensorMap, + (oindA, cindA)::Index2Tuple, + B::BraidingTensor, + (cindB, oindB)::Index2Tuple, + (p1, p2)::Index2Tuple, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂,N₃} + backend::Backend...) + # special case only defined for contracting 2 indices + length(oindB) == length(cindB) == 2 || + return planarcontract!(C, A, (oindA, cindA), TensorMap(B), (cindB, oindB), (p1, p2), + α, β, backend...) + codA, domA = codomainind(A), domainind(A) codB, domB = codomainind(B), domainind(B) oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, @@ -317,42 +332,20 @@ function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, end return C end -function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, - A::BraidingTensor{S}, - (oindA, cindA)::Index2Tuple{2,2}, - B::BraidingTensor{S}, - (cindB, oindB)::Index2Tuple{2,2}, - (p1, p2)::Index2Tuple{N₁,N₂}, - α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} - return planarcontract!(C, copy(A), (oindA, cindA), B, (cindB, oindB), (p1, p2), α, β, - backend...) -end -# Fallback cases for planarcontract! -# TODO: implement specialised cases for contracting 0, 1, 3 and 4 indices -function planarcontract!(C::AbstractTensorMap{S}, A::BraidingTensor{S}, pA::Index2Tuple, - B::BraidingTensor{S}, pB::Index2Tuple, α::Number, β::Number, - backend::Backend...) where {S} - return planarcontract!(C, copy(A), pA, copy(B), pB, α, β, backend...) -end -function planarcontract!(C::AbstractTensorMap{S}, A::BraidingTensor{S}, pA::Index2Tuple, - B::AbstractTensorMap{S}, pB::Index2Tuple, α::Number, β::Number, - backend::Backend...) where {S} - return planarcontract!(C, copy(A), pA, B, pB, α, β, backend...) -end -function planarcontract!(C::AbstractTensorMap{S}, A::AbstractTensorMap{S}, pA::Index2Tuple, - B::BraidingTensor{S}, pB::Index2Tuple, α::Number, β::Number, - backend::Backend...) where {S} - return planarcontract!(C, A, pA, copy(B), pB, α, β, backend...) +# ambiguity fix: +function planarcontract!(C::AbstractTensorMap, A::BraidingTensor, pA::Index2Tuple, + B::BraidingTensor, pB::Index2Tuple, pC::Index2Tuple, + α::Number, β::Number, backend::Backend...) + return planarcontract!(C, TensorMap(A), pA, TensorMap(B), pB, pC, α, β, backend...) end -function planartrace!(C::AbstractTensorMap{S,N₁,N₂}, - A::BraidingTensor{S}, - p::Index2Tuple{N₁,N₂}, q::Index2Tuple{N₃,N₃}, +function planartrace!(C::AbstractTensorMap, + A::BraidingTensor, + p::Index2Tuple, q::Index2Tuple, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂,N₃} - return planartrace!(C, copy(A), p, q, α, β, backend...) + backend::Backend...) + return planartrace!(C, TensorMap(A), p, q, α, β, backend...) end # function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl index dccb38a2..d6ebd3ea 100644 --- a/src/tensors/indexmanipulations.jl +++ b/src/tensors/indexmanipulations.jl @@ -1,8 +1,7 @@ # Index manipulations #--------------------- """ - permute!(tdst::AbstractTensorMap{S,N₁,N₂}, tsrc::AbstractTensorMap{S}, - (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}}) where {S,N₁,N₂} + permute!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple) -> tdst Write into `tdst` the result of permuting the indices of `tsrc`. @@ -10,16 +9,16 @@ The codomain and domain of `tdst` correspond to the indices in `p₁` and `p₂` See [`permute`](@ref) for creating a new tensor and [`add_permute!`](@ref) for a more general version. """ -@propagate_inbounds function Base.permute!(tdst::AbstractTensorMap{S,N₁,N₂}, - tsrc::AbstractTensorMap{S}, - p::Index2Tuple{N₁,N₂}) where {S,N₁,N₂} +@propagate_inbounds function Base.permute!(tdst::AbstractTensorMap, + tsrc::AbstractTensorMap, + p::Index2Tuple) return add_permute!(tdst, tsrc, p, true, false) end """ - permute(tsrc::AbstractTensorMap{S}, (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}}; - copy::Bool=false) where {S,N₁,N₂} - -> tdst::TensorMap{S,N₁,N₂} + permute(tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple; + copy::Bool=false) + -> tdst::TensorMap Return tensor `tdst` obtained by permuting the indices of `tsrc`. The codomain and domain of `tdst` correspond to the indices in `p₁` and `p₂` of `tsrc` respectively. @@ -28,10 +27,10 @@ If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwis To permute into an existing destination, see [permute!](@ref) and [`add_permute!`](@ref) """ -function permute(t::AbstractTensorMap{S}, (p₁, p₂)::Index2Tuple{N₁,N₂}; - copy::Bool=false) where {S,N₁,N₂} - cod = ProductSpace{S,N₁}(map(n -> space(t, n), p₁)) - dom = ProductSpace{S,N₂}(map(n -> dual(space(t, n)), p₂)) +function permute(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂}; + copy::Bool=false) where {N₁,N₂} + cod = ProductSpace{spacetype(t),N₁}(map(n -> space(t, n), p₁)) + dom = ProductSpace{spacetype(t),N₂}(map(n -> dual(space(t, n)), p₂)) # share data if possible if !copy if p₁ === codomainind(t) && p₂ === domainind(t) @@ -45,14 +44,13 @@ function permute(t::AbstractTensorMap{S}, (p₁, p₂)::Index2Tuple{N₁,N₂}; return permute!(similar(t, cod ← dom), t, (p₁, p₂)) end end -function permute(t::AdjointTensorMap{S}, (p₁, p₂)::Index2Tuple; - copy::Bool=false) where {S} +function permute(t::AdjointTensorMap, (p₁, p₂)::Index2Tuple; copy::Bool=false) p₁′ = adjointtensorindices(t, p₂) p₂′ = adjointtensorindices(t, p₁) - return adjoint(permute(adjoint(t), (p₁′, p₂′); copy=copy)) + return adjoint(permute(adjoint(t), (p₁′, p₂′); copy)) end function permute(t::AbstractTensorMap, p::IndexTuple; copy::Bool=false) - return permute(t, (p, ()); copy=copy) + return permute(t, (p, ()); copy) end function has_shared_permute(t::TensorMap, (p₁, p₂)::Index2Tuple) @@ -76,8 +74,8 @@ end # Braid """ - braid!(tdst::AbstractTensorMap{S,N₁,N₂}, tsrc::AbstractTensorMap{S}, - (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}}, levels::Tuple) where {S,N₁,N₂} + braid!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, + (p₁, p₂)::Index2Tuple, levels::Tuple) -> tdst Write into `tdst` the result of braiding the indices of `tsrc`. @@ -87,17 +85,17 @@ which determines whether they will braid over or under any other index with whic See [`braid`](@ref) for creating a new tensor and [`add_braid!`](@ref) for a more general version. """ -@propagate_inbounds function braid!(tdst::AbstractTensorMap{S,N₁,N₂}, - tsrc::AbstractTensorMap{S}, - (p₁, p₂)::Index2Tuple{N₁,N₂}, - levels::IndexTuple) where {S,N₁,N₂} - return add_braid!(tdst, tsrc, (p₁, p₂), levels, true, false) +@propagate_inbounds function braid!(tdst::AbstractTensorMap, + tsrc::AbstractTensorMap, + p::Index2Tuple, + levels::IndexTuple) + return add_braid!(tdst, tsrc, p, levels, true, false) end """ - braid(tsrc::AbstractTensorMap{S}, (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}}, levels::Tuple; - copy::Bool = false) where {S,N₁,N₂} - -> tdst::TensorMap{S,N₁,N₂} + braid(tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple, levels::IndexTuple; + copy::Bool = false) + -> tdst::TensorMap Return tensor `tdst` obtained by braiding the indices of `tsrc`. The codomain and domain of `tdst` correspond to the indices in `p₁` and `p₂` of `tsrc` respectively. @@ -108,18 +106,18 @@ If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwis To braid into an existing destination, see [braid!](@ref) and [`add_braid!`](@ref) """ -function braid(t::AbstractTensorMap{S}, (p₁, p₂)::Index2Tuple, levels::IndexTuple; - copy::Bool=false) where {S} +function braid(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple, levels::IndexTuple; + copy::Bool=false) @assert length(levels) == numind(t) - if BraidingStyle(sectortype(S)) isa SymmetricBraiding + if BraidingStyle(sectortype(t)) isa SymmetricBraiding return permute(t, (p₁, p₂); copy=copy) end if !copy && p₁ == codomainind(t) && p₂ == domainind(t) return t end # general case - cod = ProductSpace{S}(map(n -> space(t, n), p₁)) - dom = ProductSpace{S}(map(n -> dual(space(t, n)), p₂)) + cod = ProductSpace{spacetype(t)}(map(n -> space(t, n), p₁)) + dom = ProductSpace{spacetype(t)}(map(n -> dual(space(t, n)), p₂)) @inbounds begin return braid!(similar(t, cod ← dom), t, (p₁, p₂), levels) end @@ -130,8 +128,8 @@ end _transpose_indices(t::AbstractTensorMap) = (reverse(domainind(t)), reverse(codomainind(t))) """ - transpose!(tdst::AbstractTensorMap{S,N₁,N₂}, tsrc::AbstractTensorMap{S}, - (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}}) where {S,N₁,N₂} + transpose!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, + (p₁, p₂)::Index2Tuple) -> tdst Write into `tdst` the result of transposing the indices of `tsrc`. @@ -148,9 +146,9 @@ function LinearAlgebra.transpose!(tdst::AbstractTensorMap, end """ - transpose(tsrc::AbstractTensorMap{S}, (p₁, p₂)::Tuple{IndexTuple{N₁},IndexTuple{N₂}}; - copy::Bool=false) where {S,N₁,N₂} - -> tdst::TensorMap{S,N₁,N₂} + transpose(tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple; + copy::Bool=false) + -> tdst::TensorMap Return tensor `tdst` obtained by transposing the indices of `tsrc`. The codomain and domain of `tdst` correspond to the indices in `p₁` and `p₂` of `tsrc` respectively. @@ -161,26 +159,26 @@ If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwis To permute into an existing destination, see [permute!](@ref) and [`add_permute!`](@ref) """ -function LinearAlgebra.transpose(t::AbstractTensorMap{S}, +function LinearAlgebra.transpose(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple=_transpose_indices(t); - copy::Bool=false) where {S} - if sectortype(S) === Trivial + copy::Bool=false) + if sectortype(t) === Trivial return permute(t, (p₁, p₂); copy=copy) end if !copy && p₁ == codomainind(t) && p₂ == domainind(t) return t end # general case - cod = ProductSpace{S}(map(n -> space(t, n), p₁)) - dom = ProductSpace{S}(map(n -> dual(space(t, n)), p₂)) + cod = ProductSpace{spacetype(t)}(map(n -> space(t, n), p₁)) + dom = ProductSpace{spacetype(t)}(map(n -> dual(space(t, n)), p₂)) @inbounds begin return transpose!(similar(t, cod ← dom), t, (p₁, p₂)) end end -function LinearAlgebra.transpose(t::AdjointTensorMap{S}, +function LinearAlgebra.transpose(t::AdjointTensorMap, (p₁, p₂)::Index2Tuple=_transpose_indices(t); - copy::Bool=false) where {S} + copy::Bool=false) p₁′ = map(n -> adjointtensorindex(t, n), p₂) p₂′ = map(n -> adjointtensorindex(t, n), p₁) return adjoint(transpose(adjoint(t), (p₁′, p₂′); copy=copy)) @@ -268,23 +266,23 @@ twist(t::AbstractTensorMap, i::Int; inv::Bool=false) = twist!(copy(t), i; inv=in #------------------------------------- # Full implementations based on `add` #------------------------------------- -@propagate_inbounds function add_permute!(tdst::AbstractTensorMap{S,N₁,N₂}, +@propagate_inbounds function add_permute!(tdst::AbstractTensorMap{E,S,N₁,N₂}, tsrc::AbstractTensorMap, p::Index2Tuple{N₁,N₂}, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} + backend::Backend...) where {E,S,N₁,N₂} treepermuter(f₁, f₂) = permute(f₁, f₂, p[1], p[2]) return add_transform!(tdst, tsrc, p, treepermuter, α, β, backend...) end -@propagate_inbounds function add_braid!(tdst::AbstractTensorMap{S,N₁,N₂}, +@propagate_inbounds function add_braid!(tdst::AbstractTensorMap{E,S,N₁,N₂}, tsrc::AbstractTensorMap, p::Index2Tuple{N₁,N₂}, levels::IndexTuple, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} + backend::Backend...) where {E,S,N₁,N₂} length(levels) == numind(tsrc) || throw(ArgumentError("incorrect levels $levels for tensor map $(codomain(tsrc)) ← $(domain(tsrc))")) @@ -295,23 +293,23 @@ end return add_transform!(tdst, tsrc, p, treebraider, α, β, backend...) end -@propagate_inbounds function add_transpose!(tdst::AbstractTensorMap{S,N₁,N₂}, +@propagate_inbounds function add_transpose!(tdst::AbstractTensorMap{E,S,N₁,N₂}, tsrc::AbstractTensorMap, p::Index2Tuple{N₁,N₂}, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} + backend::Backend...) where {E,S,N₁,N₂} treetransposer(f₁, f₂) = transpose(f₁, f₂, p[1], p[2]) return add_transform!(tdst, tsrc, p, treetransposer, α, β, backend...) end -function add_transform!(tdst::AbstractTensorMap{S,N₁,N₂}, +function add_transform!(tdst::AbstractTensorMap{E,S,N₁,N₂}, tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂}, fusiontreetransform, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} + backend::Backend...) where {E,S,N₁,N₂} @boundscheck begin all(i -> space(tsrc, p₁[i]) == space(tdst, i), 1:N₁) || throw(SpaceMismatch("source = $(codomain(tsrc))←$(domain(tsrc)), diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index 4678b780..5b79d6e2 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -164,8 +164,8 @@ function Base.fill!(t::AbstractTensorMap, value::Number) end return t end -function LinearAlgebra.adjoint!(tdst::AbstractTensorMap{S}, - tsrc::AbstractTensorMap{S}) where {S} +function LinearAlgebra.adjoint!(tdst::AbstractTensorMap, + tsrc::AbstractTensorMap) InnerProductStyle(tdst) === EuclideanProduct() || throw_invalid_innerproduct(:adjoint!) space(tdst) == adjoint(space(tsrc)) || throw(SpaceMismatch("$(space(tdst)) ≠ adjoint($(space(tsrc)))")) @@ -383,8 +383,7 @@ for f in (:sqrt, :log, :asin, :acos, :acosh, :atanh, :acoth) end # concatenate tensors -function catdomain(t1::AbstractTensorMap{S,N₁,1}, - t2::AbstractTensorMap{S,N₁,1}) where {S,N₁} +function catdomain(t1::T, t2::T) where {S,N₁,T<:AbstractTensorMap{<:Any,S,N₁,1}} codomain(t1) == codomain(t2) || throw(SpaceMismatch("codomains of tensors to concatenate must match:\n" * "$(codomain(t1)) ≠ $(codomain(t2))")) @@ -401,8 +400,7 @@ function catdomain(t1::AbstractTensorMap{S,N₁,1}, end return t end -function catcodomain(t1::AbstractTensorMap{S,1,N₂}, - t2::AbstractTensorMap{S,1,N₂}) where {S,N₂} +function catcodomain(t1::T, t2::T) where {S,N₂,T<:AbstractTensorMap{<:Any,S,1,N₂}} domain(t1) == domain(t2) || throw(SpaceMismatch("domains of tensors to concatenate must match:\n" * "$(domain(t1)) ≠ $(domain(t2))")) @@ -422,14 +420,16 @@ end # tensor product of tensors """ - ⊗(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}, ...) -> TensorMap{S} - otimes(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}, ...) -> TensorMap{S} + ⊗(t1::AbstractTensorMap, t2::AbstractTensorMap, ...) -> TensorMap + otimes(t1::AbstractTensorMap, t2::AbstractTensorMap, ...) -> TensorMap Compute the tensor product between two `AbstractTensorMap` instances, which results in a new `TensorMap` instance whose codomain is `codomain(t1) ⊗ codomain(t2)` and whose domain is `domain(t1) ⊗ domain(t2)`. """ -function ⊗(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}) where {S} +function ⊗(t1::AbstractTensorMap, t2::AbstractTensorMap) + (S = spacetype(t1)) === spacetype(t2) || + throw(SpaceMismatch("spacetype(t1) ≠ spacetype(t2)")) cod1, cod2 = codomain(t1), codomain(t2) dom1, dom2 = domain(t1), domain(t2) cod = cod1 ⊗ cod2 diff --git a/src/tensors/tensor.jl b/src/tensors/tensor.jl index d0538949..82daa4bf 100644 --- a/src/tensors/tensor.jl +++ b/src/tensors/tensor.jl @@ -3,65 +3,67 @@ #==========================================================# #! format: off """ - struct TensorMap{S<:IndexSpace, N₁, N₂, ...} <: AbstractTensorMap{S, N₁, N₂} + struct TensorMap{E, S<:IndexSpace, N₁, N₂, ...} <: AbstractTensorMap{E, S, N₁, N₂} Specific subtype of [`AbstractTensorMap`](@ref) for representing tensor maps (morphisms in a tensor category) whose data is stored in blocks of some subtype of `DenseMatrix`. """ -struct TensorMap{S<:IndexSpace, N₁, N₂, I<:Sector, A<:Union{<:DenseMatrix,SectorDict{I,<:DenseMatrix}}, F₁, F₂} <: AbstractTensorMap{S, N₁, N₂} +struct TensorMap{E, S<:IndexSpace, N₁, N₂, I<:Sector, A<:Union{<:DenseMatrix{E},SectorDict{I,<:DenseMatrix{E}}}, + F₁, F₂} <: AbstractTensorMap{E, S, N₁, N₂} data::A codom::ProductSpace{S,N₁} dom::ProductSpace{S,N₂} rowr::SectorDict{I,FusionTreeDict{F₁,UnitRange{Int}}} colr::SectorDict{I,FusionTreeDict{F₂,UnitRange{Int}}} - function TensorMap{S, N₁, N₂, I, A, F₁, F₂}(data::A, + function TensorMap{E, S, N₁, N₂, I, A, F₁, F₂}(data::A, codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂}, rowr::SectorDict{I,FusionTreeDict{F₁,UnitRange{Int}}}, colr::SectorDict{I,FusionTreeDict{F₂,UnitRange{Int}}}) where - {S<:IndexSpace, N₁, N₂, I<:Sector, A<:SectorDict{I,<:DenseMatrix}, + {E,S<:IndexSpace, N₁, N₂, I<:Sector, A<:SectorDict{I,<:DenseMatrix{E}}, F₁<:FusionTree{I,N₁}, F₂<:FusionTree{I,N₂}} - T = scalartype(valtype(data)) - T ⊆ field(S) || @warn("scalartype(data) = $T ⊈ $(field(S)))", maxlog = 1) - return new{S,N₁,N₂,I,A,F₁,F₂}(data, codom, dom, rowr, colr) + E ⊆ field(S) || @warn("scalartype(data) = $T ⊈ $(field(S)))", maxlog = 1) + return new{E,S,N₁,N₂,I,A,F₁,F₂}(data, codom, dom, rowr, colr) end - function TensorMap{S,N₁,N₂,Trivial,A,Nothing,Nothing}(data::A, + function TensorMap{E,S,N₁,N₂,Trivial,A,Nothing,Nothing}(data::A, codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂}) where - {S<:IndexSpace,N₁,N₂,A<:DenseMatrix} - T = scalartype(data) - T ⊆ field(S) || - @warn("scalartype(data) = $T ⊈ $(field(S)))", maxlog = 1) - return new{S,N₁,N₂,Trivial,A,Nothing,Nothing}(data, codom, dom) + {E,S<:IndexSpace,N₁,N₂,A<:DenseMatrix{E}} + E ⊆ field(S) || @warn("scalartype(data) = $T ⊈ $(field(S)))", maxlog = 1) + return new{E,S,N₁,N₂,Trivial,A,Nothing,Nothing}(data, codom, dom) end end #! format: on """ - Tensor{S, N, I, A, F₁, F₂} = TensorMap{S, N, 0, I, A, F₁, F₂} + Tensor{E, S, N, I, A, F₁, F₂} = TensorMap{E, S, N, 0, I, A, F₁, F₂} Specific subtype of [`AbstractTensor`](@ref) for representing tensors whose data is stored in blocks of some subtype of `DenseMatrix`. -A `Tensor{S, N, I, A, F₁, F₂}` is actually a special case `TensorMap{S, N, 0, I, A, F₁, F₂}`, +A `Tensor{E, S, N, I, A, F₁, F₂}` is actually a special case `TensorMap{E, S, N, 0, I, A, F₁, F₂}`, i.e. a tensor map with only a non-trivial output space. """ -const Tensor{S,N,I,A,F₁,F₂} = TensorMap{S,N,0,I,A,F₁,F₂} +const Tensor{E,S,N,I,A,F₁,F₂} = TensorMap{E,S,N,0,I,A,F₁,F₂} + """ - TrivialTensorMap{S<:IndexSpace, N₁, N₂, A<:DenseMatrix} = TensorMap{S, N₁, N₂, Trivial, + TrivialTensorMap{E, S, N₁, N₂, A<:DenseMatrix} = TensorMap{E, S, N₁, N₂, Trivial, A, Nothing, Nothing} A special case of [`TensorMap`](@ref) for representing tensor maps with trivial symmetry, i.e., whose `sectortype` is `Trivial`. """ -const TrivialTensorMap{S,N₁,N₂,A<:DenseMatrix} = TensorMap{S,N₁,N₂,Trivial,A, - Nothing,Nothing} +const TrivialTensorMap{E,S,N₁,N₂,A<:DenseMatrix} = TensorMap{E,S,N₁,N₂,Trivial,A,Nothing, + Nothing} + """ - TrivialTensor{S, N, A} = TrivialTensorMap{S, N, 0, A} + TrivialTensor{E, S, N, A} = TrivialTensorMap{E, S, N, 0, A} A special case of [`Tensor`](@ref) for representing tensors with trivial symmetry, i.e., whose `sectortype` is `Trivial`. """ -const TrivialTensor{S,N,A} = TrivialTensorMap{S,N,0,A} +const TrivialTensor{E,S,N,A} = TrivialTensorMap{E,S,N,0,A} + +# TODO: check if argument order should change """ tensormaptype(::Type{S}, N₁::Int, N₂::Int, [::Type{T}]) where {S<:IndexSpace,T} -> ::Type{<:TensorMap} @@ -72,17 +74,19 @@ function tensormaptype(::Type{S}, N₁::Int, N₂::Int, ::Type{T}) where {S,T} I = sectortype(S) if T <: DenseMatrix M = T + E = scalartype(T) elseif T <: Number M = Matrix{T} + E = T else throw(ArgumentError("the final argument of `tensormaptype` should either be the scalar or the storage type, i.e. a subtype of `Number` or of `DenseMatrix`")) end if I === Trivial - return TensorMap{S,N₁,N₂,I,M,Nothing,Nothing} + return TensorMap{E,S,N₁,N₂,I,M,Nothing,Nothing} else F₁ = fusiontreetype(I, N₁) F₂ = fusiontreetype(I, N₂) - return TensorMap{S,N₁,N₂,I,SectorDict{I,M},F₁,F₂} + return TensorMap{E,S,N₁,N₂,I,SectorDict{I,M},F₁,F₂} end end tensormaptype(S, N₁, N₂=0) = tensormaptype(S, N₁, N₂, Float64) @@ -100,12 +104,12 @@ blocksectors(t::TensorMap) = keys(t.data) Return the type of the storage `A` of the tensor map. """ -function storagetype(::Type{<:TensorMap{<:IndexSpace,N₁,N₂,Trivial,A}}) where - {N₁,N₂,A<:DenseMatrix} +function storagetype(::Type{<:TrivialTensorMap{E,S,N₁,N₂,A}}) where + {E,S,N₁,N₂,A} return A end -function storagetype(::Type{<:TensorMap{<:IndexSpace,N₁,N₂,I,<:SectorDict{I,A}}}) where - {N₁,N₂,I<:Sector,A<:DenseMatrix} +function storagetype(::Type{<:TensorMap{E,S,N₁,N₂,I,<:SectorDict{I,A}}}) where + {E,S,N₁,N₂,I<:Sector,A<:DenseMatrix} return A end @@ -159,14 +163,15 @@ function TensorMap(data::AbstractDict{<:Sector,<:DenseMatrix}, codom::ProductSpa end F₁ = fusiontreetype(I, N₁) F₂ = fusiontreetype(I, N₂) + E = scalartype(valtype(data)) if !isreal(I) data2 = SectorDict(c => complex(data[c]) for c in blocksectoriterator) A = typeof(data2) - return TensorMap{S,N₁,N₂,I,A,F₁,F₂}(data2, codom, dom, rowr, colr) + return TensorMap{E,S,N₁,N₂,I,A,F₁,F₂}(data2, codom, dom, rowr, colr) else data2 = SectorDict(c => data[c] for c in blocksectoriterator) A = typeof(data2) - return TensorMap{S,N₁,N₂,I,A,F₁,F₂}(data2, codom, dom, rowr, colr) + return TensorMap{E,S,N₁,N₂,I,A,F₁,F₂}(data2, codom, dom, rowr, colr) end end @@ -179,7 +184,8 @@ function TensorMap(f, codom::ProductSpace{S,N₁}, d2 = dim(dom) data = f((d1, d2)) A = typeof(data) - return TensorMap{S,N₁,N₂,Trivial,A,Nothing,Nothing}(data, codom, dom) + E = scalartype(A) + return TensorMap{E,S,N₁,N₂,Trivial,A,Nothing,Nothing}(data, codom, dom) end blocksectoriterator = blocksectors(codom ← dom) rowr, rowdims = _buildblockstructure(codom, blocksectoriterator) @@ -193,7 +199,8 @@ function TensorMap(f, codom::ProductSpace{S,N₁}, F₁ = fusiontreetype(I, N₁) F₂ = fusiontreetype(I, N₂) A = typeof(data) - return TensorMap{S,N₁,N₂,I,A,F₁,F₂}(data, codom, dom, rowr, colr) + E = scalartype(valtype(A)) + return TensorMap{E,S,N₁,N₂,I,A,F₁,F₂}(data, codom, dom, rowr, colr) end # auxiliary function @@ -353,7 +360,8 @@ function TensorMap(data::DenseArray, codom::ProductSpace{S,N₁}, dom::ProductSp if sectortype(S) === Trivial data2 = reshape(data, (d1, d2)) A = typeof(data2) - return TensorMap{S,N₁,N₂,Trivial,A,Nothing,Nothing}(data2, codom, dom) + E = scalartype(A) + return TensorMap{E,S,N₁,N₂,Trivial,A,Nothing,Nothing}(data2, codom, dom) else t = TensorMap(zeros, eltype(data), codom, dom) ta = convert(Array, t) @@ -425,7 +433,7 @@ Base.similar(t::AbstractTensorMap, T::Type) = similar(t, T, space(t)) Base.similar(t::AbstractTensorMap) = similar(t, scalartype(t), space(t)) # actual implementation -function Base.similar(t::TensorMap{S}, ::Type{T}, P::TensorMapSpace{S}) where {T,S} +function Base.similar(t::TensorMap, ::Type{T}, P::TensorMapSpace{S}) where {T,S} N₁ = length(codomain(P)) N₂ = length(domain(P)) I = sectortype(S) @@ -433,14 +441,14 @@ function Base.similar(t::TensorMap{S}, ::Type{T}, P::TensorMapSpace{S}) where {T if I === Trivial data = similar(t.data, T, (dim(codomain(P)), dim(domain(P)))) A = typeof(data) - return TrivialTensorMap{S,N₁,N₂,A}(data, codomain(P), domain(P)) + return TrivialTensorMap{T,S,N₁,N₂,A}(data, codomain(P), domain(P)) end F₁ = fusiontreetype(I, N₁) F₂ = fusiontreetype(I, N₂) if space(t) == P data = SectorDict(c => similar(b, T) for (c, b) in blocks(t)) A = typeof(data) - return TensorMap{S,N₁,N₂,I,A,F₁,F₂}(data, codomain(P), domain(P), t.rowr, t.colr) + return TensorMap{T,S,N₁,N₂,I,A,F₁,F₂}(data, codomain(P), domain(P), t.rowr, t.colr) end blocksectoriterator = blocksectors(P) @@ -484,7 +492,7 @@ function Base.similar(t::TensorMap{S}, ::Type{T}, P::TensorMapSpace{S}) where {T data = SectorDict{I,M}(c => M(undef, (rowdims[c], coldims[c])) for c in blocksectoriterator) A = typeof(data) - return TensorMap{S,N₁,N₂,I,A,F₁,F₂}(data, codomain(P), domain(P), rowr, colr) + return TensorMap{T,S,N₁,N₂,I,A,F₁,F₂}(data, codomain(P), domain(P), rowr, colr) end function Base.complex(t::AbstractTensorMap) @@ -538,7 +546,7 @@ function block(t::TensorMap, s::Sector) end end -function blocks(t::TensorMap{<:IndexSpace,N₁,N₂,Trivial}) where {N₁,N₂} +function blocks(t::TrivialTensorMap) return SingletonDict(Trivial() => t.data) end blocks(t::TensorMap) = t.data @@ -547,7 +555,7 @@ fusiontrees(t::TrivialTensorMap) = ((nothing, nothing),) fusiontrees(t::TensorMap) = TensorKeyIterator(t.rowr, t.colr) """ - Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I}, + Base.getindex(t::TensorMap sectors::NTuple{N₁+N₂,I}) where {N₁,N₂,I<:Sector} -> StridedViews.StridedView t[sectors] @@ -561,8 +569,10 @@ respectively, then a `StridedViews.StridedView` of size This method is only available for the case where `FusionStyle(I) isa UniqueFusion`, since it assumes a uniquely defined coupled charge. """ -@inline function Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I}, - sectors::Tuple{Vararg{I}}) where {N₁,N₂,I<:Sector} +@inline function Base.getindex(t::TensorMap, sectors::Tuple{I,Vararg{I}}) where {I<:Sector} + I === sectortype(t) || throw(SectorMismatch("Not a valid sectortype for this tensor.")) + length(sectors) == numind(t) || + throw(ArgumentError("Number of sectors does not match.")) FusionStyle(I) isa UniqueFusion || throw(SectorMismatch("Indexing with sectors only possible if unique fusion")) s1 = TupleTools.getindices(sectors, codomainind(t)) @@ -584,9 +594,9 @@ end end """ - Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I}, + Base.getindex(t::TensorMap{E,S,N₁,N₂,I}, f₁::FusionTree{I,N₁}, - f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector} + f₂::FusionTree{I,N₂}) where {E,SN₁,N₂,I<:Sector} -> StridedViews.StridedView t[f₁, f₂] @@ -597,9 +607,9 @@ Return a view into the data slice of `t` corresponding to the splitting - fusion represents the slice of `block(t, c)` whose row indices correspond to `f₁.uncoupled` and column indices correspond to `f₂.uncoupled`. """ -@inline function Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I}, +@inline function Base.getindex(t::TensorMap{E,S,N₁,N₂,I}, f₁::FusionTree{I,N₁}, - f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector} + f₂::FusionTree{I,N₂}) where {E,S,N₁,N₂,I<:Sector} c = f₁.coupled @boundscheck begin c == f₂.coupled || throw(SectorMismatch()) @@ -613,10 +623,10 @@ column indices correspond to `f₂.uncoupled`. end """ - Base.setindex!(t::TensorMap{<:IndexSpace,N₁,N₂,I}, + Base.setindex!(t::TensorMap{E,S,N₁,N₂,I}, v, f₁::FusionTree{I,N₁}, - f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector} + f₂::FusionTree{I,N₂}) where {E,S,N₁,N₂,I<:Sector} t[f₁, f₂] = v Copies `v` into the data slice of `t` corresponding to the splitting - fusion tree pair @@ -624,12 +634,13 @@ Copies `v` into the data slice of `t` corresponding to the splitting - fusion t of size `(dims(codomain(t), f₁.uncoupled)..., dims(domain(t), f₂.uncoupled))` using `Base.copy!`. -See also [`Base.getindex(::TensorMap{<:IndexSpace,N₁,N₂,I<:Sector}, ::FusionTree{I<:Sector,N₁}, ::FusionTree{I<:Sector,N₂})`](@ref) +See also [`Base.getindex(::TensorMap{E,S,N₁,N₂,I<:Sector}, ::FusionTree{I<:Sector,N₁}, ::FusionTree{I<:Sector,N₂})`](@ref) """ -@propagate_inbounds function Base.setindex!(t::TensorMap{<:IndexSpace,N₁,N₂,I}, +@propagate_inbounds function Base.setindex!(t::TensorMap{E,S,N₁,N₂,I}, v, f₁::FusionTree{I,N₁}, - f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector} + f₂::FusionTree{I,N₂}) where {E,S,N₁,N₂, + I<:Sector} return copy!(getindex(t, f₁, f₂), v) end @@ -682,16 +693,16 @@ end function Base.summary(t::TensorMap) return print("TensorMap(", space(t), ")") end -function Base.show(io::IO, t::TensorMap{S}) where {S<:IndexSpace} +function Base.show(io::IO, t::TensorMap) if get(io, :compact, false) print(io, "TensorMap(", space(t), ")") return end println(io, "TensorMap(", space(t), "):") - if sectortype(S) == Trivial + if sectortype(t) == Trivial Base.print_array(io, t[]) println(io) - elseif FusionStyle(sectortype(S)) isa UniqueFusion + elseif FusionStyle(sectortype(t)) isa UniqueFusion for (f₁, f₂) in fusiontrees(t) println(io, "* Data for sector ", f₁.uncoupled, " ← ", f₂.uncoupled, ":") Base.print_array(io, t[f₁, f₂]) @@ -708,28 +719,28 @@ end # Real and imaginary parts #--------------------------- -function Base.real(t::AbstractTensorMap{S}) where {S} +function Base.real(t::AbstractTensorMap) # `isreal` for a `Sector` returns true iff the F and R symbols are real. This guarantees # that the real/imaginary part of a tensor `t` can be obtained by just taking # real/imaginary part of the degeneracy data. - if isreal(sectortype(S)) + if isreal(sectortype(t)) realdata = Dict(k => real(v) for (k, v) in blocks(t)) return TensorMap(realdata, codomain(t), domain(t)) else - msg = "`real` has not been implemented for `AbstractTensorMap{$(S)}`." + msg = "`real` has not been implemented for `$(typeof(t))`." throw(ArgumentError(msg)) end end -function Base.imag(t::AbstractTensorMap{S}) where {S} +function Base.imag(t::AbstractTensorMap) # `isreal` for a `Sector` returns true iff the F and R symbols are real. This guarantees # that the real/imaginary part of a tensor `t` can be obtained by just taking # real/imaginary part of the degeneracy data. - if isreal(sectortype(S)) + if isreal(sectortype(t)) imagdata = Dict(k => imag(v) for (k, v) in blocks(t)) return TensorMap(imagdata, codomain(t), domain(t)) else - msg = "`imag` has not been implemented for `AbstractTensorMap{$(S)}`." + msg = "`imag` has not been implemented for `$(typeof(t))`." throw(ArgumentError(msg)) end end @@ -741,19 +752,20 @@ function Base.convert(::Type{TensorMap}, t::AbstractTensorMap) return copy!(TensorMap(undef, scalartype(t), codomain(t), domain(t)), t) end -function Base.convert(T::Type{TensorMap{S,N₁,N₂,I,A,F₁,F₂}}, - t::AbstractTensorMap{S,N₁,N₂}) where {S,N₁,N₂,I,A,F₁,F₂} - if typeof(t) == T +function Base.convert(T::Type{<:TensorMap{E,S,N₁,N₂}}, + t::AbstractTensorMap{<:Any,S,N₁,N₂}) where {E,S,N₁,N₂} + if typeof(t) === T return t else - data = Dict{I,storagetype(T)}(c => convert(storagetype(T), b) - for (c, b) in blocks(t)) + data = Dict{sectortype(T),storagetype(T)}(c => convert(storagetype(T), b) + for (c, b) in blocks(t)) return TensorMap(data, codomain(t), domain(t)) end end -function Base.promote_rule(::Type{<:T1}, - t2::Type{<:T2}) where {S,N₁,N₂,T1<:TensorMap{S,N₁,N₂}, - T2<:TensorMap{S,N₁,N₂}} - return tensormaptype(S, N₁, N₂, promote_type(storagetype(T1), storagetype(T2))) +function Base.promote_rule(::Type{<:T₁}, + ::Type{<:T₂}) where {S,N₁,N₂, + T₁<:TensorMap{<:Any,S,N₁,N₂}, + T₂<:TensorMap{<:Any,S,N₁,N₂}} + return tensormaptype(S, N₁, N₂, promote_type(storagetype(T₁), storagetype(T₂))) end diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl index 9f8ed5a9..9114f767 100644 --- a/src/tensors/tensoroperations.jl +++ b/src/tensors/tensoroperations.jl @@ -34,9 +34,9 @@ function _canonicalize(p::IndexTuple, t::AbstractTensorMap) end # tensoradd! -function TO.tensoradd!(C::AbstractTensorMap{S}, pC::Index2Tuple, - A::AbstractTensorMap{S}, conjA::Symbol, - α::Number, β::Number, backend::Backend...) where {S} +function TO.tensoradd!(C::AbstractTensorMap, pC::Index2Tuple, + A::AbstractTensorMap, conjA::Symbol, + α::Number, β::Number, backend::Backend...) if conjA == :N A′ = A pC′ = _canonicalize(pC, C) @@ -50,17 +50,17 @@ function TO.tensoradd!(C::AbstractTensorMap{S}, pC::Index2Tuple, return C end -function TO.tensoradd_type(TC, ::Index2Tuple{N₁,N₂}, A::AbstractTensorMap{S}, - ::Symbol) where {S,N₁,N₂} +function TO.tensoradd_type(TC, ::Index2Tuple{N₁,N₂}, A::AbstractTensorMap, + ::Symbol) where {N₁,N₂} M = similarstoragetype(A, TC) - return tensormaptype(S, N₁, N₂, M) + return tensormaptype(spacetype(A), N₁, N₂, M) end function TO.tensoradd_structure(pC::Index2Tuple{N₁,N₂}, - A::AbstractTensorMap{S}, conjA::Symbol) where {S,N₁,N₂} + A::AbstractTensorMap, conjA::Symbol) where {N₁,N₂} if conjA == :N - cod = ProductSpace{S,N₁}(space.(Ref(A), pC[1])) - dom = ProductSpace{S,N₂}(dual.(space.(Ref(A), pC[2]))) + cod = ProductSpace{spacetype(A),N₁}(space.(Ref(A), pC[1])) + dom = ProductSpace{spacetype(A),N₂}(dual.(space.(Ref(A), pC[2]))) return dom → cod else return TO.tensoradd_structure(adjointtensorindices(A, pC), adjoint(A), :N) @@ -68,9 +68,9 @@ function TO.tensoradd_structure(pC::Index2Tuple{N₁,N₂}, end # tensortrace! -function TO.tensortrace!(C::AbstractTensorMap{S}, p::Index2Tuple, - A::AbstractTensorMap{S}, q::Index2Tuple, conjA::Symbol, - α::Number, β::Number, backend::Backend...) where {S} +function TO.tensortrace!(C::AbstractTensorMap, p::Index2Tuple, + A::AbstractTensorMap, q::Index2Tuple, conjA::Symbol, + α::Number, β::Number, backend::Backend...) if conjA == :N A′ = A p′ = _canonicalize(p, C) @@ -89,10 +89,10 @@ function TO.tensortrace!(C::AbstractTensorMap{S}, p::Index2Tuple, end # tensorcontract! -function TO.tensorcontract!(C::AbstractTensorMap{S,N₁,N₂}, pAB::Index2Tuple, - A::AbstractTensorMap{S}, pA::Index2Tuple, conjA::Symbol, - B::AbstractTensorMap{S}, pB::Index2Tuple, conjB::Symbol, - α::Number, β::Number, backend::Backend...) where {S,N₁,N₂} +function TO.tensorcontract!(C::AbstractTensorMap, pAB::Index2Tuple, + A::AbstractTensorMap, pA::Index2Tuple, conjA::Symbol, + B::AbstractTensorMap, pB::Index2Tuple, conjB::Symbol, + α::Number, β::Number, backend::Backend...) pAB′ = _canonicalize(pAB, C) if conjA == :N A′ = A @@ -117,28 +117,30 @@ function TO.tensorcontract!(C::AbstractTensorMap{S,N₁,N₂}, pAB::Index2Tuple, end function TO.tensorcontract_type(TC, ::Index2Tuple{N₁,N₂}, - A::AbstractTensorMap{S}, pA, conjA, - B::AbstractTensorMap{S}, pB, conjB) where {S,N₁,N₂} + A::AbstractTensorMap, pA, conjA, + B::AbstractTensorMap, pB, conjB) where {N₁,N₂} M = similarstoragetype(A, TC) M == similarstoragetype(B, TC) || throw(ArgumentError("incompatible storage types")) - return tensormaptype(S, N₁, N₂, M) + spacetype(A) == spacetype(B) || throw(SpaceMismatch("incompatible space types")) + return tensormaptype(spacetype(A), N₁, N₂, M) end function TO.tensorcontract_structure(pC::Index2Tuple{N₁,N₂}, - A::AbstractTensorMap{S}, pA::Index2Tuple, conjA, - B::AbstractTensorMap{S}, pB::Index2Tuple, - conjB) where {S,N₁,N₂} + A::AbstractTensorMap, pA::Index2Tuple, conjA, + B::AbstractTensorMap, pB::Index2Tuple, + conjB) where {N₁,N₂} + spacetype(A) == spacetype(B) || throw(SpaceMismatch("incompatible space types")) spaces1 = TO.flag2op(conjA).(space.(Ref(A), pA[1])) spaces2 = TO.flag2op(conjB).(space.(Ref(B), pB[2])) spaces = (spaces1..., spaces2...) - cod = ProductSpace{S,N₁}(getindex.(Ref(spaces), pC[1])) - dom = ProductSpace{S,N₂}(dual.(getindex.(Ref(spaces), pC[2]))) + cod = ProductSpace{spacetype(A),N₁}(getindex.(Ref(spaces), pC[1])) + dom = ProductSpace{spacetype(A),N₂}(dual.(getindex.(Ref(spaces), pC[2]))) return dom → cod end -function TO.checkcontractible(tA::AbstractTensorMap{S}, iA::Int, conjA::Symbol, - tB::AbstractTensorMap{S}, iB::Int, conjB::Symbol, - label) where {S} +function TO.checkcontractible(tA::AbstractTensorMap, iA::Int, conjA::Symbol, + tB::AbstractTensorMap, iB::Int, conjB::Symbol, + label) sA = TO.tensorstructure(tA, iA, conjA)' sB = TO.tensorstructure(tB, iB, conjB) sA == sB || @@ -154,17 +156,27 @@ TO.tensorcost(t::AbstractTensorMap, i::Int) = dim(space(t, i)) # Trace implementation #---------------------- -function trace_permute!(tdst::AbstractTensorMap{S,N₁,N₂}, - tsrc::AbstractTensorMap{S}, - (p₁, p₂)::Index2Tuple{N₁,N₂}, - (q₁, q₂)::Index2Tuple{N₃,N₃}, +function trace_permute!(tdst::AbstractTensorMap, + tsrc::AbstractTensorMap, + (p₁, p₂)::Index2Tuple, + (q₁, q₂)::Index2Tuple, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂,N₃} + backend::Backend...) + # some input checks + (S = spacetype(tdst)) == spacetype(tsrc) || + throw(SpaceMismatch("incompatible spacetypes")) if !(BraidingStyle(sectortype(S)) isa SymmetricBraiding) throw(SectorMismatch("only tensors with symmetric braiding rules can be contracted; try `@planar` instead")) end + (N₃ = length(q₁)) == length(q₂) || + throw(IndexError("number of trace indices does not match")) + + N₁, N₂ = length(p₁), length(p₂) + @boundscheck begin + numout(tdst) == N₁ || throw(IndexError("number of output indices does not match")) + numin(tdst) == N₂ || throw(IndexError("number of input indices does not match")) all(i -> space(tsrc, p₁[i]) == space(tdst, i), 1:N₁) || throw(SpaceMismatch("trace: tsrc = $(codomain(tsrc))←$(domain(tsrc)), tdst = $(codomain(tdst))←$(domain(tdst)), p₁ = $(p₁), p₂ = $(p₂)")) @@ -217,15 +229,15 @@ end # TODO: contraction with either A or B a rank (1, 1) tensor does not require to # permute the fusion tree and should therefore be special cased. This will speed # up MPS algorithms -function contract!(C::AbstractTensorMap{S}, - A::AbstractTensorMap{S}, - (oindA, cindA)::Index2Tuple{N₁,N₃}, - B::AbstractTensorMap{S}, - (cindB, oindB)::Index2Tuple{N₃,N₂}, +function contract!(C::AbstractTensorMap, + A::AbstractTensorMap, (oindA, cindA)::Index2Tuple, + B::AbstractTensorMap, (cindB, oindB)::Index2Tuple, (p₁, p₂)::Index2Tuple, - α::Number, - β::Number, - backend::Backend...) where {S,N₁,N₂,N₃} + α::Number, β::Number, + backend::Backend...) + length(cindA) == length(cindB) || + throw(IndexError("number of contracted indices does not match")) + N₁, N₂ = length(oindA), length(oindB) # find optimal contraction scheme hsp = has_shared_permute @@ -275,16 +287,17 @@ function contract!(C::AbstractTensorMap{S}, end # TODO: also transform _contract! into new interface, and add backend support -function _contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S}, - β, C::AbstractTensorMap{S}, - oindA::IndexTuple{N₁}, cindA::IndexTuple, - oindB::IndexTuple{N₂}, cindB::IndexTuple, - p₁::IndexTuple, p₂::IndexTuple) where {S,N₁,N₂} - if !(BraidingStyle(sectortype(S)) isa SymmetricBraiding) +function _contract!(α, A::AbstractTensorMap, B::AbstractTensorMap, + β, C::AbstractTensorMap, + oindA::IndexTuple, cindA::IndexTuple, + oindB::IndexTuple, cindB::IndexTuple, + p₁::IndexTuple, p₂::IndexTuple) + if !(BraidingStyle(sectortype(C)) isa SymmetricBraiding) throw(SectorMismatch("only tensors with symmetric braiding rules can be contracted; try `@planar` instead")) end + N₁, N₂ = length(oindA), length(oindB) copyA = false - if BraidingStyle(sectortype(S)) isa Fermionic + if BraidingStyle(sectortype(A)) isa Fermionic for i in cindA if !isdual(space(A, i)) copyA = true @@ -293,7 +306,7 @@ function _contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S}, end A′ = permute(A, (oindA, cindA); copy=copyA) B′ = permute(B, (cindB, oindB)) - if BraidingStyle(sectortype(S)) isa Fermionic + if BraidingStyle(sectortype(A)) isa Fermionic for i in domainind(A′) if !isdual(space(A′, i)) A′ = twist!(A′, i) @@ -315,7 +328,7 @@ end # Scalar implementation #----------------------- -function scalar(t::AbstractTensorMap{S}) where {S<:IndexSpace} +function scalar(t::AbstractTensorMap) return dim(codomain(t)) == dim(domain(t)) == 1 ? first(blocks(t))[2][1, 1] : throw(DimensionMismatch()) end diff --git a/test/braidingtensor.jl b/test/braidingtensor.jl index b61abc61..726a28c4 100644 --- a/test/braidingtensor.jl +++ b/test/braidingtensor.jl @@ -13,119 +13,119 @@ for V in (V1, V2, V3) t = TensorMap(randn, V * V' * V' * V, V * V') - ττ = copy(BraidingTensor(V, V')) + ττ = TensorMap(BraidingTensor(V, V')) @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-1 -2; 1 2] * t[1 2 -3 -4; -5 -6] @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-1 -2; 1 2] * t[1 2 -3 -4; -5 -6] @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -2 -1] * t[1 2 -3 -4; -5 -6] @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-2 2; -1 1] * t[1 2 -3 -4; -5 -6] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V', V')) + ττ = TensorMap(BraidingTensor(V', V')) @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6] @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6] @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -3 -2] * t[-1 1 2 -4; -5 -6] @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-3 2; -2 1] * t[-1 1 2 -4; -5 -6] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V', V)) + ττ = TensorMap(BraidingTensor(V', V)) @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6] @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6] @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[-3 2; -2 1] * t[-1 1 2 -4; -5 -6] @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V, V')) + ττ = TensorMap(BraidingTensor(V, V')) @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-3 2; -2 1] * t[-1 1 2 -4; -5 -6] @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-3 2; -2 1] * t[-1 1 2 -4; -5 -6] @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6] @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[2 1; -3 -2] * t[-1 1 2 -4; -5 -6] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V, V)) + ττ = TensorMap(BraidingTensor(V, V)) @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[2 1; -3 -2] * t[-1 1 2 -4; -5 -6] @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[2 1; -3 -2] * t[-1 1 2 -4; -5 -6] @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6] @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V', V)) + ττ = TensorMap(BraidingTensor(V', V)) @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-3 -4; 1 2] * t[-1 -2 1 2; -5 -6] @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-3 -4; 1 2] * t[-1 -2 1 2; -5 -6] @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -4 -3] * t[-1 -2 1 2; -5 -6] @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-4 2; -3 1] * t[-1 -2 1 2; -5 -6] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V', V)) + ττ = TensorMap(BraidingTensor(V', V)) @planar2 t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ[1 2; -5 -6] @planar2 t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * ττ[1 2; -5 -6] @planar2 t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ[-6 -5; 2 1] @planar2 t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ'[2 -6; 1 -5] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V, V')) + ττ = TensorMap(BraidingTensor(V, V')) @planar2 t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ'[1 2; -5 -6] @planar2 t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * ττ'[1 2; -5 -6] @planar2 t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ'[-6 -5; 2 1] @planar2 t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ[2 -6; 1 -5] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V, V)) + ττ = TensorMap(BraidingTensor(V, V)) @planar2 t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * τ[-4 -6; 1 2] @planar2 t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * ττ[-4 -6; 1 2] @planar2 t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * τ[2 1; -6 -4] @planar2 t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * τ'[-6 2; -4 1] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V', V)) + ττ = TensorMap(BraidingTensor(V', V)) @planar2 t1[(); (-1, -2)] := τ[2 1; 3 4] * t[1 2 3 4; -1 -2] @planar2 t2[(); (-1, -2)] := ττ[2 1; 3 4] * t[1 2 3 4; -1 -2] @planar2 t3[(); (-1, -2)] := τ[4 3; 1 2] * t[1 2 3 4; -1 -2] @planar2 t4[(); (-1, -2)] := τ'[1 4; 2 3] * t[1 2 3 4; -1 -2] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V, V)) + ττ = TensorMap(BraidingTensor(V, V)) @planar2 t1[-1; -2] := τ[2 1; 3 4] * t[-1 1 2 3; -2 4] @planar2 t2[-1; -2] := ττ[2 1; 3 4] * t[-1 1 2 3; -2 4] @planar2 t3[-1; -2] := τ[4 3; 1 2] * t[-1 1 2 3; -2 4] @planar2 t4[-1; -2] := τ'[1 4; 2 3] * t[-1 1 2 3; -2 4] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V, V')) + ττ = TensorMap(BraidingTensor(V, V')) @planar2 t1[-1 -2] := τ[2 1; 3 4] * t[-1 -2 1 2; 4 3] @planar2 t2[-1 -2] := ττ[2 1; 3 4] * t[-1 -2 1 2; 4 3] @planar2 t3[-1 -2] := τ[4 3; 1 2] * t[-1 -2 1 2; 4 3] @planar2 t4[-1 -2] := τ'[1 4; 2 3] * t[-1 -2 1 2; 4 3] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V, V')) + ττ = TensorMap(BraidingTensor(V, V')) @planar2 t1[-1 -2; -3 -4] := τ[-1 3; 1 2] * t[1 2 3 -2; -3 -4] @planar2 t2[-1 -2; -3 -4] := ττ[-1 3; 1 2] * t[1 2 3 -2; -3 -4] @planar2 t3[-1 -2; -3 -4] := τ[2 1; 3 -1] * t[1 2 3 -2; -3 -4] @planar2 t4[-1 -2; -3 -4] := τ'[3 2; -1 1] * t[1 2 3 -2; -3 -4] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V', V')) + ττ = TensorMap(BraidingTensor(V', V')) @planar2 t1[-1 -2; -3 -4] := τ'[-2 3; 1 2] * t[-1 1 2 3; -3 -4] @planar2 t2[-1 -2; -3 -4] := ττ'[-2 3; 1 2] * t[-1 1 2 3; -3 -4] @planar2 t3[-1 -2; -3 -4] := τ'[2 1; 3 -2] * t[-1 1 2 3; -3 -4] @planar2 t4[-1 -2; -3 -4] := τ[3 2; -2 1] * t[-1 1 2 3; -3 -4] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V', V)) + ττ = TensorMap(BraidingTensor(V', V)) @planar2 t1[-1 -2 -3; -4] := τ[-3 3; 1 2] * t[-1 -2 1 2; -4 3] @planar2 t2[-1 -2 -3; -4] := ττ[-3 3; 1 2] * t[-1 -2 1 2; -4 3] @planar2 t3[-1 -2 -3; -4] := τ[2 1; 3 -3] * t[-1 -2 1 2; -4 3] @planar2 t4[-1 -2 -3; -4] := τ'[3 2; -3 1] * t[-1 -2 1 2; -4 3] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V', V)) + ττ = TensorMap(BraidingTensor(V', V)) @planar2 t1[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ[1 2; -4 3] @planar2 t2[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * ττ[1 2; -4 3] @planar2 t3[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ[3 -4; 2 1] @planar2 t4[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ'[2 3; 1 -4] @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4) - ττ = copy(BraidingTensor(V, V')) + ττ = TensorMap(BraidingTensor(V, V')) @planar2 t1[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ'[1 2; -4 3] @planar2 t2[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * ττ'[1 2; -4 3] @planar2 t3[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ'[3 -4; 2 1] diff --git a/test/planar.jl b/test/planar.jl index d7edf146..c9b8cc54 100644 --- a/test/planar.jl +++ b/test/planar.jl @@ -12,13 +12,13 @@ function force_planar(V::GradedSpace) return GradedSpace((c ⊠ PlanarTrivial() => dim(V, c) for c in sectors(V))..., isdual(V)) end force_planar(V::ProductSpace) = mapreduce(force_planar, ⊗, V) -function force_planar(tsrc::TensorMap{ComplexSpace}) +function force_planar(tsrc::TensorMap{<:Any,ComplexSpace}) tdst = TensorMap(undef, scalartype(tsrc), force_planar(codomain(tsrc)) ← force_planar(domain(tsrc))) copyto!(blocks(tdst)[PlanarTrivial()], blocks(tsrc)[Trivial()]) return tdst end -function force_planar(tsrc::TensorMap{<:GradedSpace}) +function force_planar(tsrc::TensorMap{<:Any,<:GradedSpace}) tdst = TensorMap(undef, scalartype(tsrc), force_planar(codomain(tsrc)) ← force_planar(domain(tsrc))) for (c, b) in blocks(tsrc) diff --git a/test/runtests.jl b/test/runtests.jl index 21760665..34037d27 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,9 +10,6 @@ using Base.Iterators: take, product # const SU3Irrep = SUNIrrep{3} using LinearAlgebra: LinearAlgebra -include("newsectors.jl") -using .NewSectors - const TK = TensorKit Random.seed!(1234) @@ -49,15 +46,14 @@ function hasfusiontensor(I::Type{<:Sector}) end end -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, +sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, Z3Irrep ⊠ Z4Irrep, + U1Irrep, CU1Irrep, SU2Irrep, + FermionParity, FermionParity ⊠ FermionParity, + FermionParity ⊠ U1Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, # Hubbard-like + FibonacciAnyon, IsingAnyon, Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon) Ti = time() -include("sectors.jl") include("fusiontrees.jl") include("spaces.jl") include("tensors.jl") diff --git a/test/sectors.jl b/test/sectors.jl deleted file mode 100644 index a0702552..00000000 --- a/test/sectors.jl +++ /dev/null @@ -1,134 +0,0 @@ -println("------------------------------------") -println("Sectors") -println("------------------------------------") -@timedtestset "Sectors" verbose = true begin - @timedtestset "Sector properties of $(TensorKit.type_repr(I))" for I in sectorlist - Istr = TensorKit.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(TensorKit.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] - @test eltype(it) === I - @test collect(it) isa Array{I} - @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 - if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector - @test_throws ArgumentError values(I)[i] - @test_throws ArgumentError TensorKit.findindex(values(I), s) - elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int}) - @test s == @constinferred (values(I)[i]) - @test TensorKit.findindex(values(I), s) == i - end - sprev = s - i >= 10 && break - end - @test one(I) == first(values(I)) - if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector - @test_throws ArgumentError TensorKit.findindex(values(I), one(I)) - elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int}) - @test (@constinferred TensorKit.findindex(values(I), one(I))) == 1 - for s in smallset(I) - @test (@constinferred values(I)[TensorKit.findindex(values(I), s)]) == s - end - end - end - if BraidingStyle(I) isa Bosonic && hasfusiontensor(I) - @testset "Sector $I: 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 - 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 -end