diff --git a/.github/workflows/CI-Sectors.yml b/.github/workflows/CI-Sectors.yml new file mode 100644 index 00000000..1d3dc754 --- /dev/null +++ b/.github/workflows/CI-Sectors.yml @@ -0,0 +1,53 @@ +name: CI - TensorKitSectors +on: + push: + branches: + - 'master' + - 'main' + - 'release-' + paths: + - TensorKitSectors/ + tags: '*' + pull_request: + paths: + - TensorKitSectors/ + workflow_dispatch: + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + test: + name: TensorKitSectors - ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1' # automatically expands to the latest stable 1.x release of Julia + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@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/.github/workflows/CI.yml b/.github/workflows/CI.yml index f56ac38f..9959b2f9 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -21,7 +21,7 @@ jobs: fail-fast: false matrix: version: - - '1.6' # LTS version + - '1.8' # LTS version - '1' # automatically expands to the latest stable 1.x release of Julia os: - ubuntu-latest 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..bc30dd82 --- /dev/null +++ b/Changelog.md @@ -0,0 +1,44 @@ +# Planned changes for v1.0 + +Features that are planned to be implemented before the release of v1.0.0, in no particular order. + +- [x] 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 +- [x] 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. + +### Refactor tensormap constructors + +This PR refactors the constructors for `TensorMap` to be more in line with `Array` +constructors in Julia. In particular, this means that the default way to create +uninitialized tensors is now `TensorMap{E}(undef, codomain ← domain)`, reminiscent of +`Array{E}(undef, dims)`. Several convenience constructors are also added: `ones`, `zeros`, +`rand` and `randn` construct tensors when `dims` is replaced by `domain ← codomain`. + +### TensorOperations v5 + +This PR bumps the compatibility of `TensorOperations` to v5. This is a breaking change +as there are some changes in the API. diff --git a/Project.toml b/Project.toml index 8ac0760e..4990f8fd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,14 @@ name = "TensorKit" uuid = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" authors = ["Jutho Haegeman"] -version = "0.12.5" +version = "1.0.0-DEV" [deps] HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" @@ -34,7 +35,7 @@ LinearAlgebra = "1" PackageExtensionCompat = "1" Random = "1" Strided = "2" -TensorOperations = "4.1" +TensorOperations = "5" Test = "1" TestExtras = "0.2" TupleTools = "1.1" @@ -50,11 +51,10 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a" WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" [targets] -test = ["Aqua", "Combinatorics", "HalfIntegers", "LinearAlgebra", "Random", "TensorOperations", "Test", "TestExtras", "WignerSymbols", "ChainRulesCore", "ChainRulesTestUtils", "FiniteDifferences"] +test = ["Aqua", "Combinatorics", "HalfIntegers", "LinearAlgebra", "TensorOperations", "Test", "TestExtras", "WignerSymbols", "ChainRulesCore", "ChainRulesTestUtils", "FiniteDifferences"] 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..89d91339 --- /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, invsqrtdim, frobeniusschur, twist, fusiontensor, dual +export otimes, deligneproduct, times +export FusionStyle, UniqueFusion, MultipleFusion, SimpleFusion, GenericFusion, + MultiplicityFreeFusion +export BraidingStyle, NoBraiding, SymmetricBraiding, Bosonic, Fermionic, Anyonic +export SectorSet, SectorValues, findindex, vertex_ind2label, vertex_labeltype + +export pentagon_equation, hexagon_equation + +export Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep +export ProductSector +export FermionParity, FermionNumber, FermionSpin +export PlanarTrivial, FibonacciAnyon, IsingAnyon + +# unicode exports +# --------------- +export ⊠, ⊗, × +export ℤ, ℤ₂, ℤ₃, ℤ₄, U₁, SU, SU₂, CU₁ +export fℤ₂, fU₁, fSU₂ + +# imports +# ------- +using Base: SizeUnknown, HasLength, IsInfinite +using Base: HasEltype, EltypeUnknown +using Base.Iterators: product, filter +using Base: tuple_type_head, tuple_type_tail + +using LinearAlgebra: tr +using TensorOperations +using HalfIntegers +using WignerSymbols + +# includes +# -------- +include("auxiliary.jl") +include("sectors.jl") +include("trivial.jl") +include("groups.jl") +include("irreps.jl") # irreps of symmetry groups, with bosonic braiding +include("product.jl") # direct product of different sectors +include("fermions.jl") # irreps with defined fermionparity and fermionic braiding +include("anyons.jl") # non-group sectors + +# precompile +# ---------- +include("precompile.jl") + +function __precompile__() + for I in (Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep, + FermionParity, FermionNumber, FermionSpin, PlanarTrivial, FibonacciAnyon, + IsingAnyon) + precompile_sector(I) + end +end + +end # module TensorKitSectors diff --git a/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..9e216017 --- /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(invsqrtdim, (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 86% rename from src/sectors/sectors.jl rename to TensorKitSectors/src/sectors.jl index bfb83312..51425285 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,51 +128,40 @@ 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 -abstract type MultipleFusion <: FusionStyle end -struct SimpleFusion <: MultipleFusion # multiple fusion but multiplicity free -end -struct GenericFusion <: MultipleFusion # multiple fusion with multiplicities -end -const MultiplicityFreeFusion = Union{UniqueFusion,SimpleFusion} - """ - FusionStyle(a::Sector) -> ::FusionStyle - FusionStyle(I::Type{<:Sector}) -> ::FusionStyle + FusionStyle(::Sector) + FusionStyle(I::Type{<:Sector}) -Return the type of fusion behavior of sectors of type I, which can be either +Trait to describe the fusion behavior of sectors of type `I`, which can be either * `UniqueFusion()`: single fusion output when fusing two sectors; * `SimpleFusion()`: multiple outputs, but every output occurs at most one, - also known as multiplicity free (e.g. irreps of ``SU(2)``); + also known as multiplicity-free (e.g. irreps of ``SU(2)``); * `GenericFusion()`: multiple outputs that can occur more than once (e.g. irreps of ``SU(3)``). + There is an abstract supertype `MultipleFusion` of which both `SimpleFusion` and -`GenericFusion` are subtypes. Furthermore, there is a type alias `MultiplicityFreeFusion` for those fusion types which do not require muliplicity labels, i.e. +`GenericFusion` are subtypes. Furthermore, there is a type alias `MultiplicityFreeFusion` +for those fusion types which do not require muliplicity labels, i.e. `MultiplicityFreeFusion = Union{UniqueFusion,SimpleFusion}`. """ +abstract type FusionStyle end FusionStyle(a::Sector) = FusionStyle(typeof(a)) -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 +struct UniqueFusion <: FusionStyle end # unique fusion output when fusing two sectors +abstract type MultipleFusion <: FusionStyle end +struct SimpleFusion <: MultipleFusion end # multiple fusion but multiplicity free +struct GenericFusion <: MultipleFusion end # multiple fusion with multiplicities +const MultiplicityFreeFusion = Union{UniqueFusion,SimpleFusion} + +# combine fusion properties of tensor products of sectors +Base.:&(f::F, ::F) where {F<:FusionStyle} = f +Base.:&(f₁::FusionStyle, f₂::FusionStyle) = f₂ & f₁ + +Base.:&(::SimpleFusion, ::UniqueFusion) = SimpleFusion() +Base.:&(::GenericFusion, ::UniqueFusion) = GenericFusion() +Base.:&(::GenericFusion, ::SimpleFusion) = GenericFusion() """ Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I<:Sector} @@ -204,23 +183,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} - -Returns the R-symbol ``R^{ab}_c`` that maps between ``c → a ⊗ b`` and ``c → b ⊗ a`` as in -``` -a -<-μ-<- c b -<-ν-<- c - ∨ -> Rsymbol(a,b,c)[μ,ν] v - b a -``` -If `FusionStyle(I)` is `UniqueFusion()` or `SimpleFusion()`, the R-symbol is a -number. Otherwise it is a square matrix with row and column size -`Nsymbol(a,b,c) == Nsymbol(b,a,c)`. -""" -function Rsymbol end -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 @@ -247,14 +209,6 @@ Return the type of labels for the fusion vertices of sectors of type `I`. """ vertex_labeltype(I::Type{<:Sector}) = typeof(vertex_ind2label(1, one(I), one(I), one(I))) -# combine fusion properties of tensor products of sectors -Base.:&(f::F, ::F) where {F<:FusionStyle} = f -Base.:&(f₁::FusionStyle, f₂::FusionStyle) = f₂ & f₁ - -Base.:&(::SimpleFusion, ::UniqueFusion) = SimpleFusion() -Base.:&(::GenericFusion, ::UniqueFusion) = GenericFusion() -Base.:&(::GenericFusion, ::SimpleFusion) = GenericFusion() - # properties that can be determined in terms of the F symbol # TODO: find mechanism for returning these numbers with custom type T<:AbstractFloat """ @@ -272,7 +226,7 @@ function dim(a::Sector) end end sqrtdim(a::Sector) = (FusionStyle(a) isa UniqueFusion) ? 1 : sqrt(dim(a)) -isqrtdim(a::Sector) = (FusionStyle(a) isa UniqueFusion) ? 1 : inv(sqrt(dim(a))) +invsqrtdim(a::Sector) = (FusionStyle(a) isa UniqueFusion) ? 1 : inv(sqrt(dim(a))) """ frobeniusschur(a::Sector) @@ -287,12 +241,17 @@ function frobeniusschur(a::Sector) end end -""" - twist(a::Sector) - -Return the twist of a sector `a` -""" -twist(a::Sector) = sum(dim(b) / dim(a) * tr(Rsymbol(a, a, b)) for b in a ⊗ a) +# Not necessary +function Asymbol(a::I, b::I, c::I) where {I<:Sector} + if FusionStyle(I) isa UniqueFusion || FusionStyle(I) isa SimpleFusion + (sqrtdim(a) * sqrtdim(b) * invsqrtdim(c)) * + conj(frobeniusschur(a) * Fsymbol(dual(a), a, b, b, one(a), c)) + else + reshape((sqrtdim(a) * sqrtdim(b) * invsqrtdim(c)) * + conj(frobeniusschur(a) * Fsymbol(dual(a), a, b, b, one(a), c)), + (Nsymbol(a, b, c), Nsymbol(dual(a), c, b))) + end +end """ Bsymbol(a::I, b::I, c::I) where {I<:Sector} @@ -310,31 +269,35 @@ number. Otherwise it is a square matrix with row and column size """ function Bsymbol(a::I, b::I, c::I) where {I<:Sector} if FusionStyle(I) isa UniqueFusion || FusionStyle(I) isa SimpleFusion - (sqrtdim(a) * sqrtdim(b) * isqrtdim(c)) * Fsymbol(a, b, dual(b), a, c, one(a)) + (sqrtdim(a) * sqrtdim(b) * invsqrtdim(c)) * Fsymbol(a, b, dual(b), a, c, one(a)) else - reshape((sqrtdim(a) * sqrtdim(b) * isqrtdim(c)) * + reshape((sqrtdim(a) * sqrtdim(b) * invsqrtdim(c)) * Fsymbol(a, b, dual(b), a, c, one(a)), (Nsymbol(a, b, c), Nsymbol(c, dual(b), a))) end end -# Not necessary -function Asymbol(a::I, b::I, c::I) where {I<:Sector} - if FusionStyle(I) isa UniqueFusion || FusionStyle(I) isa SimpleFusion - (sqrtdim(a) * sqrtdim(b) * isqrtdim(c)) * - conj(frobeniusschur(a) * Fsymbol(dual(a), a, b, b, one(a), c)) - else - reshape((sqrtdim(a) * sqrtdim(b) * isqrtdim(c)) * - conj(frobeniusschur(a) * Fsymbol(dual(a), a, b, b, one(a), c)), - (Nsymbol(a, b, c), Nsymbol(dual(a), c, b))) - end -end - # Braiding: #------------------------------------------------- # trait to describe type to denote how the elementary spaces in a tensor product space # interact under permutations or actions of the braid group -abstract type BraidingStyle end # generic braiding +""" + BraidingStyle(::Sector) -> ::BraidingStyle + BraidingStyle(I::Type{<:Sector}) -> ::BraidingStyle + +Return the type of braiding and twist behavior of sectors of type `I`, which can be either +* `Bosonic()`: symmetric braiding with trivial twist (i.e. identity) +* `Fermionic()`: symmetric braiding with non-trivial twist (squares to identity) +* `Anyonic()`: general ``R_(a,b)^c`` phase or matrix (depending on `SimpleFusion` or + `GenericFusion` fusion) and arbitrary twists + +Note that `Bosonic` and `Fermionic` are subtypes of `SymmetricBraiding`, which means that +braids are in fact equivalent to crossings (i.e. braiding twice is an identity: +`isone(Rsymbol(b,a,c)*Rsymbol(a,b,c)) == true`) and permutations are uniquely defined. +""" +abstract type BraidingStyle end +BraidingStyle(a::Sector) = BraidingStyle(typeof(a)) + abstract type HasBraiding <: BraidingStyle end struct NoBraiding <: BraidingStyle end abstract type SymmetricBraiding <: HasBraiding end # symmetric braiding => actions of permutation group are well defined @@ -352,21 +315,28 @@ Base.:&(::Fermionic, ::NoBraiding) = NoBraiding() Base.:&(::Anyonic, ::NoBraiding) = NoBraiding() """ - BraidingStyle(::Sector) -> ::BraidingStyle - BraidingStyle(I::Type{<:Sector}) -> ::BraidingStyle + Rsymbol(a::I, b::I, c::I) where {I<:Sector} -Return the type of braiding and twist behavior of sectors of type `I`, which can be either -* `Bosonic()`: symmetric braiding with trivial twist (i.e. identity) -* `Fermionic()`: symmetric braiding with non-trivial twist (squares to identity) -* `Anyonic()`: general ``R_(a,b)^c`` phase or matrix (depending on `SimpleFusion` or - `GenericFusion` fusion) and arbitrary twists +Returns the R-symbol ``R^{ab}_c`` that maps between ``c → a ⊗ b`` and ``c → b ⊗ a`` as in +``` +a -<-μ-<- c b -<-ν-<- c + ∨ -> Rsymbol(a,b,c)[μ,ν] v + b a +``` +If `FusionStyle(I)` is `UniqueFusion()` or `SimpleFusion()`, the R-symbol is a +number. Otherwise it is a square matrix with row and column size +`Nsymbol(a,b,c) == Nsymbol(b,a,c)`. +""" +function Rsymbol end + +# properties that can be determined in terms of the R symbol -Note that `Bosonic` and `Fermionic` are subtypes of `SymmetricBraiding`, which means that -braids are in fact equivalent to crossings (i.e. braiding twice is an identity: -`isone(Rsymbol(b,a,c)*Rsymbol(a,b,c)) == true`) and permutations are uniquely defined. """ -BraidingStyle(a::Sector) = BraidingStyle(typeof(a)) -BraidingStyle(::Type{Trivial}) = Bosonic() + twist(a::Sector) + +Return the twist of a sector `a` +""" +twist(a::Sector) = sum(dim(b) / dim(a) * tr(Rsymbol(a, a, b)) for b in a ⊗ a) # Pentagon and Hexagon equations #------------------------------------------------------------------------------- @@ -457,10 +427,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/TensorKitChainRulesCoreExt.jl b/ext/TensorKitChainRulesCoreExt/TensorKitChainRulesCoreExt.jl index da202178..1044b4ba 100644 --- a/ext/TensorKitChainRulesCoreExt/TensorKitChainRulesCoreExt.jl +++ b/ext/TensorKitChainRulesCoreExt/TensorKitChainRulesCoreExt.jl @@ -8,7 +8,7 @@ using LinearAlgebra using TupleTools import TensorOperations as TO -using TensorOperations: Backend, promote_contract +using TensorOperations: promote_contract using VectorInterface: promote_scale, promote_add ext = @static if isdefined(Base, :get_extension) @@ -16,7 +16,6 @@ ext = @static if isdefined(Base, :get_extension) else TensorOperations.TensorOperationsChainRulesCoreExt end -const _conj = ext._conj const trivtuple = ext.trivtuple include("utility.jl") diff --git a/ext/TensorKitChainRulesCoreExt/linalg.jl b/ext/TensorKitChainRulesCoreExt/linalg.jl index a7bd64fa..a814b3ad 100644 --- a/ext/TensorKitChainRulesCoreExt/linalg.jl +++ b/ext/TensorKitChainRulesCoreExt/linalg.jl @@ -42,7 +42,7 @@ function ChainRulesCore.rrule(::typeof(⊗), A::AbstractTensorMap, B::AbstractTe pB = (allind(B), ()) dA = zerovector(A, promote_contract(scalartype(ΔC), scalartype(B))) tB = twist(B, filter(x -> isdual(space(B, x)), allind(B))) - dA = tensorcontract!(dA, ipA, ΔC, pΔC, :N, tB, pB, :C) + dA = tensorcontract!(dA, ΔC, pΔC, false, tB, pB, true, ipA) return projectA(dA) end dB_ = @thunk begin @@ -50,7 +50,7 @@ function ChainRulesCore.rrule(::typeof(⊗), A::AbstractTensorMap, B::AbstractTe pA = ((), allind(A)) dB = zerovector(B, promote_contract(scalartype(ΔC), scalartype(A))) tA = twist(A, filter(x -> isdual(space(A, x)), allind(A))) - dB = tensorcontract!(dB, ipB, tA, pA, :C, ΔC, pΔC, :N) + dB = tensorcontract!(dB, tA, pA, true, ΔC, pΔC, false, ipB) return projectB(dB) end return NoTangent(), dA_, dB_ diff --git a/ext/TensorKitChainRulesCoreExt/tensoroperations.jl b/ext/TensorKitChainRulesCoreExt/tensoroperations.jl index 72c25a5a..43911e43 100644 --- a/ext/TensorKitChainRulesCoreExt/tensoroperations.jl +++ b/ext/TensorKitChainRulesCoreExt/tensoroperations.jl @@ -1,10 +1,45 @@ -function ChainRulesCore.rrule(::typeof(TO.tensorcontract!), - C::AbstractTensorMap{S}, pC::Index2Tuple, - A::AbstractTensorMap{S}, pA::Index2Tuple, conjA::Symbol, - B::AbstractTensorMap{S}, pB::Index2Tuple, conjB::Symbol, - α::Number, β::Number, - backend::Backend...) where {S} - C′ = tensorcontract!(copy(C), pC, A, pA, conjA, B, pB, conjB, α, β, backend...) +function ext._rrule_tensoradd!(C::AbstractTensorMap, A::AbstractTensorMap, pA::Index2Tuple, + conjA::Bool, α::Number, β::Number, ba) + C′ = tensoradd!(copy(C), A, pA, conjA, α, β, ba...) + + projectA = ProjectTo(A) + projectC = ProjectTo(C) + projectα = ProjectTo(α) + projectβ = ProjectTo(β) + + function pullback(ΔC′) + ΔC = unthunk(ΔC′) + dC = @thunk projectC(scale(ΔC, conj(β))) + dA = @thunk begin + ipA = invperm(linearize(pA)) + _dA = zerovector(A, promote_add(ΔC, α)) + _dA = tensoradd!(_dA, ΔC, (ipA, ()), conjA, conjA ? α : conj(α), Zero(), ba...) + return projectA(_dA) + end + dα = @thunk begin + # TODO: this is an inner product implemented as a contraction + # for non-symmetric tensors this might be more efficient like this, + # but for symmetric tensors an intermediate object will anyways be created + # and then it might be more efficient to use an addition and inner product + tΔC = twist(ΔC, filter(x -> isdual(space(ΔC, x)), allind(ΔC))) + _dα = tensorscalar(tensorcontract(A, ((), linearize(pA)), + !conjA, tΔC, + (trivtuple(TO.numind(pA)), ()), false, + ((), ()), One(), ba...)) + return projectα(_dα) + end + dβ = @thunk projectβ(inner(C, ΔC)) + return NoTangent(), dC, dA, NoTangent(), NoTangent(), dα, dβ + end + + return C′, pullback +end + +function ext._rrule_tensorcontract!(C::AbstractTensorMap, A::AbstractTensorMap, + pA::Index2Tuple, + conjA::Bool, B::AbstractTensorMap, pB::Index2Tuple, + conjB::Bool, pAB::Index2Tuple, α::Number, β::Number, ba) + C′ = tensorcontract!(copy(C), A, pA, conjA, B, pB, conjB, pAB, α, β, ba...) projectA = ProjectTo(A) projectB = ProjectTo(B) @@ -14,101 +49,58 @@ function ChainRulesCore.rrule(::typeof(TO.tensorcontract!), function pullback(ΔC′) ΔC = unthunk(ΔC′) - ipC = invperm(linearize(pC)) - pΔC = (TupleTools.getindices(ipC, trivtuple(TO.numout(pA))), - TupleTools.getindices(ipC, TO.numout(pA) .+ trivtuple(TO.numin(pB)))) + ipAB = invperm(linearize(pAB)) + pΔC = (TupleTools.getindices(ipAB, trivtuple(TO.numout(pA))), + TupleTools.getindices(ipAB, TO.numout(pA) .+ trivtuple(TO.numin(pB)))) dC = @thunk projectC(scale(ΔC, conj(β))) dA = @thunk begin ipA = (invperm(linearize(pA)), ()) - conjΔC = conjA == :C ? :C : :N - conjB′ = conjA == :C ? conjB : _conj(conjB) + conjΔC = conjA + conjB′ = conjA ? conjB : !conjB _dA = zerovector(A, promote_contract(scalartype(ΔC), scalartype(B), scalartype(α))) tB = twist(B, TupleTools.vcat(filter(x -> !isdual(space(B, x)), pB[1]), filter(x -> isdual(space(B, x)), pB[2]))) - _dA = tensorcontract!(_dA, ipA, + _dA = tensorcontract!(_dA, ΔC, pΔC, conjΔC, - tB, reverse(pB), conjB′, - conjA == :C ? α : conj(α), Zero(), backend...) + tB, reverse(pB), conjB′, ipA, + conjA ? α : conj(α), Zero(), ba...) return projectA(_dA) end dB = @thunk begin ipB = (invperm(linearize(pB)), ()) - conjΔC = conjB == :C ? :C : :N - conjA′ = conjB == :C ? conjA : _conj(conjA) + conjΔC = conjB + conjA′ = conjB ? conjA : !conjA _dB = zerovector(B, promote_contract(scalartype(ΔC), scalartype(A), scalartype(α))) tA = twist(A, TupleTools.vcat(filter(x -> isdual(space(A, x)), pA[1]), filter(x -> !isdual(space(A, x)), pA[2]))) - _dB = tensorcontract!(_dB, ipB, + _dB = tensorcontract!(_dB, tA, reverse(pA), conjA′, - ΔC, pΔC, conjΔC, - conjB == :C ? α : conj(α), Zero(), backend...) + ΔC, pΔC, conjΔC, ipB, + conjB ? α : conj(α), Zero(), ba...) return projectB(_dB) end dα = @thunk begin # TODO: this result should be AB = (C′ - βC) / α as C′ = βC + αAB - AB = tensorcontract(pC, A, pA, conjA, B, pB, conjB) + AB = tensorcontract(A, pA, conjA, B, pB, conjB, pAB, One(), ba...) return projectα(inner(AB, ΔC)) end dβ = @thunk projectβ(inner(C, ΔC)) - dbackend = map(x -> NoTangent(), backend) - return NoTangent(), dC, NoTangent(), - dA, NoTangent(), NoTangent(), dB, NoTangent(), NoTangent(), dα, dβ, - dbackend... + return NoTangent(), dC, + dA, NoTangent(), NoTangent(), + dB, NoTangent(), NoTangent(), + NoTangent(), dα, dβ end return C′, pullback end -function ChainRulesCore.rrule(::typeof(TO.tensoradd!), - C::AbstractTensorMap{S}, pC::Index2Tuple, - A::AbstractTensorMap{S}, conjA::Symbol, - α::Number, β::Number, backend::Backend...) where {S} - C′ = tensoradd!(copy(C), pC, A, conjA, α, β, backend...) - - projectA = ProjectTo(A) - projectC = ProjectTo(C) - projectα = ProjectTo(α) - projectβ = ProjectTo(β) - - function pullback(ΔC′) - ΔC = unthunk(ΔC′) - dC = @thunk projectC(scale(ΔC, conj(β))) - dA = @thunk begin - ipC = invperm(linearize(pC)) - _dA = zerovector(A, promote_add(ΔC, α)) - _dA = tensoradd!(_dA, (ipC, ()), ΔC, conjA, conjA == :N ? conj(α) : α, Zero(), - backend...) - return projectA(_dA) - end - dα = @thunk begin - # TODO: this is an inner product implemented as a contraction - # for non-symmetric tensors this might be more efficient like this, - # but for symmetric tensors an intermediate object will anyways be created - # and then it might be more efficient to use an addition and inner product - tΔC = twist(ΔC, filter(x -> isdual(space(ΔC, x)), allind(ΔC))) - _dα = tensorscalar(tensorcontract(((), ()), A, ((), linearize(pC)), - _conj(conjA), tΔC, - (trivtuple(TO.numind(pC)), - ()), :N, One(), backend...)) - return projectα(_dα) - end - dβ = @thunk projectβ(inner(C, ΔC)) - dbackend = map(x -> NoTangent(), backend) - return NoTangent(), dC, NoTangent(), dA, NoTangent(), dα, dβ, dbackend... - end - - return C′, pullback -end - -function ChainRulesCore.rrule(::typeof(tensortrace!), C::AbstractTensorMap{S}, - pC::Index2Tuple, A::AbstractTensorMap{S}, - pA::Index2Tuple, conjA::Symbol, α::Number, β::Number, - backend::Backend...) where {S} - C′ = tensortrace!(copy(C), pC, A, pA, conjA, α, β, backend...) +function ext._rrule_tensortrace!(C::AbstractTensorMap, A::AbstractTensorMap, p::Index2Tuple, + q::Index2Tuple, conjA::Bool, α::Number, β::Number, ba) + C′ = tensortrace!(copy(C), A, p, q, conjA, α, β, ba...) projectA = ProjectTo(A) projectC = ProjectTo(C) @@ -119,26 +111,24 @@ function ChainRulesCore.rrule(::typeof(tensortrace!), C::AbstractTensorMap{S}, ΔC = unthunk(ΔC′) dC = @thunk projectC(scale(ΔC, conj(β))) dA = @thunk begin - ipC = invperm((linearize(pC)..., pA[1]..., pA[2]...)) - E = one!(TO.tensoralloc_add(scalartype(A), pA, A, conjA)) + ip = invperm((linearize(p)..., q[1]..., q[2]...)) + E = one!(TO.tensoralloc_add(scalartype(A), A, q, conjA)) twist!(E, filter(x -> !isdual(space(E, x)), codomainind(E))) _dA = zerovector(A, promote_scale(ΔC, α)) - _dA = tensorproduct!(_dA, (ipC, ()), ΔC, - (trivtuple(TO.numind(pC)), ()), conjA, E, - ((), trivtuple(TO.numind(pA))), conjA, - conjA == :N ? conj(α) : α, Zero(), backend...) + _dA = tensorproduct!(_dA, ΔC, + (trivtuple(TO.numind(p)), ()), conjA, E, + ((), trivtuple(TO.numind(q))), conjA, (ip, ()), + conjA ? α : conj(α), Zero(), ba...) return projectA(_dA) end dα = @thunk begin # TODO: this result might be easier to compute as: # C′ = βC + α * trace(A) ⟹ At = (C′ - βC) / α - At = tensortrace(pC, A, pA, conjA) + At = tensortrace(A, p, q, conjA) return projectα(inner(At, ΔC)) end dβ = @thunk projectβ(inner(C, ΔC)) - dbackend = map(x -> NoTangent(), backend) - return NoTangent(), dC, NoTangent(), dA, NoTangent(), NoTangent(), dα, dβ, - dbackend... + return NoTangent(), dC, dA, NoTangent(), NoTangent(), NoTangent(), dα, dβ end return C′, pullback diff --git a/ext/TensorKitChainRulesCoreExt/utility.jl b/ext/TensorKitChainRulesCoreExt/utility.jl index 170ed09f..789286f5 100644 --- a/ext/TensorKitChainRulesCoreExt/utility.jl +++ b/ext/TensorKitChainRulesCoreExt/utility.jl @@ -17,8 +17,8 @@ end TensorKit.block(t::ZeroTangent, c::Sector) = t 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/ext/TensorKitFiniteDifferencesExt.jl b/ext/TensorKitFiniteDifferencesExt.jl index d4104a1a..bccaca75 100644 --- a/ext/TensorKitFiniteDifferencesExt.jl +++ b/ext/TensorKitFiniteDifferencesExt.jl @@ -1,7 +1,7 @@ module TensorKitFiniteDifferencesExt using TensorKit -using TensorKit: sqrtdim, isqrtdim +using TensorKit: sqrtdim, invsqrtdim using VectorInterface: scale! using FiniteDifferences @@ -18,7 +18,7 @@ function FiniteDifferences.to_vec(t::AbstractTensorMap) t′ = similar(t) xvec_of_vecs = back(x) for (i, (c, b)) in enumerate(blocks(t′)) - scale!(b, xvec_of_vecs[i], isqrtdim(c)) + scale!(b, xvec_of_vecs[i], invsqrtdim(c)) end return t′ end diff --git a/src/TensorKit.jl b/src/TensorKit.jl index a796c9e5..e1f62415 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -95,7 +95,7 @@ using Strided using VectorInterface using TensorOperations: TensorOperations, @tensor, @tensoropt, @ncon, ncon -using TensorOperations: IndexTuple, Index2Tuple, linearize, Backend +using TensorOperations: IndexTuple, Index2Tuple, linearize, AbstractBackend const TO = TensorOperations using LRUCache @@ -118,6 +118,8 @@ using LinearAlgebra: norm, dot, normalize, normalize!, tr, Diagonal, Hermitian import Base.Meta +using Random: Random + using PackageExtensionCompat # Auxiliary files @@ -163,11 +165,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/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 3cf0da46..4dbec263 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -50,6 +50,8 @@ else using Base: @constprop end +const MatOrNumber{T<:Number} = Union{DenseMatrix{T},T} + """ _interleave(a::NTuple{N}, b::NTuple{N}) -> NTuple{2N} diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl index f6ebe84f..c470eb6f 100644 --- a/src/fusiontrees/fusiontrees.jl +++ b/src/fusiontrees/fusiontrees.jl @@ -169,23 +169,16 @@ function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I,2}) where {I} a, b = f.uncoupled isduala, isdualb = f.isdual c = f.coupled - da, db, dc = dim.((a, b, c)) μ = (FusionStyle(I) isa GenericFusion) ? f.vertices[1] : 1 C = convert(A, fusiontensor(a, b, c))[:, :, :, μ] X = C if isduala - Xtemp = X - X = similar(Xtemp) Za = convert(A, FusionTree((a,), a, (isduala,), ())) - TO.tensorcontract!(X, ((1, 2, 3), ()), Za, ((1,), (2,)), :N, Xtemp, ((1,), (2, 3)), - :N, true, false) + @tensor X[a′, b, c] := Za[a′, a] * X[a, b, c] end if isdualb - Xtemp = X - X = similar(Xtemp) Zb = convert(A, FusionTree((b,), b, (isdualb,), ())) - TO.tensorcontract!(X, ((2, 1, 3), ()), Zb, ((1,), (2,)), :N, Xtemp, ((2,), (1, 3)), - :N, true, false) + @tensor X[a, b′, c] := Zb[b′, b] * X[a, b, c] end return X end @@ -203,8 +196,10 @@ function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I,N}) where {I,N} d1 = size(C1) X = similar(C1, (d1[1], d1[2], Base.tail(dtail)...)) trivialtuple = ntuple(identity, Val(N)) - return TO.tensorcontract!(X, ((trivialtuple..., N + 1), ()), C1, ((1, 2), (3,)), :N, - Ctail, ((1,), Base.tail(trivialtuple)), :N, true, false) + return TO.tensorcontract!(X, + C1, ((1, 2), (3,)), false, + Ctail, ((1,), Base.tail(trivialtuple)), false, + ((trivialtuple..., N + 1), ())) end # TODO: is this piracy? diff --git a/src/fusiontrees/manipulations.jl b/src/fusiontrees/manipulations.jl index dc857025..47051909 100644 --- a/src/fusiontrees/manipulations.jl +++ b/src/fusiontrees/manipulations.jl @@ -253,7 +253,7 @@ function bendright(f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {I< isdual2 = (f₂.isdual..., !(f₁.isdual[N₁])) inner2 = N₂ > 1 ? (f₂.innerlines..., c) : () - coeff₀ = sqrtdim(c) * isqrtdim(a) + coeff₀ = sqrtdim(c) * invsqrtdim(a) if f₁.isdual[N₁] coeff₀ *= conj(frobeniusschur(dual(b))) end diff --git a/src/planar/planaroperations.jl b/src/planar/planaroperations.jl index ee23bcd7..ea3161b7 100644 --- a/src/planar/planaroperations.jl +++ b/src/planar/planaroperations.jl @@ -1,34 +1,96 @@ # 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₂} - return add_transpose!(C, A, p, α, β, backend...) + +# insert default backend +function planaradd!(C, A, p::Index2Tuple, α::Number, β::Number) + return planaradd!(C, A, p, α, β, TO.DefaultBackend()) +end +# insert default allocator +function planaradd!(C, A, p::Index2Tuple, α::Number, β::Number, backend::AbstractBackend) + return planaradd!(C, A, p, α, β, backend, TO.DefaultAllocator()) +end +# replace default backend with select_backend mechanism +function planaradd!(C, A, p::Index2Tuple, α::Number, β::Number, backend::AbstractBackend, + allocator) + if backend isa TO.DefaultBackend + backend = TO.select_backend(planaradd!, C, A) + return planaradd!(C, A, p, α, β, backend, allocator) + elseif backend isa TO.NoBackend + # error for missing backend + TC = typeof(C) + TA = typeof(A) + throw(ArgumentError("No suitable backend found for planaradd! and tensor types $TC and $TA")) + else + # error for unknown backend + TC = typeof(C) + TA = typeof(A) + throw(ArgumentError("Unknown backend $backend for planaradd! and tensor types $TC and $TA")) + end +end +# implementation +function planaradd!(C::AbstractTensorMap, + A::AbstractTensorMap, + p::Index2Tuple, + α::Number, β::Number, + backend::AbstractBackend, allocator) + 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₃}, +# insert default backend +function planartrace!(C, A, p::Index2Tuple, q::Index2Tuple, α::Number, β::Number) + return planartrace!(C, A, p, q, α, β, TO.DefaultBackend()) +end +# insert default allocator +function planartrace!(C, A, p::Index2Tuple, q::Index2Tuple, α::Number, β::Number, + backend::AbstractBackend) + return planartrace!(C, A, p, q, α, β, backend, TO.DefaultAllocator()) +end +# replace default backend with select_backend mechanism +function planartrace!(C, A, p::Index2Tuple, q::Index2Tuple, α::Number, β::Number, + backend::AbstractBackend, allocator) + if backend isa TO.DefaultBackend + backend = TO.select_backend(planartrace!, C, A) + return planartrace!(C, A, p, q, α, β, backend, allocator) + elseif backend isa TO.NoBackend + # error for missing backend + TC = typeof(C) + TA = typeof(A) + throw(ArgumentError("No suitable backend found for planartrace! and tensor types $TC and $TA")) + else + # error for unknown backend + TC = typeof(C) + TA = typeof(A) + throw(ArgumentError("Unknown backend $backend for planartrace! and tensor types $TC and $TA")) + end +end +# implementation +function planartrace!(C::AbstractTensorMap, + A::AbstractTensorMap, + (p₁, p₂)::Index2Tuple, + (q₁, q₂)::Index2Tuple, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂,N₃} + backend::AbstractBackend, allocator) + (S = spacetype(C)) == spacetype(A) || + throw(SpaceMismatch("incompatible spacetypes")) if BraidingStyle(sectortype(S)) == Bosonic() - return trace_permute!(C, A, p, q, α, β, backend...) + 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(β) @@ -36,25 +98,60 @@ function planartrace!(C::AbstractTensorMap{S,N₁,N₂}, elseif !isone(β) rmul!(C, β) end + β′ = One() 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₂′], A[f₁, f₂], (p₁, p₂), (q₁, q₂), false, α * coeff, + β′, + backend, allocator) end end return C end -function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, - A::AbstractTensorMap{S}, +# insert default backend +function planarcontract!(C, A, pA::Index2Tuple, B, pB::Index2Tuple, pAB::Index2Tuple, + α::Number, β::Number) + return planarcontract!(C, A, pA, B, pB, pAB, α, β, TO.DefaultBackend()) +end +# insert default allocator +function planarcontract!(C, A, pA::Index2Tuple, B, pB::Index2Tuple, pAB::Index2Tuple, + α::Number, β::Number, backend::AbstractBackend) + return planarcontract!(C, A, pA, B, pB, pAB, α, β, backend, TO.DefaultAllocator()) +end +# replace default backend with select_backend mechanism +function planarcontract!(C, A, pA::Index2Tuple, B, pB::Index2Tuple, pAB::Index2Tuple, + α::Number, β::Number, backend::AbstractBackend, allocator) + if backend isa TO.DefaultBackend + backend = TO.select_backend(planarcontract!, C, A, B) + return planarcontract!(C, A, pA, B, pB, pAB, α, β, backend, allocator) + elseif backend isa TO.NoBackend + # error for missing backend + TC = typeof(C) + TA = typeof(A) + TB = typeof(B) + throw(ArgumentError("No suitable backend found for planarcontract! and tensor types $TC, $TA and $TB")) + else + # error for unknown backend + TC = typeof(C) + TA = typeof(A) + TB = typeof(B) + throw(ArgumentError("Unknown backend $backend for planarcontract! and tensor types $TC, $TA and $TB")) + end +end +# implementation +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() - return contract!(C, A, pA, B, pB, pAB, α, β, backend...) + backend::AbstractBackend, + allocator) + if BraidingStyle(sectortype(C)) == Bosonic() + return contract!(C, A, pA, B, pB, pAB, α, β, backend, allocator) end codA, domA = codomainind(A), domainind(A) @@ -67,19 +164,21 @@ function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, if oindA == codA && cindA == domA A′ = A else - A′ = TO.tensoralloc_add(scalartype(A), (oindA, cindA), A, :N, true) - add_transpose!(A′, A, (oindA, cindA), true, false, backend...) + A′ = TO.tensoralloc_add(scalartype(A), A, (oindA, cindA), false, Val(true), + allocator) + add_transpose!(A′, A, (oindA, cindA), One(), Zero(), backend) end if cindB == codB && oindB == domB B′ = B else - B′ = TensorOperations.tensoralloc_add(scalartype(B), (cindB, oindB), B, :N, true) - add_transpose!(B′, B, (cindB, oindB), true, false, backend...) + B′ = TensorOperations.tensoralloc_add(scalartype(B), B, (cindB, oindB), false, + Val(true), allocator) + add_transpose!(B′, B, (cindB, oindB), One(), Zero(), backend) end mul!(C, A′, B′, α, β) - (oindA == codA && cindA == domA) || TO.tensorfree!(A′) - (cindB == codB && oindB == domB) || TO.tensorfree!(B′) + (oindA == codA && cindA == domA) || TO.tensorfree!(A′, allocator) + (cindB == codB && oindB == domB) || TO.tensorfree!(B′, allocator) return C end diff --git a/src/planar/postprocessors.jl b/src/planar/postprocessors.jl index 69100f1a..37eff61d 100644 --- a/src/planar/postprocessors.jl +++ b/src/planar/postprocessors.jl @@ -12,7 +12,7 @@ function _annotate_temporaries(ex, temporaries) if i !== nothing rhs = ex.args[2] # add `istemp = true` flag - newrhs = Expr(:call, rhs.args[1:(end - 1)]..., true) + newrhs = Expr(:call, rhs.args[1:(end - 1)]..., Val(true)) return Expr(:(=), lhs, newrhs) end elseif ex isa Expr @@ -51,32 +51,26 @@ end # Replace the tensor operations created by `TO.instantiate` with the corresponding # planar operations, immediately inserting them with `GlobalRef`. -# NOTE: work around a somewhat unfortunate interface choice in TensorOperations, which we will correct in the future. -_planaradd!(C, p, A, α, β, backend...) = planaradd!(C, A, p, α, β, backend...) -_planartrace!(C, p, A, q, α, β, backend...) = planartrace!(C, A, p, q, α, β, backend...) -function _planarcontract!(C, pAB, A, pA, B, pB, α, β, backend...) - return planarcontract!(C, A, pA, B, pB, pAB, α, β, backend...) -end # TODO: replace _planarmethod with planarmethod in everything below -const _PLANAR_OPERATIONS = (:_planaradd!, :_planartrace!, :_planarcontract!) +const _PLANAR_OPERATIONS = (:planaradd!, :planartrace!, :planarcontract!) function _insert_planar_operations(ex) if isexpr(ex, :call) if ex.args[1] == GlobalRef(TensorOperations, :tensoradd!) conjA = popat!(ex.args, 5) - @assert conjA == :(:N) "conj flag should be `:N` ($conjA)" - return Expr(ex.head, GlobalRef(TensorKit, Symbol(:_planaradd!)), + @assert !conjA "conj flag should be disabled" + return Expr(ex.head, GlobalRef(TensorKit, Symbol(:planaradd!)), map(_insert_planar_operations, ex.args[2:end])...) elseif ex.args[1] == GlobalRef(TensorOperations, :tensorcontract!) - conjB = popat!(ex.args, 9) - conjA = popat!(ex.args, 6) - @assert conjA == conjB == :(:N) "conj flag should be `:N` ($conjA), ($conjB)" - return Expr(ex.head, GlobalRef(TensorKit, Symbol(:_planarcontract!)), + conjB = popat!(ex.args, 8) + conjA = popat!(ex.args, 5) + @assert !conjA && !conjB "conj flags should be disabled ($conjA), ($conjB)" + return Expr(ex.head, GlobalRef(TensorKit, Symbol(:planarcontract!)), map(_insert_planar_operations, ex.args[2:end])...) elseif ex.args[1] == GlobalRef(TensorOperations, :tensortrace!) conjA = popat!(ex.args, 6) - @assert conjA == :(:N) "conj flag should be `:N` ($conjA)" - return Expr(ex.head, GlobalRef(TensorKit, Symbol(:_planartrace!)), + @assert !conjA "conj flag should be disabled" + return Expr(ex.head, GlobalRef(TensorKit, Symbol(:planartrace!)), map(_insert_planar_operations, ex.args[2:end])...) elseif ex.args[1] in TensorOperations.tensoroperationsfunctions return Expr(ex.head, GlobalRef(TensorOperations, ex.args[1]), @@ -88,16 +82,24 @@ function _insert_planar_operations(ex) return ex end -# Mimick `TO.insert_operationbackend` for planar operations. -function insert_operationbackend(ex, backend) - if isexpr(ex, :call) && ex.args[1] isa GlobalRef && - ex.args[1].mod == TensorKit && - ex.args[1].name ∈ _PLANAR_OPERATIONS - b = Backend{backend}() - return Expr(:call, ex.args..., b) - elseif isa(ex, Expr) - return Expr(ex.head, (insert_operationbackend(e, backend) for e in ex.args)...) - else - return ex - end +""" + insertplanarbackend(ex, backend) + +Insert the backend argument into the tensor operation methods `planaradd!`, `planartrace!`, and `planarcontract!`. + +See also: [`TensorOperations.insertbackend`](@ref). +""" +function insertplanarbackend(ex, backend) + return TO.insertargument(ex, backend, (:planaradd!, :planartrace!, :planarcontract!)) +end + +""" + insertplanarallocator(ex, allocator) + +Insert the allocator argument into the tensor operation methods `planaradd!`, `planartrace!`, and `planarcontract!`. + +See also: [`TensorOperations.insertallocator`](@ref). +""" +function insertplanarallocator(ex, allocator) + return TO.insertargument(ex, allocator, (:planaradd!, :planartrace!, :planarcontract!)) 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 38997c94..1ee401c2 100644 --- a/src/tensors/abstracttensor.jl +++ b/src/tensors/abstracttensor.jl @@ -3,71 +3,71 @@ # Abstract Tensor type #---------------------- """ - abstract type AbstractTensorMap{S<:IndexSpace, N₁, N₂} end + abstract type AbstractTensorMap{T<: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 `T`. 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{T<:Number,S<:IndexSpace,N₁,N₂} end + """ - AbstractTensor{S<:IndexSpace, N} = AbstractTensorMap{S, N, 0} + AbstractTensor{T,S,N} = AbstractTensorMap{T,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 `T`. -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{T, S, N}` is actually a special case `AbstractTensorMap{T, 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{T,S,N} = AbstractTensorMap{T,S,N,0} # tensor characteristics #------------------------ -Base.eltype(::Union{T,Type{T}}) where {T<:AbstractTensorMap} = scalartype(T) +Base.eltype(::Type{<:AbstractTensorMap{T}}) where {T} = T """ - spacetype(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Type{S<:IndexSpace} + spacetype(::Union{TT,Type{TT}}) where {TT<: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{<:Any,S}}) where {S} = S """ - sectortype(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Type{I<:Sector} + sectortype(::Union{TT,Type{TT}}) where {TT<:AbstractTensorMap} -> Type{I<:Sector} Return the type of sector `I` of a tensor. """ -sectortype(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = sectortype(S) +sectortype(::Type{TT}) where {TT<:AbstractTensorMap} = sectortype(spacetype(TT)) -function InnerProductStyle(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} - return InnerProductStyle(S) +function InnerProductStyle(::Type{TT}) where {TT<:AbstractTensorMap} + return InnerProductStyle(spacetype(TT)) end """ - field(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Type{𝕂<:Field} + field(::Union{TT,Type{TT}}) where {TT<:AbstractTensorMap} -> Type{𝕂<:Field} Return the type of field `𝕂` of a tensor. """ -field(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = field(S) +field(::Type{TT}) where {TT<:AbstractTensorMap} = field(spacetype(TT)) """ - numout(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Int + numout(::Union{TT,Type{TT}}) where {TT<:AbstractTensorMap} -> Int Return the number of output spaces of a tensor. This is equivalent to the number of spaces in the codomain of that tensor. See also [`numin`](@ref) and [`numind`](@ref). """ -numout(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} = N₁ +numout(::Type{<:AbstractTensorMap{T,S,N₁}}) where {T,S,N₁} = N₁ """ - numin(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Int + numin(::Union{TT,Type{TT}}) where {TT<:AbstractTensorMap} -> Int Return the number of input spaces of a tensor. This is equivalent to the number of spaces in the domain of that tensor. See also [`numout`](@ref) and [`numind`](@ref). """ -numin(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} = N₂ +numin(::Type{<:AbstractTensorMap{T,S,N₁,N₂}}) where {T,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{TT}) where {TT<:AbstractTensorMap} = numin(TT) + numout(TT) function similarstoragetype(TT::Type{<:AbstractTensorMap}, ::Type{T}) where {T} return Core.Compiler.return_type(similar, Tuple{storagetype(TT),Type{T}}) @@ -92,12 +92,13 @@ numin(t::AbstractTensorMap) = numin(typeof(t)) numind(t::AbstractTensorMap) = numind(typeof(t)) storagetype(t::AbstractTensorMap) = storagetype(typeof(t)) -similarstoragetype(t::AbstractTensorMap, T) = similarstoragetype(typeof(t), T) +similarstoragetype(t::AbstractTensorMap, TT) = similarstoragetype(typeof(t), TT) const order = numind @doc """ - codomain(t::AbstractTensorMap{S,N₁,N₂}, [i::Int]) -> ProductSpace{S,N₁} + codomain(t::AbstractTensorMap{T,S,N₁,N₂}) -> ProductSpace{S,N₁} + codomain(t::AbstractTensorMap{T,S,N₁,N₂}, i::Int) -> S 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 +110,8 @@ 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{T,S,N₁,N₂}) -> ProductSpace{S,N₂} + domain(t::AbstractTensorMap{T,S,N₁,N₂}, i::Int) -> S 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 +123,8 @@ 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{T,S,N₁,N₂}) -> HomSpace{S,N₁,N₂} + space(t::AbstractTensorMap{T,S,N₁,N₂}, i::Int) -> S 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. """ @@ -137,43 +140,43 @@ symmetry. This is also the dimension of the `HomSpace` on which the `TensorMap` dim(t::AbstractTensorMap) = dim(space(t)) """ - codomainind(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Tuple{Int} + codomainind(::Union{TT,Type{TT}}) where {TT<:AbstractTensorMap} -> Tuple{Int} 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{TT}) where {TT<:AbstractTensorMap} + return ntuple(identity, numout(TT)) end codomainind(t::AbstractTensorMap) = codomainind(typeof(t)) """ - domainind(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Tuple{Int} + domainind(::Union{TT,Type{TT}}) where {TT<:AbstractTensorMap} -> Tuple{Int} 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{TT}) where {TT<:AbstractTensorMap} + return ntuple(n -> numout(TT) + n, numin(TT)) end domainind(t::AbstractTensorMap) = domainind(typeof(t)) """ - allind(::Union{T,Type{T}}) where {T<:AbstractTensorMap} -> Tuple{Int} + allind(::Union{TT,Type{TT}}) where {TT<:AbstractTensorMap} -> Tuple{Int} 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{TT}) where {TT<:AbstractTensorMap} + return ntuple(identity, numind(TT)) 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) @@ -226,6 +229,55 @@ Return the dimensions of the block of a tensor corresponding to a coupled sector Return an iterator over all splitting - fusion tree pairs of a tensor. """ fusiontrees(::AbstractTensorMap) +# Similar +#--------- +# The implementation is written for similar(t, TorA, V::TensorMapSpace) -> TensorMap +# and all other methods are just filling in default arguments +# 4 arguments +@doc """ + similar(t::AbstractTensorMap, [AorT=storagetype(t)], [V=space(t)]) + similar(t::AbstractTensorMap, [AorT=storagetype(t)], codomain, domain) + +Creates an uninitialized mutable tensor with the given scalar or storagetype `AorT` and +structure `V` or `codomain ← domain`, based on the source tensormap. The second and third +arguments are both optional, defaulting to the given tensor's `storagetype` and `space`. +The structure may be specified either as a single `HomSpace` argument or as `codomain` and +`domain`. + +By default, this will result in `TensorMap{T}(undef, V)` when custom objects do not +specialize this method. +""" Base.similar(::AbstractTensorMap, args...) + +function Base.similar(t::AbstractTensorMap, ::Type{T}, codomain::TensorSpace{S}, + domain::TensorSpace{S}) where {T,S} + return similar(t, T, codomain ← domain) +end +# 3 arguments +function Base.similar(t::AbstractTensorMap, codomain::TensorSpace{S}, + domain::TensorSpace{S}) where {S} + return similar(t, storagetype(t), codomain ← domain) +end +function Base.similar(t::AbstractTensorMap, ::Type{T}, codomain::TensorSpace) where {T} + return similar(t, T, codomain ← one(codomain)) +end +# 2 arguments +function Base.similar(t::AbstractTensorMap, codomain::TensorSpace) + return similar(t, storagetype(t), codomain ← one(codomain)) +end +Base.similar(t::AbstractTensorMap, P::TensorMapSpace) = similar(t, storagetype(t), P) +Base.similar(t::AbstractTensorMap, ::Type{T}) where {T} = similar(t, T, space(t)) +# 1 argument +Base.similar(t::AbstractTensorMap) = similar(t, storagetype(t), space(t)) + +# generic implementation for AbstractTensorMap -> returns `TensorMap` +function Base.similar(::AbstractTensorMap, ::Type{TorA}, + P::TensorMapSpace{S}) where {TorA<:MatOrNumber,S} + N₁ = length(codomain(P)) + N₂ = length(domain(P)) + TT = tensormaptype(S, N₁, N₂, TorA) + return TT(undef, codomain(P), domain(P)) +end + # Equality and approximality #---------------------------- function Base.:(==)(t1::AbstractTensorMap, t2::AbstractTensorMap) @@ -258,7 +310,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 14776038..62f36e75 100644 --- a/src/tensors/adjoint.jl +++ b/src/tensors/adjoint.jl @@ -1,26 +1,29 @@ # AdjointTensorMap: lazy adjoint #==========================================================# """ - struct AdjointTensorMap{S<:IndexSpace, N₁, N₂, ...} <: AbstractTensorMap{S, N₁, N₂} + struct AdjointTensorMap{T, S, N₁, N₂, ...} <: AbstractTensorMap{T, 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{T,S,N₁,N₂,I,A,F₁,F₂} <: + AbstractTensorMap{T,S,N₁,N₂} + parent::TensorMap{T,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{T,S,N₁,N₂,A<:DenseMatrix} = + AdjointTensorMap{T,S,N₁,N₂,Trivial,A,Nothing,Nothing} #! format: on # Constructor: construct from taking adjoint of a tensor Base.adjoint(t::TensorMap) = AdjointTensorMap(t) Base.adjoint(t::AdjointTensorMap) = t.parent -Base.similar(t::AdjointTensorMap, T::Type, P::TensorMapSpace) = similar(t', T, P) +function Base.similar(t::AdjointTensorMap, ::Type{TorA}, + P::TensorMapSpace) where {TorA<:MatOrNumber} + return similar(t', TorA, P) +end # Properties codomain(t::AdjointTensorMap) = domain(t.parent) @@ -29,8 +32,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{T,S,N₁,N₂,A}}) where {T,S,N₁,N₂,A<:DenseMatrix} = A +storagetype(::Type{<:AdjointTensorMap{T,S,N₁,N₂,I,<:SectorDict{I,A}}}) where {T,S,N₁,N₂,I<:Sector,A<:DenseMatrix} = A #! format: on dim(t::AdjointTensorMap) = dim(t.parent) @@ -44,8 +47,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{T,S,N₁,N₂,I}, + f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {T,S,N₁,N₂,I} c = f₁.coupled @boundscheck begin c == f₂.coupled || throw(SectorMismatch()) @@ -55,9 +58,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{T,S,N₁,N₂,I}, v, f₁::FusionTree{I,N₁}, - f₂::FusionTree{I,N₂}) where {S,N₁,N₂,I} + f₂::FusionTree{I,N₂}) where {T,S,N₁,N₂,I} return copy!(getindex(t, f₁, f₂), v) end @@ -90,16 +93,16 @@ end function Base.summary(io::IO, t::AdjointTensorMap) return print(io, "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..1ba27d49 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{T,S<:IndexSpace} <: AbstractTensorMap{T, 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{T,S} <: AbstractTensorMap{T,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{T,S}(V1::S, V2::S, adjoint::Bool=false) where {T,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{T,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{T,S}) where {T,S} + return BraidingTensor{T,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{T,S}}) where {T,S} = Matrix{T} 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) @@ -125,9 +125,9 @@ end end end -Base.similar(::BraidingTensor, T::Type, P::TensorMapSpace) = TensorMap(undef, T, P) +# efficient copy constructor +Base.copy(b::BraidingTensor) = b -Base.copy(b::BraidingTensor) = copy!(similar(b), b) function Base.copy!(t::TensorMap, b::BraidingTensor) space(t) == space(b) || throw(SectorMismatch()) fill!(t, zero(scalartype(t))) @@ -148,8 +148,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 +192,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::AbstractBackend...) + return add_transform!(tdst, TensorMap(tsrc), (p₁, p₂), fusiontreetransform, α, β, + backend...) end # VectorInterface @@ -209,30 +210,40 @@ 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, + A::BraidingTensor, pA::Index2Tuple, conjA::Symbol, + α::Number, β::Number, backend::AbstractBackend=TO.DefaultBackend(), + allocator=TO.DefaultAllocator()) + return TO.tensoradd!(C, TensorMap(A), pA, conjA, α, β, backend, allocator) 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::AbstractBackend, + allocator) + return planaradd!(C, TensorMap(A), p, α, β, backend, allocator) 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::AbstractBackend, + allocator) + # 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, allocator) + codA, domA = codomainind(A), domainind(A) codB, domB = codomainind(B), domainind(B) oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, @@ -244,7 +255,7 @@ function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, end if BraidingStyle(sectortype(B)) isa Bosonic - return add_permute!(C, B, (reverse(cindB), oindB), α, β, backend...) + return add_permute!(C, B, (reverse(cindB), oindB), α, β, backend) end τ_levels = A.adjoint ? (1, 2, 2, 1) : (2, 1, 1, 2) @@ -265,20 +276,26 @@ function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, end end for ((f₁′, f₂′), coeff) in newtrees - TO.tensoradd!(C[f₁′, f₂′], (reverse(cindB), oindB), B[f₁, f₂], :N, α * coeff, - One(), backend...) + TO.tensoradd!(C[f₁′, f₂′], B[f₁, f₂], (reverse(cindB), oindB), false, α * coeff, + One(), backend, allocator) end 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::AbstractBackend, + allocator) + # 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, allocator) + codA, domA = codomainind(A), domainind(A) codB, domB = codomainind(B), domainind(B) oindA, cindA, oindB, cindB = reorder_indices(codA, domA, codB, domB, oindA, cindA, @@ -290,7 +307,7 @@ function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, end if BraidingStyle(sectortype(A)) isa Bosonic - return add_permute!(C, A, (oindA, reverse(cindA)), α, β, backend...) + return add_permute!(C, A, (oindA, reverse(cindA)), α, β, backend) end scale!(C, β) @@ -311,48 +328,29 @@ function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, end end for ((f₁′, f₂′), coeff) in newtrees - TO.tensoradd!(C[f₁′, f₂′], (oindA, reverse(cindA)), A[f₁, f₂], :N, α * coeff, - One(), backend...) + TO.tensoradd!(C[f₁′, f₂′], A[f₁, f₂], (oindA, reverse(cindA)), false, α * coeff, + One(), backend, allocator) end 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, pAB::Index2Tuple, + α::Number, β::Number, backend::AbstractBackend, + allocator) + return planarcontract!(C, TensorMap(A), pA, TensorMap(B), pB, pAB, α, β, backend, + allocator) 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::AbstractBackend, + allocator) + return planartrace!(C, TensorMap(A), p, q, α, β, backend, allocator) end # function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl index afca5b0d..6d8e791b 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₂} - return add_permute!(tdst, tsrc, p, true, false) +@propagate_inbounds function Base.permute!(tdst::AbstractTensorMap, + tsrc::AbstractTensorMap, + p::Index2Tuple) + return add_permute!(tdst, tsrc, p, One(), Zero()) 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,8 +27,8 @@ 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₂} +function permute(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂}; + copy::Bool=false) where {N₁,N₂} space′ = permute(space(t), (p₁, p₂)) # share data if possible if !copy @@ -45,14 +44,13 @@ function permute(t::AbstractTensorMap{S}, (p₁, p₂)::Index2Tuple{N₁,N₂}; return permute!(similar(t, space′), 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, One(), Zero()) 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,10 +106,10 @@ 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) @@ -129,8 +127,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`. @@ -143,13 +141,13 @@ See [`transpose`](@ref) for creating a new tensor and [`add_transpose!`](@ref) f function LinearAlgebra.transpose!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple=_transpose_indices(t)) - return add_transpose!(tdst, tsrc, (p₁, p₂), true, false) + return add_transpose!(tdst, tsrc, (p₁, p₂), One(), Zero()) 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. @@ -160,10 +158,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 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) @@ -176,9 +174,9 @@ function LinearAlgebra.transpose(t::AbstractTensorMap{S}, 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)) @@ -266,23 +264,23 @@ twist(t::AbstractTensorMap, i; inv::Bool=false) = twist!(copy(t), i; inv) #------------------------------------- # Full implementations based on `add` #------------------------------------- -@propagate_inbounds function add_permute!(tdst::AbstractTensorMap{S,N₁,N₂}, +@propagate_inbounds function add_permute!(tdst::AbstractTensorMap{T,S,N₁,N₂}, tsrc::AbstractTensorMap, p::Index2Tuple{N₁,N₂}, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} + backend::AbstractBackend...) where {T,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{T,S,N₁,N₂}, tsrc::AbstractTensorMap, p::Index2Tuple{N₁,N₂}, levels::IndexTuple, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} + backend::AbstractBackend...) where {T,S,N₁,N₂} length(levels) == numind(tsrc) || throw(ArgumentError("incorrect levels $levels for tensor map $(codomain(tsrc)) ← $(domain(tsrc))")) @@ -293,23 +291,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{T,S,N₁,N₂}, tsrc::AbstractTensorMap, p::Index2Tuple{N₁,N₂}, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} + backend::AbstractBackend...) where {T,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{T,S,N₁,N₂}, tsrc::AbstractTensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂}, fusiontreetransform, α::Number, β::Number, - backend::Backend...) where {S,N₁,N₂} + backend::AbstractBackend...) where {T,S,N₁,N₂} @boundscheck begin permute(space(tsrc), (p₁, p₂)) == space(tdst) || throw(SpaceMismatch("source = $(codomain(tsrc))←$(domain(tsrc)), @@ -331,7 +329,7 @@ end # internal methods: no argument types function _add_trivial_kernel!(tdst, tsrc, p, fusiontreetransform, α, β, backend...) - TO.tensoradd!(tdst[], p, tsrc[], :N, α, β, backend...) + TO.tensoradd!(tdst[], tsrc[], p, false, α, β, backend...) return nothing end @@ -352,7 +350,8 @@ end function _add_abelian_block!(tdst, tsrc, p, fusiontreetransform, f₁, f₂, α, β, backend...) (f₁′, f₂′), coeff = first(fusiontreetransform(f₁, f₂)) - @inbounds TO.tensoradd!(tdst[f₁′, f₂′], p, tsrc[f₁, f₂], :N, α * coeff, β, backend...) + @inbounds TO.tensoradd!(tdst[f₁′, f₂′], tsrc[f₁, f₂], p, false, α * coeff, β, + backend...) return nothing end @@ -362,28 +361,30 @@ function _add_general_kernel!(tdst, tsrc, p, fusiontreetransform, α, β, backen elseif β != 1 tdst = scale!(tdst, β) end + β′ = One() if Threads.nthreads() > 1 Threads.@sync for s₁ in sectors(codomain(tsrc)), s₂ in sectors(domain(tsrc)) Threads.@spawn _add_nonabelian_sector!(tdst, tsrc, p, fusiontreetransform, s₁, - s₂, α, β, backend...) + s₂, α, β′, backend...) end else for (f₁, f₂) in fusiontrees(tsrc) for ((f₁′, f₂′), coeff) in fusiontreetransform(f₁, f₂) - @inbounds TO.tensoradd!(tdst[f₁′, f₂′], p, tsrc[f₁, f₂], :N, α * coeff, - true, backend...) + @inbounds TO.tensoradd!(tdst[f₁′, f₂′], tsrc[f₁, f₂], p, false, α * coeff, + β′, backend...) end end end return nothing end +# TODO: β argument is weird here because it has to be 1 function _add_nonabelian_sector!(tdst, tsrc, p, fusiontreetransform, s₁, s₂, α, β, backend...) for (f₁, f₂) in fusiontrees(tsrc) (f₁.uncoupled == s₁ && f₂.uncoupled == s₂) || continue for ((f₁′, f₂′), coeff) in fusiontreetransform(f₁, f₂) - @inbounds TO.tensoradd!(tdst[f₁′, f₂′], p, tsrc[f₁, f₂], :N, α * coeff, true, + @inbounds TO.tensoradd!(tdst[f₁′, f₂′], tsrc[f₁, f₂], p, false, α * coeff, β, backend...) end end diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index 235c9507..88199850 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -53,106 +53,104 @@ function one!(t::AbstractTensorMap) end """ - id([A::Type{<:DenseMatrix} = Matrix{Float64},] space::VectorSpace) -> TensorMap + id([T::Type=Float64,] V::TensorSpace) -> TensorMap -Construct the identity endomorphism on space `space`, i.e. return a `t::TensorMap` with `domain(t) == codomain(t) == V`, where `storagetype(t) = A` can be specified. +Construct the identity endomorphism on space `V`, i.e. return a `t::TensorMap` with +`domain(t) == codomain(t) == V`, where either `scalartype(t) = T` if `T` is a `Number` type +or `storagetype(t) = T` if `T` is a `DenseMatrix` type. """ -id(A, V::ElementarySpace) = id(A, ProductSpace(V)) -id(V::VectorSpace) = id(Matrix{Float64}, V) -function id(::Type{A}, P::ProductSpace) where {A<:DenseMatrix} - return one!(TensorMap(s -> A(undef, s), P, P)) +id(V::TensorSpace) = id(Float64, V) +function id(::Type{A}, V::TensorSpace{S}) where {A<:MatOrNumber,S} + W = V ← V + N = length(domain(W)) + t = tensormaptype(S, N, N, A)(undef, codomain(W), domain(W)) + return one!(t) end """ - isomorphism([A::Type{<:DenseMatrix} = Matrix{Float64},] - cod::VectorSpace, dom::VectorSpace) - -> TensorMap + isomorphism([T::Type=Float64,] codomain::TensorSpace, domain::TensorSpace) -> TensorMap + isomorphism([T::Type=Float64,] codomain ← domain) -> TensorMap + isomorphism([T::Type=Float64,] domain → codomain) -> TensorMap -Return a `t::TensorMap` that implements a specific isomorphism between the codomain `cod` -and the domain `dom`, and for which `storagetype(t)` can optionally be chosen to be of type -`A`. If the two spaces do not allow for such an isomorphism, and are thus not isomorphic, -and error will be thrown. When they are isomorphic, there is no canonical choice for a -specific isomorphism, but the current choice is such that -`isomorphism(cod, dom) == inv(isomorphism(dom, cod))`. +Construct a specific isomorphism between the codomain and the domain, i.e. return a +`t::TensorMap` where either `scalartype(t) = T` if `T` is a `Number` type or +`storagetype(t) = T` if `T` is a `DenseMatrix` type. If the spaces are not isomorphic, an +error will be thrown. + +!!! note + There is no canonical choice for a specific isomorphism, but the current choice is such + that `isomorphism(cod, dom) == inv(isomorphism(dom, cod))`. See also [`unitary`](@ref) when `InnerProductStyle(cod) === EuclideanProduct()`. """ -isomorphism(cod::TensorSpace, dom::TensorSpace) = isomorphism(Matrix{Float64}, cod, dom) -isomorphism(P::TensorMapSpace) = isomorphism(codomain(P), domain(P)) -function isomorphism(A::Type{<:DenseMatrix}, P::TensorMapSpace) - return isomorphism(A, codomain(P), domain(P)) -end -function isomorphism(A::Type{<:DenseMatrix}, cod::TensorSpace, dom::TensorSpace) - return isomorphism(A, convert(ProductSpace, cod), convert(ProductSpace, dom)) -end -function isomorphism(::Type{A}, cod::ProductSpace, dom::ProductSpace) where {A<:DenseMatrix} - cod ≅ dom || throw(SpaceMismatch("codomain $cod and domain $dom are not isomorphic")) - t = TensorMap(s -> A(undef, s), cod, dom) - for (c, b) in blocks(t) +function isomorphism(::Type{A}, V::TensorMapSpace{S,N₁,N₂}) where {A<:MatOrNumber,S,N₁,N₂} + codomain(V) ≅ domain(V) || + throw(SpaceMismatch("codomain and domain are not isomorphic: $V")) + t = tensormaptype(S, N₁, N₂, A)(undef, codomain(V), domain(V)) + for (_, b) in blocks(t) MatrixAlgebra.one!(b) end return t end """ - unitary([A::Type{<:DenseMatrix} = Matrix{Float64},] cod::VectorSpace, dom::VectorSpace) - -> TensorMap - -Return a `t::TensorMap` that implements a specific unitary isomorphism between the codomain -`cod` and the domain `dom`, for which `spacetype(dom)` (`== spacetype(cod)`) must have an -inner product. Furthermore, `storagetype(t)` can optionally be chosen to be -of type `A`. If the two spaces do not allow for such an isomorphism, and are thus not -isomorphic, and error will be thrown. When they are isomorphic, there is no canonical choice -for a specific isomorphism, but the current choice is such that -`unitary(cod, dom) == inv(unitary(dom, cod)) = adjoint(unitary(dom, cod))`. + unitary([T::Type=Float64,] codomain::TensorSpace, domain::TensorSpace) -> TensorMap + unitary([T::Type=Float64,] codomain ← domain) -> TensorMap + unitary([T::Type=Float64,] domain → codomain) -> TensorMap + +Construct a specific unitary morphism between the codomain and the domain, i.e. return a +`t::TensorMap` where either `scalartype(t) = T` if `T` is a `Number` type or +`storagetype(t) = T` if `T` is a `DenseMatrix` type. If the spaces are not isomorphic, or +the spacetype does not have a Euclidean inner product, an error will be thrown. + +!!! note + There is no canonical choice for a specific unitary, but the current choice is such that + `unitary(cod, dom) == inv(unitary(dom, cod)) = adjoint(unitary(dom, cod))`. + +See also [`isomorphism`](@ref) and [`isometry`](@ref). """ -function unitary(cod::TensorSpace{S}, dom::TensorSpace{S}) where {S} - InnerProductStyle(S) === EuclideanProduct() || throw_invalid_innerproduct(:unitary) - return isomorphism(cod, dom) -end -function unitary(P::TensorMapSpace{S}) where {S} - InnerProductStyle(S) === EuclideanProduct() || throw_invalid_innerproduct(:unitary) - return isomorphism(P) -end -function unitary(A::Type{<:DenseMatrix}, P::TensorMapSpace{S}) where {S} - InnerProductStyle(S) === EuclideanProduct() || throw_invalid_innerproduct(:unitary) - return isomorphism(A, P) -end -function unitary(A::Type{<:DenseMatrix}, cod::TensorSpace{S}, dom::TensorSpace{S}) where {S} +function unitary(::Type{A}, V::TensorMapSpace{S,N₁,N₂}) where {A<:MatOrNumber,S,N₁,N₂} InnerProductStyle(S) === EuclideanProduct() || throw_invalid_innerproduct(:unitary) - return isomorphism(A, cod, dom) + return isomorphism(A, V) end """ - isometry([A::Type{<:DenseMatrix} = Matrix{Float64},] cod::VectorSpace, dom::VectorSpace) - -> TensorMap - -Return a `t::TensorMap` that implements a specific isometry that embeds the domain `dom` -into the codomain `cod`, and which requires that `spacetype(dom)` (`== spacetype(cod)`) has -an Euclidean inner product. An isometry `t` is such that its adjoint `t'` is the left -inverse of `t`, i.e. `t'*t = id(dom)`, while `t*t'` is some idempotent endomorphism of -`cod`, i.e. it squares to itself. When `dom` and `cod` do not allow for such an isometric -inclusion, an error will be thrown. + isometry([T::Type=Float64,] codomain::TensorSpace, domain::TensorSpace) -> TensorMap + isometry([T::Type=Float64,] codomain ← domain) -> TensorMap + isometry([T::Type=Float64,] domain → codomain) -> TensorMap + +Construct a specific isometry between the codomain and the domain, i.e. return a +`t::TensorMap` where either `scalartype(t) = T` if `T` is a `Number` type or +`storagetype(t) = T` if `T` is a `DenseMatrix` type. The isometry `t` then satisfies +`t' * t = id(domain)` and `(t * t')^2 = t * t'`. If the spaces do not allow for such an +isometric inclusion, an error will be thrown. + +See also [`isomorphism`](@ref) and [`unitary`](@ref). """ -isometry(cod::TensorSpace, dom::TensorSpace) = isometry(Matrix{Float64}, cod, dom) -isometry(P::TensorMapSpace) = isometry(codomain(P), domain(P)) -isometry(A::Type{<:DenseMatrix}, P::TensorMapSpace) = isometry(A, codomain(P), domain(P)) -function isometry(A::Type{<:DenseMatrix}, cod::TensorSpace, dom::TensorSpace) - return isometry(A, convert(ProductSpace, cod), convert(ProductSpace, dom)) -end -function isometry(::Type{A}, - cod::ProductSpace{S}, - dom::ProductSpace{S}) where {A<:DenseMatrix,S<:ElementarySpace} +function isometry(::Type{A}, V::TensorMapSpace{S,N₁,N₂}) where {A<:MatOrNumber,S,N₁,N₂} InnerProductStyle(S) === EuclideanProduct() || throw_invalid_innerproduct(:isometry) - dom ≾ cod || - throw(SpaceMismatch("codomain $cod and domain $dom do not allow for an isometric mapping")) - t = TensorMap(s -> A(undef, s), cod, dom) - for (c, b) in blocks(t) + domain(V) ≾ codomain(V) || + throw(SpaceMismatch("$V does not allow for an isometric inclusion")) + t = tensormaptype(S, N₁, N₂, A)(undef, codomain(V), domain(V)) + for (_, b) in blocks(t) MatrixAlgebra.one!(b) end return t end +# expand methods with default arguments +for morphism in (:isomorphism, :unitary, :isometry) + @eval begin + $morphism(V::TensorMapSpace) = $morphism(Float64, V) + $morphism(codomain::TensorSpace, domain::TensorSpace) = $morphism(codomain ← domain) + function $morphism(::Type{T}, codomain::TensorSpace, + domain::TensorSpace) where {T<:MatOrNumber} + return $morphism(T, codomain ← domain) + end + $morphism(t::AbstractTensorMap) = $morphism(storagetype(t), space(t)) + end +end + # Diagonal tensors # ---------------- # TODO: consider adding a specialised DiagonalTensorMap type @@ -185,8 +183,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)))")) @@ -403,8 +401,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))")) @@ -421,8 +418,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))")) @@ -442,19 +438,21 @@ 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 dom = dom1 ⊗ dom2 - t = TensorMap(zeros, promote_type(scalartype(t1), scalartype(t2)), cod, dom) + t = zeros(promote_type(scalartype(t1), scalartype(t2)), cod ← dom) if sectortype(S) === Trivial d1 = dim(cod1) d2 = dim(cod2) diff --git a/src/tensors/tensor.jl b/src/tensors/tensor.jl index 85e46858..679ae3a1 100644 --- a/src/tensors/tensor.jl +++ b/src/tensors/tensor.jl @@ -3,86 +3,109 @@ #==========================================================# #! format: off """ - struct TensorMap{S<:IndexSpace, N₁, N₂, ...} <: AbstractTensorMap{S, N₁, N₂} + struct TensorMap{T, S<:IndexSpace, N₁, N₂, ...} <: AbstractTensorMap{T, 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{T, S<:IndexSpace, N₁, N₂, I<:Sector, A<:Union{<:DenseMatrix{T},SectorDict{I,<:DenseMatrix{T}}}, + F₁, F₂} <: AbstractTensorMap{T, 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, + + # uninitialized constructors + function TensorMap{T,S,N₁,N₂,Trivial,A,Nothing,Nothing}(::UndefInitializer, codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂}) where {T,S<:IndexSpace,N₁,N₂,A<:DenseMatrix{T}} + data = A(undef, dim(codom), dim(dom)) + return TensorMap{T,S,N₁,N₂,Trivial,A,Nothing,Nothing}(data, codom, dom) + end + function TensorMap{T,S,N₁,N₂,I,A,F₁,F₂}(::UndefInitializer, codom::TensorSpace{S}, + dom::TensorSpace{S}) where {T,S<:IndexSpace,N₁,N₂, + I<:Sector,A<:SectorDict{I,<:DenseMatrix{T}},F₁,F₂} + I === sectortype(S) || throw(SectorMismatch()) + blocksectoriterator = blocksectors(codom ← dom) + rowr, rowdims = _buildblockstructure(codom, blocksectoriterator) + colr, coldims = _buildblockstructure(dom, blocksectoriterator) + data = SectorDict(c => valtype(A)(undef, rowdims[c], coldims[c]) for c in blocksectoriterator) + return TensorMap{T,S,N₁,N₂,I,A,F₁,F₂}(data, codom, dom, rowr, colr) + end + + # constructors from data + function TensorMap{T,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}, + {T,S<:IndexSpace, N₁, N₂, I<:Sector, A<:SectorDict{I,<:DenseMatrix{T}}, 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) + return new{T,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{T,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) + {T,S<:IndexSpace,N₁,N₂,A<:DenseMatrix{T}} + T ⊆ field(S) || @warn("scalartype(data) = $T ⊈ $(field(S)))", maxlog = 1) + return new{T,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{T, S, N, I, A, F₁, F₂} = TensorMap{T, 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{T, S, N, I, A, F₁, F₂}` is actually a special case `TensorMap{T, 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{T,S,N,I,A,F₁,F₂} = TensorMap{T,S,N,0,I,A,F₁,F₂} + """ - TrivialTensorMap{S<:IndexSpace, N₁, N₂, A<:DenseMatrix} = TensorMap{S, N₁, N₂, Trivial, + TrivialTensorMap{T, S, N₁, N₂, A<:DenseMatrix} = TensorMap{T, 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{T,S,N₁,N₂,A<:DenseMatrix} = TensorMap{T,S,N₁,N₂,Trivial,A,Nothing, + Nothing} + """ - TrivialTensor{S, N, A} = TrivialTensorMap{S, N, 0, A} + TrivialTensor{T, S, N, A} = TrivialTensorMap{T, 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{T,S,N,A} = TrivialTensorMap{T,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} Return the fully specified type of a tensor map with elementary space `S`, `N₁` output spaces and `N₂` input spaces, either with scalar type `T` or with storage type `T`. """ -function tensormaptype(::Type{S}, N₁::Int, N₂::Int, ::Type{T}) where {S,T} +function tensormaptype(::Type{S}, N₁::Int, N₂::Int, + ::Type{TorA}) where {S,TorA<:MatOrNumber} I = sectortype(S) - if T <: DenseMatrix - M = T - elseif T <: Number - M = Matrix{T} + if TorA <: DenseMatrix + M = TorA + T = scalartype(TorA) + elseif TorA <: Number + M = Matrix{TorA} + T = TorA 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{T,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{T,S,N₁,N₂,I,SectorDict{I,M},F₁,F₂} end end tensormaptype(S, N₁, N₂=0) = tensormaptype(S, N₁, N₂, Float64) @@ -100,12 +123,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{T,S,N₁,N₂,A}}) where + {T,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{T,S,N₁,N₂,I,<:SectorDict{I,A}}}) where + {T,S,N₁,N₂,I<:Sector,A<:DenseMatrix} return A end @@ -134,8 +157,10 @@ Construct a `TensorMap` by explicitly specifying its block data. Alternatively, the domain and codomain can be specified by passing a [`HomSpace`](@ref) using the syntax `codomain ← domain` or `domain → codomain`. """ -function TensorMap(data::AbstractDict{<:Sector,<:DenseMatrix}, codom::ProductSpace{S,N₁}, - dom::ProductSpace{S,N₂}) where {S<:IndexSpace,N₁,N₂} +function TensorMap(data::AbstractDict{<:Sector,<:DenseMatrix}, + V::TensorMapSpace{S,N₁,N₂}) where {S,N₁,N₂} + dom = domain(V) + codom = codomain(V) I = sectortype(S) I == keytype(data) || throw(SectorMismatch()) if I == Trivial @@ -157,43 +182,39 @@ function TensorMap(data::AbstractDict{<:Sector,<:DenseMatrix}, codom::ProductSpa isempty(b) || size(b) == (rowdims[c], coldims[c]) || throw(DimensionMismatch("wrong size of block for sector $c")) end - F₁ = fusiontreetype(I, N₁) - F₂ = fusiontreetype(I, N₂) - 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) + data2 = if isreal(I) + SectorDict(c => data[c] for c in blocksectoriterator) 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) + SectorDict(c => complex(data[c]) for c in blocksectoriterator) end + TT = tensormaptype(S, N₁, N₂, valtype(data)) + return TT(data2, codom, dom, rowr, colr) +end +function TensorMap(data::AbstractDict{<:Sector,<:DenseMatrix}, codom::TensorSpace{S}, + dom::TensorSpace{S}) where {S} + return TensorMap(data, codom ← dom) end -# constructor from general callable that produces block data -function TensorMap(f, codom::ProductSpace{S,N₁}, - dom::ProductSpace{S,N₂}) where {S<:IndexSpace,N₁,N₂} - I = sectortype(S) - if I == Trivial - d1 = dim(codom) - d2 = dim(dom) - data = f((d1, d2)) - A = typeof(data) - return TensorMap{S,N₁,N₂,Trivial,A,Nothing,Nothing}(data, codom, dom) - end - blocksectoriterator = blocksectors(codom ← dom) - rowr, rowdims = _buildblockstructure(codom, blocksectoriterator) - colr, coldims = _buildblockstructure(dom, blocksectoriterator) - if !isreal(I) - data = SectorDict(c => complex(f((rowdims[c], coldims[c]))) - for c in blocksectoriterator) - else - data = SectorDict(c => f((rowdims[c], coldims[c])) for c in blocksectoriterator) - end - 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) +""" + TensorMap{T}(undef, codomain::ProductSpace{S,N₁}, domain::ProductSpace{S,N₂}) + where {T,S,N₁,N₂} + TensorMap{T}(undef, codomain ← domain) + TensorMap{T}(undef, domain → codomain) + # expert mode: select storage type `A` + TensorMap{T,S,N₁,N₂,I,A}(undef, codomain::ProductSpace{S,N₁}, domain::ProductSpace{S,N₂}) + +Construct a `TensorMap` with uninitialized data. +""" +function TensorMap{T}(::UndefInitializer, V::TensorMapSpace{S,N₁,N₂}) where {T,S,N₁,N₂} + TT = tensormaptype(S, N₁, N₂, T) # construct full type + return TT(undef, codomain(V), domain(V)) +end +function TensorMap{T}(::UndefInitializer, codomain::TensorSpace{S}, + domain::TensorSpace{S}) where {T,S} + return TensorMap{T}(undef, codomain ← domain) +end +function Tensor{T}(::UndefInitializer, V::TensorSpace{S}) where {T,S} + return TensorMap{T}(undef, V ← one(V)) end # auxiliary function @@ -219,101 +240,120 @@ function _buildblockstructure(P::ProductSpace{S,N}, blocksectors) where {S<:Inde return treeranges, blockdims end -""" - TensorMap([f, eltype,] codomain::ProductSpace{S,N₁}, domain::ProductSpace{S,N₂}) - where {S<:ElementarySpace,N₁,N₂} - TensorMap([f, eltype,], codomain ← domain) - TensorMap([f, eltype,], domain → codomain) - -Construct a `TensorMap` from a general callable that produces block data for each coupled -sector. - -## Arguments -- `f`: callable object that returns a `DenseMatrix`, or `UndefInitializer`. -- `eltype::Type{<:Number}`: element type of the data. -- `codomain::ProductSpace{S,N₁}`: the codomain as a `ProductSpace` of `N₁` spaces of type - `S<:ElementarySpace`. -- `domain::ProductSpace{S,N₂}`: the domain as a `ProductSpace` of `N₂` spaces of type - `S<:ElementarySpace`. - -If `eltype` is left unspecified, `f` should support the calling syntax `f(::Tuple{Int,Int})` -such that `f((m, n))` returns a `DenseMatrix` with `size(f((m, n))) == (m, n)`. If `eltype` is -specified, `f` is instead called as `f(eltype, (m, n))`. In the case where `f` is left -unspecified or `undef` is passed explicitly, a `TensorMap` with uninitialized data is -generated. +@doc """ + zeros([T=Float64,], codomain::ProductSpace{S,N₁}, domain::ProductSpace{S,N₂}) where {S,N₁,N₂,T} + zeros([T=Float64,], codomain ← domain) -Alternatively, the domain and codomain can be specified by passing a [`HomSpace`](@ref) -using the syntax `codomain ← domain` or `domain → codomain`. +Create a `TensorMap` with element type `T`, of all zeros with spaces specified by `codomain` and `domain`. """ -function TensorMap(f, ::Type{T}, codom::ProductSpace{S}, - dom::ProductSpace{S}) where {S<:IndexSpace,T<:Number} - return TensorMap(d -> f(T, d), codom, dom) -end - -function TensorMap(::Type{T}, codom::ProductSpace{S}, - dom::ProductSpace{S}) where {S<:IndexSpace,T<:Number} - return TensorMap(d -> Array{T}(undef, d), codom, dom) -end +Base.zeros -function TensorMap(::UndefInitializer, ::Type{T}, codom::ProductSpace{S}, - dom::ProductSpace{S}) where {S<:IndexSpace,T<:Number} - return TensorMap(d -> Array{T}(undef, d), codom, dom) -end - -function TensorMap(::UndefInitializer, codom::ProductSpace{S}, - dom::ProductSpace{S}) where {S<:IndexSpace} - return TensorMap(undef, Float64, codom, dom) -end - -function TensorMap(::Type{T}, codom::TensorSpace{S}, - dom::TensorSpace{S}) where {T<:Number,S<:IndexSpace} - return TensorMap(T, convert(ProductSpace, codom), convert(ProductSpace, dom)) -end - -function TensorMap(dataorf, codom::TensorSpace{S}, - dom::TensorSpace{S}) where {S<:IndexSpace} - return TensorMap(dataorf, convert(ProductSpace, codom), convert(ProductSpace, dom)) -end - -function TensorMap(dataorf, ::Type{T}, codom::TensorSpace{S}, - dom::TensorSpace{S}) where {T<:Number,S<:IndexSpace} - return TensorMap(dataorf, T, convert(ProductSpace, codom), convert(ProductSpace, dom)) -end +@doc """ + ones([T=Float64,], codomain::ProductSpace{S,N₁}, domain::ProductSpace{S,N₂}) where {S,N₁,N₂,T} + ones([T=Float64,], codomain ← domain) + +Create a `TensorMap` with element type `T`, of all ones with spaces specified by `codomain` and `domain`. +""" +Base.ones -function TensorMap(codom::TensorSpace{S}, dom::TensorSpace{S}) where {S<:IndexSpace} - return TensorMap(Float64, convert(ProductSpace, codom), convert(ProductSpace, dom)) +for (fname, felt) in ((:zeros, :zero), (:ones, :one)) + @eval begin + function Base.$fname(codomain::TensorSpace{S}, + domain::TensorSpace{S}=one(codomain)) where {S<:IndexSpace} + return Base.$fname(codomain ← domain) + end + function Base.$fname(::Type{T}, codomain::TensorSpace{S}, + domain::TensorSpace{S}=one(codomain)) where {T,S<:IndexSpace} + return Base.$fname(T, codomain ← domain) + end + Base.$fname(V::TensorMapSpace) = Base.$fname(Float64, V) + function Base.$fname(::Type{T}, V::TensorMapSpace) where {T} + t = TensorMap{T}(undef, V) + fill!(t, $felt(T)) + return t + end + end end -function TensorMap(dataorf, T::Type{<:Number}, P::TensorMapSpace{S}) where {S<:IndexSpace} - return TensorMap(dataorf, T, codomain(P), domain(P)) -end +for randfun in (:rand, :randn, :randexp) + randfun! = Symbol(randfun, :!) + _docstr = """ + $randfun([rng=default_rng()], [T=Float64], codomain::ProductSpace{S,N₁}, + domain::ProductSpace{S,N₂}) where {S,N₁,N₂,T} -> t + $randfun([rng=default_rng()], [T=Float64], codomain ← domain) -> t + + Generate a tensor `t` with entries generated by `$randfun)`. + """ + _docstr! = """ + $randfun!([rng=default_rng()], t::AbstractTensorMap) -> t + + Fill the tensor `t` with entries generated by `$randfun!`. + """ + + @eval begin + @doc $_docstr Random.$randfun + @doc $_docstr! Random.$randfun! + + # converting `codomain` and `domain` into `HomSpace` + function Random.$randfun(codomain::TensorSpace{S}, + domain::TensorSpace{S}) where {S<:IndexSpace} + return Random.$randfun(codomain ← domain) + end + function Random.$randfun(::Type{T}, codomain::TensorSpace{S}, + domain::TensorSpace{S}) where {T,S<:IndexSpace} + return Random.$randfun(T, codomain ← domain) + end + function Random.$randfun(rng::Random.AbstractRNG, ::Type{T}, + codomain::TensorSpace{S}, + domain::TensorSpace{S}) where {T,S<:IndexSpace} + return Random.$randfun(rng, T, codomain ← domain) + end -function TensorMap(dataorf, P::TensorMapSpace{S}) where {S<:IndexSpace} - return TensorMap(dataorf, codomain(P), domain(P)) -end + # accepting single `TensorSpace` + Random.$randfun(codomain::TensorSpace) = Random.$randfun(codomain ← one(codomain)) + function Random.$randfun(::Type{T}, codomain::TensorSpace) where {T} + return Random.$randfun(T, codomain ← one(codomain)) + end + function Random.$randfun(rng::Random.AbstractRNG, ::Type{T}, + codomain::TensorSpace) where {T} + return Random.$randfun(rng, T, codomain ← one(domain)) + end -function TensorMap(T::Type{<:Number}, P::TensorMapSpace{S}) where {S<:IndexSpace} - return TensorMap(T, codomain(P), domain(P)) -end + # filling in default eltype + Random.$randfun(V::TensorMapSpace) = Random.$randfun(Float64, V) + function Random.$randfun(rng::Random.AbstractRNG, V::TensorMapSpace) + return Random.$randfun(rng, Float64, V) + end -TensorMap(P::TensorMapSpace{S}) where {S<:IndexSpace} = TensorMap(codomain(P), domain(P)) + # filling in default rng + function Random.$randfun(::Type{T}, V::TensorMapSpace) where {T} + return Random.$randfun(Random.default_rng(), T, V) + end + Random.$randfun!(t::AbstractTensorMap) = Random.$randfun!(Random.default_rng(), t) + + # implementation + function Random.$randfun(rng::Random.AbstractRNG, ::Type{T}, + V::TensorMapSpace) where {T} + t = TensorMap{T}(undef, V) + Random.$randfun!(rng, t) + return t + end -function Tensor(dataorf, T::Type{<:Number}, P::TensorSpace{S}) where {S<:IndexSpace} - return TensorMap(dataorf, T, P, one(P)) + function Random.$randfun!(rng::Random.AbstractRNG, t::AbstractTensorMap) + for (_, b) in blocks(t) + Random.$randfun!(rng, b) + end + return t + end + end end -Tensor(dataorf, P::TensorSpace{S}) where {S<:IndexSpace} = TensorMap(dataorf, P, one(P)) - -Tensor(T::Type{<:Number}, P::TensorSpace{S}) where {S<:IndexSpace} = TensorMap(T, P, one(P)) - -Tensor(P::TensorSpace{S}) where {S<:IndexSpace} = TensorMap(P, one(P)) - # constructor starting from a dense array """ TensorMap(data::DenseArray, codomain::ProductSpace{S,N₁}, domain::ProductSpace{S,N₂}; tol=sqrt(eps(real(float(eltype(data)))))) where {S<:ElementarySpace,N₁,N₂} - TensorMap(data, codomain ← domain) - TensorMap(data, domain → codomain) + TensorMap(data, codomain ← domain; tol=sqrt(eps(real(float(eltype(data)))))) + TensorMap(data, domain → codomain; tol=sqrt(eps(real(float(eltype(data)))))) Construct a `TensorMap` from a plain multidimensional array. @@ -343,8 +383,10 @@ using the syntax `codomain ← domain` or `domain → codomain`. array is possible, and only in the case where the `data` actually respects the specified symmetry structure. """ -function TensorMap(data::DenseArray, codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂}; +function TensorMap(data::DenseArray, V::TensorMapSpace{S,N₁,N₂}; tol=sqrt(eps(real(float(eltype(data)))))) where {S<:IndexSpace,N₁,N₂} + codom = codomain(V) + dom = domain(V) (d1, d2) = (dim(codom), dim(dom)) if !(length(data) == d1 * d2 || size(data) == (d1, d2) || size(data) == (dims(codom)..., dims(dom)...)) @@ -354,10 +396,10 @@ 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) + return tensormaptype(S, N₁, N₂, A)(data2, codom, dom) end - t = TensorMap(undef, eltype(data), codom, dom) + t = TensorMap{eltype(data)}(undef, codom, dom) project_symmetric!(t, data) if !isapprox(data, convert(Array, t); atol=tol) @@ -398,52 +440,34 @@ function project_symmetric!(t::TensorMap, data::DenseArray) return t end +function TensorMap(data::DenseArray, codom::TensorSpace{S}, + dom::TensorSpace{S}; kwargs...) where {S} + return TensorMap(data, codom ← dom; kwargs...) +end # Efficient copy constructors #----------------------------- Base.copy(t::TrivialTensorMap) = typeof(t)(copy(t.data), t.codom, t.dom) Base.copy(t::TensorMap) = typeof(t)(deepcopy(t.data), t.codom, t.dom, t.rowr, t.colr) -# Similar -#--------- -# 4 arguments -function Base.similar(t::AbstractTensorMap, T::Type, codomain::VectorSpace, - domain::VectorSpace) - return similar(t, T, codomain ← domain) -end -# 3 arguments -function Base.similar(t::AbstractTensorMap, codomain::VectorSpace, domain::VectorSpace) - return similar(t, scalartype(t), codomain ← domain) -end -function Base.similar(t::AbstractTensorMap, T::Type, codomain::VectorSpace) - return similar(t, T, codomain ← one(codomain)) -end -# 2 arguments -function Base.similar(t::AbstractTensorMap, codomain::VectorSpace) - return similar(t, scalartype(t), codomain ← one(codomain)) -end -Base.similar(t::AbstractTensorMap, P::TensorMapSpace) = similar(t, scalartype(t), P) -Base.similar(t::AbstractTensorMap, T::Type) = similar(t, T, space(t)) -# 1 argument -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} +# specializations when data can be re-used +function Base.similar(t::TensorMap, ::Type{TorA}, + P::TensorMapSpace{S}) where {TorA<:MatOrNumber,S} N₁ = length(codomain(P)) N₂ = length(domain(P)) I = sectortype(S) + # speed up specialized cases - 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)) - end - F₁ = fusiontreetype(I, N₁) - F₂ = fusiontreetype(I, N₂) + TT = tensormaptype(S, N₁, N₂, TorA) + I === Trivial && return TT(undef, codomain(P), domain(P)) + 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) + data = if TorA <: Number + SectorDict(c => similar(b, TorA) for (c, b) in blocks(t)) + else + SectorDict(c => similar(TorA, size(b)) for (c, b) in blocks(t)) + end + return TT(data, codomain(P), domain(P), t.rowr, t.colr) end blocksectoriterator = blocksectors(P) @@ -465,7 +489,7 @@ function Base.similar(t::TensorMap{S}, ::Type{T}, P::TensorMapSpace{S}) where {T else rowr, rowdims = _buildblockstructure(codomain(P), blocksectoriterator) end - # try to recylce colr + # try to recycle colr if domain(P) == codomain(t) && all(c -> haskey(t.rowr, c), blocksectoriterator) if length(t.rowr) == length(blocksectoriterator) colr = t.rowr @@ -483,11 +507,10 @@ function Base.similar(t::TensorMap{S}, ::Type{T}, P::TensorMapSpace{S}) where {T else colr, coldims = _buildblockstructure(domain(P), blocksectoriterator) end - M = similarstoragetype(t, T) + M = storagetype(TT) 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 TT(data, codomain(P), domain(P), rowr, colr) end function Base.complex(t::AbstractTensorMap) @@ -541,7 +564,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 @@ -550,7 +573,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] @@ -564,8 +587,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)) @@ -587,9 +612,9 @@ end end """ - Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,I}, + Base.getindex(t::TensorMap{T,S,N₁,N₂,I}, f₁::FusionTree{I,N₁}, - f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector} + f₂::FusionTree{I,N₂}) where {T,SN₁,N₂,I<:Sector} -> StridedViews.StridedView t[f₁, f₂] @@ -600,9 +625,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{T,S,N₁,N₂,I}, f₁::FusionTree{I,N₁}, - f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector} + f₂::FusionTree{I,N₂}) where {T,S,N₁,N₂,I<:Sector} c = f₁.coupled @boundscheck begin c == f₂.coupled || throw(SectorMismatch()) @@ -616,10 +641,10 @@ column indices correspond to `f₂.uncoupled`. end """ - Base.setindex!(t::TensorMap{<:IndexSpace,N₁,N₂,I}, + Base.setindex!(t::TensorMap{T,S,N₁,N₂,I}, v, f₁::FusionTree{I,N₁}, - f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector} + f₂::FusionTree{I,N₂}) where {T,S,N₁,N₂,I<:Sector} t[f₁, f₂] = v Copies `v` into the data slice of `t` corresponding to the splitting - fusion tree pair @@ -627,12 +652,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{T,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{T,S,N₁,N₂,I}, v, f₁::FusionTree{I,N₁}, - f₂::FusionTree{I,N₂}) where {N₁,N₂,I<:Sector} + f₂::FusionTree{I,N₂}) where {T,S,N₁,N₂, + I<:Sector} return copy!(getindex(t, f₁, f₂), v) end @@ -685,16 +711,16 @@ end function Base.summary(io::IO, t::TensorMap) return print(io, "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₂]) @@ -711,28 +737,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,22 +767,23 @@ end #--------------------------- Base.convert(::Type{TensorMap}, t::TensorMap) = t function Base.convert(::Type{TensorMap}, t::AbstractTensorMap) - return copy!(TensorMap(undef, scalartype(t), codomain(t), domain(t)), t) + return copy!(TensorMap{scalartype(t)}(undef, space(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(TT::Type{<:TensorMap{T,S,N₁,N₂}}, + t::AbstractTensorMap{<:Any,S,N₁,N₂}) where {T,S,N₁,N₂} + if typeof(t) === TT return t else - data = Dict{I,storagetype(T)}(c => convert(storagetype(T), b) - for (c, b) in blocks(t)) + data = Dict{sectortype(TT),storagetype(TT)}(c => convert(storagetype(TT), 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{<:TT₁}, + ::Type{<:TT₂}) where {S,N₁,N₂, + TT₁<:TensorMap{<:Any,S,N₁,N₂}, + TT₂<:TensorMap{<:Any,S,N₁,N₂}} + return tensormaptype(S, N₁, N₂, promote_type(storagetype(TT₁), storagetype(TT₂))) end diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl index 1d3c7f30..9b43d590 100644 --- a/src/tensors/tensoroperations.jl +++ b/src/tensors/tensoroperations.jl @@ -1,21 +1,37 @@ # Implement full TensorOperations.jl interface #---------------------------------------------- TO.tensorstructure(t::AbstractTensorMap) = space(t) -function TO.tensorstructure(t::AbstractTensorMap, iA::Int, conjA::Symbol) - return conjA == :N ? space(t, iA) : conj(space(t, iA)) +function TO.tensorstructure(t::AbstractTensorMap, iA::Int, conjA::Bool) + return !conjA ? space(t, iA) : conj(space(t, iA)) end -function TO.tensoralloc(::Type{TT}, structure, istemp=false, - backend::Backend...) where {TT<:AbstractTensorMap} - function blockallocator(d) - return TO.tensoralloc(storagetype(TT), d, istemp, backend...) - end - return TensorMap(blockallocator, structure) +function TO.tensoralloc(::Type{TT}, structure::TensorMapSpace{S,N₁,N₂}, istemp::Val, + allocator=TO.DefaultAllocator()) where {T,S,N₁,N₂,A, + TT<:TrivialTensorMap{T,S,N₁, + N₂, + A}} + data = TO.tensoralloc(A, (dim(codomain(structure)), dim(domain(structure))), istemp, + allocator) + return TT(data, codomain(structure), domain(structure)) +end + +function TO.tensoralloc(::Type{TT}, structure::TensorMapSpace{S,N₁,N₂}, istemp::Val, + allocator=TO.DefaultAllocator()) where {T,S,N₁,N₂, + TT<:AbstractTensorMap{T,S, + N₁, + N₂}} + blocksectoriterator = blocksectors(structure) + rowr, rowdims = _buildblockstructure(codomain(structure), blocksectoriterator) + colr, coldims = _buildblockstructure(domain(structure), blocksectoriterator) + A = storagetype(TT) + blockallocator(c) = TO.tensoralloc(A, (rowdims[c], coldims[c]), istemp, allocator) + data = SectorDict(c => blockallocator(c) for c in blocksectoriterator) + return TT(data, codomain(structure), domain(structure), rowr, colr) end -function TO.tensorfree!(t::AbstractTensorMap, backend::Backend...) +function TO.tensorfree!(t::AbstractTensorMap, allocator=TO.DefaultAllocator()) for (c, b) in blocks(t) - TO.tensorfree!(b, backend...) + TO.tensorfree!(b, allocator) end return nothing end @@ -34,106 +50,99 @@ 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} - if conjA == :N +function TO.tensoradd!(C::AbstractTensorMap, + A::AbstractTensorMap, pA::Index2Tuple, conjA::Bool, + α::Number, β::Number, backend::AbstractBackend, allocator) + if !conjA A′ = A - pC′ = _canonicalize(pC, C) - elseif conjA == :C - A′ = adjoint(A) - pC′ = adjointtensorindices(A, _canonicalize(pC, C)) + pA′ = _canonicalize(pA, C) else - throw(ArgumentError("unknown conjugation flag $conjA")) + A′ = adjoint(A) + pA′ = adjointtensorindices(A, _canonicalize(pA, C)) end - add_permute!(C, A′, pC′, α, β, backend...) + add_permute!(C, A′, pA′, α, β, backend) return C end -function TO.tensoradd_type(TC, ::Index2Tuple{N₁,N₂}, A::AbstractTensorMap{S}, - ::Symbol) where {S,N₁,N₂} +function TO.tensoradd_type(TC, A::AbstractTensorMap, ::Index2Tuple{N₁,N₂}, + ::Bool) 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₂} - if conjA == :N - return permute(space(A), pC) +function TO.tensoradd_structure(A::AbstractTensorMap, pA::Index2Tuple{N₁,N₂}, + conjA::Bool) where {N₁,N₂} + if !conjA + return permute(space(A), pA) else - return TO.tensoradd_structure(adjointtensorindices(A, pC), adjoint(A), :N) + return TO.tensoradd_structure(adjoint(A), adjointtensorindices(A, pA), false) end end # tensortrace! -function TO.tensortrace!(C::AbstractTensorMap{S}, p::Index2Tuple, - A::AbstractTensorMap{S}, q::Index2Tuple, conjA::Symbol, - α::Number, β::Number, backend::Backend...) where {S} - if conjA == :N +function TO.tensortrace!(C::AbstractTensorMap, + A::AbstractTensorMap, p::Index2Tuple, q::Index2Tuple, + conjA::Bool, + α::Number, β::Number, backend::AbstractBackend, allocator) + if !conjA A′ = A p′ = _canonicalize(p, C) q′ = q - elseif conjA == :C + else A′ = adjoint(A) p′ = adjointtensorindices(A, _canonicalize(p, C)) q′ = adjointtensorindices(A, q) - else - throw(ArgumentError("unknown conjugation flag $conjA")) end - # TODO: novel syntax for tensortrace? - # tensortrace!(C, pC′, A′, qA′, α, β, backend...) - trace_permute!(C, A′, p′, q′, α, β, backend...) + trace_permute!(C, A′, p′, q′, α, β, backend) return C 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, + A::AbstractTensorMap, pA::Index2Tuple, conjA::Bool, + B::AbstractTensorMap, pB::Index2Tuple, conjB::Bool, + pAB::Index2Tuple, α::Number, β::Number, + backend::AbstractBackend, allocator) pAB′ = _canonicalize(pAB, C) - if conjA == :N + if !conjA A′ = A pA′ = pA - elseif conjA == :C + else A′ = A' pA′ = adjointtensorindices(A, pA) - else - throw(ArgumentError("unknown conjugation flag $conjA")) end - if conjB == :N + if !conjB B′ = B pB′ = pB - elseif conjB == :C + else B′ = B' pB′ = adjointtensorindices(B, pB) - else - throw(ArgumentError("unknown conjugation flag $conjB")) end - contract!(C, A′, pA′, B′, pB′, pAB′, α, β, backend...) + contract!(C, A′, pA′, B′, pB′, pAB′, α, β, backend, allocator) return C end -function TO.tensorcontract_type(TC, ::Index2Tuple{N₁,N₂}, - A::AbstractTensorMap{S}, pA, conjA, - B::AbstractTensorMap{S}, pB, conjB) where {S,N₁,N₂} +function TO.tensorcontract_type(TC, + A::AbstractTensorMap, ::Index2Tuple, ::Bool, + B::AbstractTensorMap, ::Index2Tuple, ::Bool, + ::Index2Tuple{N₁,N₂}) 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₂} - sA = TO.tensoradd_structure(pA, A, conjA) - sB = TO.tensoradd_structure(pB, B, conjB) - return permute(compose(sA, sB), pC) +function TO.tensorcontract_structure(A::AbstractTensorMap, pA::Index2Tuple, conjA::Bool, + B::AbstractTensorMap, pB::Index2Tuple, conjB::Bool, + pAB::Index2Tuple{N₁,N₂}) where {N₁,N₂} + sA = TO.tensoradd_structure(A, pA, conjA) + sB = TO.tensoradd_structure(B, pB, conjB) + return permute(compose(sA, sB), pAB) 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::Bool, + tB::AbstractTensorMap, iB::Int, conjB::Bool, + label) sA = TO.tensorstructure(tA, iA, conjA)' sB = TO.tensorstructure(tB, iB, conjB) sA == sB || @@ -149,16 +158,24 @@ 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::AbstractBackend=TO.DefaultBackend()) + # 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 space(tdst) == permute(space(tsrc), (p₁, p₂)) || throw(SpaceMismatch("trace: tsrc = $(codomain(tsrc))←$(domain(tsrc)), @@ -174,7 +191,7 @@ function trace_permute!(tdst::AbstractTensorMap{S,N₁,N₂}, cod = codomain(tsrc) dom = domain(tsrc) n = length(cod) - TO.tensortrace!(tdst[], (p₁, p₂), tsrc[], (q₁, q₂), :N, α, β) + TO.tensortrace!(tdst[], tsrc[], (p₁, p₂), (q₁, q₂), false, α, β, backend) # elseif FusionStyle(I) isa UniqueFusion else cod = codomain(tsrc) @@ -197,7 +214,7 @@ function trace_permute!(tdst::AbstractTensorMap{S,N₁,N₂}, C = tdst[f₁′′, f₂′′] A = tsrc[f₁, f₂] α′ = α * coeff - TO.tensortrace!(C, (p₁, p₂), A, (q₁, q₂), :N, α′, One(), backend...) + TO.tensortrace!(C, A, (p₁, p₂), (q₁, q₂), false, α′, One(), backend) end end end @@ -209,21 +226,21 @@ 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::AbstractBackend, allocator) + 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 - ipC = TupleTools.invperm((p₁..., p₂...)) - oindAinC = TupleTools.getindices(ipC, ntuple(n -> n, N₁)) - oindBinC = TupleTools.getindices(ipC, ntuple(n -> n + N₁, N₂)) + ipAB = TupleTools.invperm((p₁..., p₂...)) + oindAinC = TupleTools.getindices(ipAB, ntuple(n -> n, N₁)) + oindBinC = TupleTools.getindices(ipAB, ntuple(n -> n + N₁, N₂)) qA = TupleTools.sortperm(cindA) cindA′ = TupleTools.getindices(cindA, qA) @@ -267,16 +284,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 @@ -286,9 +304,9 @@ function _contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S}, A′ = permute(A, (oindA, cindA); copy=copyA) B′ = permute(B, (cindB, oindB)) A′ = twist!(A′, filter(i -> !isdual(space(A′, i)), domainind(A′))) - ipC = TupleTools.invperm((p₁..., p₂...)) - oindAinC = TupleTools.getindices(ipC, ntuple(n -> n, N₁)) - oindBinC = TupleTools.getindices(ipC, ntuple(n -> n + N₁, N₂)) + ipAB = TupleTools.invperm((p₁..., p₂...)) + oindAinC = TupleTools.getindices(ipAB, ntuple(n -> n, N₁)) + oindBinC = TupleTools.getindices(ipAB, ntuple(n -> n + N₁, N₂)) if has_shared_permute(C, (oindAinC, oindBinC)) C′ = permute(C, (oindAinC, oindBinC)) mul!(C′, A′, B′, α, β) @@ -301,7 +319,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/ad.jl b/test/ad.jl index 020a49e1..16fa05c7 100644 --- a/test/ad.jl +++ b/test/ad.jl @@ -13,7 +13,7 @@ end # Test utility # ------------- function ChainRulesTestUtils.rand_tangent(rng::AbstractRNG, x::AbstractTensorMap) - return TensorMap(randn, scalartype(x), space(x)) + return randn!(similar(x)) end ChainRulesTestUtils.rand_tangent(::AbstractRNG, ::VectorSpace) = NoTangent() function ChainRulesTestUtils.test_approx(actual::AbstractTensorMap, @@ -44,13 +44,6 @@ function FiniteDifferences.to_vec(t::TensorKit.SectorDict) return vec_real, from_vec end -function _randomize!(a::TensorMap) - for b in values(blocks(a)) - copyto!(b, randn(size(b))) - end - return a -end - # Float32 and finite differences don't mix well precision(::Type{<:Union{Float32,Complex{Float32}}}) = 1e-2 precision(::Type{<:Union{Float64,Complex{Float64}}}) = 1e-6 @@ -134,8 +127,8 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), @timedtestset "Automatic Differentiation with spacetype $(TensorKit.type_repr(eltype(V)))" verbose = true for V in Vlist @timedtestset "Basic utility" begin - T1 = TensorMap(randn, Float64, V[1] ⊗ V[2] ← V[3] ⊗ V[4]) - T2 = TensorMap(randn, ComplexF64, V[1] ⊗ V[2] ← V[3] ⊗ V[4]) + T1 = randn(Float64, V[1] ⊗ V[2] ← V[3] ⊗ V[4]) + T2 = randn(ComplexF64, V[1] ⊗ V[2] ← V[3] ⊗ V[4]) P1 = ProjectTo(T1) @test P1(T1) == T1 @@ -150,8 +143,8 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), end @timedtestset "Basic Linear Algebra with scalartype $T" for T in (Float64, ComplexF64) - A = TensorMap(randn, T, V[1] ⊗ V[2] ← V[3] ⊗ V[4] ⊗ V[5]) - B = TensorMap(randn, T, space(A)) + A = randn(T, V[1] ⊗ V[2] ← V[3] ⊗ V[4] ⊗ V[5]) + B = randn(T, space(A)) test_rrule(+, A, B) test_rrule(-, A) @@ -161,27 +154,27 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), test_rrule(*, α, A) test_rrule(*, A, α) - C = TensorMap(randn, T, domain(A), codomain(A)) + C = randn(T, domain(A), codomain(A)) test_rrule(*, A, C) test_rrule(permute, A, ((1, 3, 2), (5, 4))) - D = TensorMap(randn, T, V[1] ⊗ V[2] ← V[3]) - E = TensorMap(randn, T, V[4] ← V[5]) + D = randn(T, V[1] ⊗ V[2] ← V[3]) + E = randn(T, V[4] ← V[5]) test_rrule(⊗, D, E) end @timedtestset "Linear Algebra part II with scalartype $T" for T in (Float64, ComplexF64) for i in 1:3 - E = TensorMap(randn, T, ⊗(V[1:i]...) ← ⊗(V[1:i]...)) + E = randn(T, ⊗(V[1:i]...) ← ⊗(V[1:i]...)) test_rrule(LinearAlgebra.tr, E) end - A = TensorMap(randn, T, V[1] ⊗ V[2] ← V[3] ⊗ V[4] ⊗ V[5]) + A = randn(T, V[1] ⊗ V[2] ← V[3] ⊗ V[4] ⊗ V[5]) test_rrule(LinearAlgebra.adjoint, A) test_rrule(LinearAlgebra.norm, A, 2) - B = TensorMap(randn, T, space(A)) + B = randn(T, space(A)) test_rrule(LinearAlgebra.dot, A, B) end @@ -200,19 +193,19 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), p = _repartition(_p, rand(0:k1)) q = _repartition(_q, k2) ip = _repartition(invperm(linearize((_p, _q))), rand(0:(k1 + 2 * k2))) - A = TensorMap(randn, T, permute(prod(V1) ⊗ prod(V2) ← prod(V2), ip)) + A = randn(T, permute(prod(V1) ⊗ prod(V2) ← prod(V2), ip)) α = randn(T) β = randn(T) - for conjA in (:N, :C) - C = _randomize!(TensorOperations.tensoralloc_add(T, p, A, conjA, false)) - test_rrule(tensortrace!, C, p, A, q, conjA, α, β; atol, rtol) + for conjA in (false, true) + C = randn!(TensorOperations.tensoralloc_add(T, A, p, conjA, Val(false))) + test_rrule(tensortrace!, C, A, p, q, conjA, α, β; atol, rtol) end end end @timedtestset "tensoradd!" begin - A = TensorMap(randn, T, V[1] ⊗ V[2] ⊗ V[3] ← V[4] ⊗ V[5]) + A = randn(T, V[1] ⊗ V[2] ⊗ V[3] ← V[4] ⊗ V[5]) α = randn(T) β = randn(T) @@ -220,11 +213,11 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), for _ in 1:5 p = randindextuple(length(V)) - C1 = _randomize!(TensorOperations.tensoralloc_add(T, p, A, :N, false)) - test_rrule(tensoradd!, C1, p, A, :N, α, β; atol, rtol) + C1 = randn!(TensorOperations.tensoralloc_add(T, A, p, false, Val(false))) + test_rrule(tensoradd!, C1, A, p, false, α, β; atol, rtol) - C2 = _randomize!(TensorOperations.tensoralloc_add(T, p, A, :C, false)) - test_rrule(tensoradd!, C2, p, A, :C, α, β; atol, rtol) + C2 = randn!(TensorOperations.tensoralloc_add(T, A, p, true, Val(false))) + test_rrule(tensoradd!, C2, A, p, true, α, β; atol, rtol) A = rand(Bool) ? C1 : C2 end @@ -255,33 +248,31 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), β = randn(T) V2_conj = prod(conj, V2; init=one(V[1])) - for conjA in (:N, :C), conjB in (:N, :C) - A = TensorMap(randn, T, - permute(V1 ← (conjA === :C ? V2_conj : V2), ipA)) - B = TensorMap(randn, T, - permute((conjB === :C ? V2_conj : V2) ← V3, ipB)) - C = _randomize!(TensorOperations.tensoralloc_contract(T, pAB, A, pA, - conjA, - B, pB, conjB, - false)) - test_rrule(tensorcontract!, C, pAB, - A, pA, conjA, B, pB, conjB, + for conjA in (false, true), conjB in (false, true) + A = randn(T, permute(V1 ← (conjA ? V2_conj : V2), ipA)) + B = randn(T, permute((conjB ? V2_conj : V2) ← V3, ipB)) + C = randn!(TensorOperations.tensoralloc_contract(T, A, pA, + conjA, + B, pB, conjB, pAB, + Val(false))) + test_rrule(tensorcontract!, C, + A, pA, conjA, B, pB, conjB, pAB, α, β; atol, rtol) end end end @timedtestset "tensorscalar" begin - A = Tensor(randn, T, ProductSpace{typeof(V[1]),0}()) + A = randn(T, ProductSpace{typeof(V[1]),0}()) test_rrule(tensorscalar, A) end end @timedtestset "Factorizations with scalartype $T" for T in (Float64, ComplexF64) - A = TensorMap(randn, T, V[1] ⊗ V[2] ← V[3] ⊗ V[4] ⊗ V[5]) - B = TensorMap(randn, T, space(A)') - C = TensorMap(randn, T, V[1] ⊗ V[2] ← V[1] ⊗ V[2]) - H = TensorMap(randn, T, V[3] ⊗ V[4] ← V[3] ⊗ V[4]) + A = randn(T, V[1] ⊗ V[2] ← V[3] ⊗ V[4] ⊗ V[5]) + B = randn(T, space(A)') + C = randn(T, V[1] ⊗ V[2] ← V[1] ⊗ V[2]) + H = randn(T, V[3] ⊗ V[4] ← V[3] ⊗ V[4]) H = (H + H') / 2 atol = precision(T) @@ -298,8 +289,8 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), end let (D, V) = eig(C) - ΔD = TensorMap(randn, scalartype(D), space(D)) - ΔV = TensorMap(randn, scalartype(V), space(V)) + ΔD = randn(scalartype(D), space(D)) + ΔV = randn(scalartype(V), space(V)) gaugepart = V' * ΔV for (c, b) in blocks(gaugepart) mul!(block(ΔV, c), inv(block(V, c))', Diagonal(diag(b)), -1, 1) @@ -308,8 +299,8 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), end let (D, U) = eigh′(H) - ΔD = TensorMap(randn, scalartype(D), space(D)) - ΔU = TensorMap(randn, scalartype(U), space(U)) + ΔD = randn(scalartype(D), space(D)) + ΔU = randn(scalartype(U), space(U)) if T <: Complex gaugepart = U' * ΔU for (c, b) in blocks(gaugepart) @@ -320,9 +311,9 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), end let (U, S, V, ϵ) = tsvd(A) - ΔU = TensorMap(randn, scalartype(U), space(U)) - ΔS = TensorMap(randn, scalartype(S), space(S)) - ΔV = TensorMap(randn, scalartype(V), space(V)) + ΔU = randn(scalartype(U), space(U)) + ΔS = randn(scalartype(S), space(S)) + ΔV = randn(scalartype(V), space(V)) if T <: Complex # remove gauge dependent components gaugepart = U' * ΔU + V * ΔV' for (c, b) in blocks(gaugepart) @@ -334,18 +325,18 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), allS = mapreduce(x -> diag(x[2]), vcat, blocks(S)) truncval = (maximum(allS) + minimum(allS)) / 2 U, S, V, ϵ = tsvd(A; trunc=truncerr(truncval)) - ΔU = TensorMap(randn, scalartype(U), space(U)) - ΔS = TensorMap(randn, scalartype(S), space(S)) - ΔV = TensorMap(randn, scalartype(V), space(V)) + ΔU = randn(scalartype(U), space(U)) + ΔS = randn(scalartype(S), space(S)) + ΔV = randn(scalartype(V), space(V)) T <: Complex && remove_svdgauge_depence!(ΔU, ΔV, U, S, V) test_rrule(tsvd, A; atol, output_tangent=(ΔU, ΔS, ΔV, 0.0), fkwargs=(; trunc=truncerr(truncval))) end let (U, S, V, ϵ) = tsvd(B) - ΔU = TensorMap(randn, scalartype(U), space(U)) - ΔS = TensorMap(randn, scalartype(S), space(S)) - ΔV = TensorMap(randn, scalartype(V), space(V)) + ΔU = randn(scalartype(U), space(U)) + ΔS = randn(scalartype(S), space(S)) + ΔV = randn(scalartype(V), space(V)) T <: Complex && remove_svdgauge_depence!(ΔU, ΔV, U, S, V) test_rrule(tsvd, B; atol, output_tangent=(ΔU, ΔS, ΔV, 0.0)) @@ -353,40 +344,40 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), for (c, b) in blocks(S))) U, S, V, ϵ = tsvd(B; trunc=truncspace(Vtrunc)) - ΔU = TensorMap(randn, scalartype(U), space(U)) - ΔS = TensorMap(randn, scalartype(S), space(S)) - ΔV = TensorMap(randn, scalartype(V), space(V)) + ΔU = randn(scalartype(U), space(U)) + ΔS = randn(scalartype(S), space(S)) + ΔV = randn(scalartype(V), space(V)) T <: Complex && remove_svdgauge_depence!(ΔU, ΔV, U, S, V) test_rrule(tsvd, B; atol, output_tangent=(ΔU, ΔS, ΔV, 0.0), fkwargs=(; trunc=truncspace(Vtrunc))) end let (U, S, V, ϵ) = tsvd(C) - ΔU = TensorMap(randn, scalartype(U), space(U)) - ΔS = TensorMap(randn, scalartype(S), space(S)) - ΔV = TensorMap(randn, scalartype(V), space(V)) + ΔU = randn(scalartype(U), space(U)) + ΔS = randn(scalartype(S), space(S)) + ΔV = randn(scalartype(V), space(V)) T <: Complex && remove_svdgauge_depence!(ΔU, ΔV, U, S, V) test_rrule(tsvd, C; atol, output_tangent=(ΔU, ΔS, ΔV, 0.0)) c, = TensorKit.MatrixAlgebra._argmax(x -> sqrt(dim(x[1])) * maximum(diag(x[2])), blocks(S)) U, S, V, ϵ = tsvd(C; trunc=truncdim(2 * dim(c))) - ΔU = TensorMap(randn, scalartype(U), space(U)) - ΔS = TensorMap(randn, scalartype(S), space(S)) - ΔV = TensorMap(randn, scalartype(V), space(V)) + ΔU = randn(scalartype(U), space(U)) + ΔS = randn(scalartype(S), space(S)) + ΔV = randn(scalartype(V), space(V)) T <: Complex && remove_svdgauge_depence!(ΔU, ΔV, U, S, V) test_rrule(tsvd, C; atol, output_tangent=(ΔU, ΔS, ΔV, 0.0), fkwargs=(; trunc=truncdim(2 * dim(c)))) end let D = LinearAlgebra.eigvals(C) - ΔD = diag(TensorMap(randn, complex(scalartype(C)), space(C))) + ΔD = diag(randn(complex(scalartype(C)), space(C))) test_rrule(LinearAlgebra.eigvals, C; atol, output_tangent=ΔD, fkwargs=(; sortby=nothing)) end let S = LinearAlgebra.svdvals(C) - ΔS = diag(TensorMap(randn, real(scalartype(C)), space(C))) + ΔS = diag(randn(real(scalartype(C)), space(C))) test_rrule(LinearAlgebra.svdvals, C; atol, output_tangent=ΔS) end end diff --git a/test/braidingtensor.jl b/test/braidingtensor.jl index b61abc61..f290bc7f 100644 --- a/test/braidingtensor.jl +++ b/test/braidingtensor.jl @@ -11,121 +11,121 @@ V3 = GradedSpace{IsingAnyon}(:I => 2, :psi => 1, :sigma => 1) for V in (V1, V2, V3) @show V - t = TensorMap(randn, V * V' * V' * V, V * V') + t = 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/bugfixes.jl b/test/bugfixes.jl index d395235f..d539c01e 100644 --- a/test/bugfixes.jl +++ b/test/bugfixes.jl @@ -1,23 +1,23 @@ @timedtestset "Bugfixes" verbose = true begin @testset "BugfixConvert" begin - v = TensorMap(randn, ComplexF64, - (Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((-3, 1 / 2, 1) => 3, - (-5, 1 / 2, 1) => 10, - (-7, 1 / 2, 1) => 13, - (-9, 1 / 2, 1) => 9, - (-11, 1 / 2, 1) => 1, - (-5, 3 / 2, 1) => 3, - (-7, 3 / 2, 1) => 3, - (-9, 3 / 2, 1) => 1) ⊗ - Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((1, 1 / 2, 1) => 1)') ← - Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((-3, 1 / 2, 1) => 3, - (-5, 1 / 2, 1) => 10, - (-7, 1 / 2, 1) => 13, - (-9, 1 / 2, 1) => 9, - (-11, 1 / 2, 1) => 1, - (-5, 3 / 2, 1) => 3, - (-7, 3 / 2, 1) => 3, - (-9, 3 / 2, 1) => 1)) + v = randn(ComplexF64, + (Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((-3, 1 / 2, 1) => 3, + (-5, 1 / 2, 1) => 10, + (-7, 1 / 2, 1) => 13, + (-9, 1 / 2, 1) => 9, + (-11, 1 / 2, 1) => 1, + (-5, 3 / 2, 1) => 3, + (-7, 3 / 2, 1) => 3, + (-9, 3 / 2, 1) => 1) ⊗ + Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((1, 1 / 2, 1) => 1)') ← + Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((-3, 1 / 2, 1) => 3, + (-5, 1 / 2, 1) => 10, + (-7, 1 / 2, 1) => 13, + (-9, 1 / 2, 1) => 9, + (-11, 1 / 2, 1) => 1, + (-5, 3 / 2, 1) => 3, + (-7, 3 / 2, 1) => 3, + (-9, 3 / 2, 1) => 1)) w = convert(typeof(real(v)), v) @test w == v @test scalartype(w) == Float64 diff --git a/test/planar.jl b/test/planar.jl index d7edf146..e071ec6f 100644 --- a/test/planar.jl +++ b/test/planar.jl @@ -12,15 +12,17 @@ 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}) - tdst = TensorMap(undef, scalartype(tsrc), - force_planar(codomain(tsrc)) ← force_planar(domain(tsrc))) +function force_planar(tsrc::TensorMap{<:Any,ComplexSpace}) + tdst = TensorMap{scalartype(tsrc)}(undef, + force_planar(codomain(tsrc)) ← + force_planar(domain(tsrc))) copyto!(blocks(tdst)[PlanarTrivial()], blocks(tsrc)[Trivial()]) return tdst end -function force_planar(tsrc::TensorMap{<:GradedSpace}) - tdst = TensorMap(undef, scalartype(tsrc), - force_planar(codomain(tsrc)) ← force_planar(domain(tsrc))) +function force_planar(tsrc::TensorMap{<:Any,<:GradedSpace}) + tdst = TensorMap{scalartype(tsrc)}(undef, + force_planar(codomain(tsrc)) ← + force_planar(domain(tsrc))) for (c, b) in blocks(tsrc) copyto!(blocks(tdst)[c ⊠ PlanarTrivial()], b) end @@ -29,32 +31,32 @@ end @testset "planar methods" verbose = true begin @testset "planaradd" begin - A = TensorMap(randn, ℂ^2 ⊗ ℂ^3 ← ℂ^6 ⊗ ℂ^5 ⊗ ℂ^4) - C = TensorMap(randn, (ℂ^5)' ⊗ (ℂ^6)' ← ℂ^4 ⊗ (ℂ^3)' ⊗ (ℂ^2)') + A = randn(ℂ^2 ⊗ ℂ^3 ← ℂ^6 ⊗ ℂ^5 ⊗ ℂ^4) + C = randn((ℂ^5)' ⊗ (ℂ^6)' ← ℂ^4 ⊗ (ℂ^3)' ⊗ (ℂ^2)') A′ = force_planar(A) C′ = force_planar(C) p = ((4, 3), (5, 2, 1)) - @test force_planar(tensoradd!(C, p, A, :N, true, true)) ≈ + @test force_planar(tensoradd!(C, A, p, false, true, true)) ≈ planaradd!(C′, A′, p, true, true) end @testset "planartrace" begin - A = TensorMap(randn, ℂ^2 ⊗ ℂ^3 ← ℂ^2 ⊗ ℂ^5 ⊗ ℂ^4) - C = TensorMap(randn, (ℂ^5)' ⊗ ℂ^3 ← ℂ^4) + A = randn(ℂ^2 ⊗ ℂ^3 ← ℂ^2 ⊗ ℂ^5 ⊗ ℂ^4) + C = randn((ℂ^5)' ⊗ ℂ^3 ← ℂ^4) A′ = force_planar(A) C′ = force_planar(C) p = ((4, 2), (5,)) q = ((1,), (3,)) - @test force_planar(tensortrace!(C, p, A, q, :N, true, true)) ≈ + @test force_planar(tensortrace!(C, A, p, q, false, true, true)) ≈ planartrace!(C′, A′, p, q, true, true) end @testset "planarcontract" begin - A = TensorMap(randn, ℂ^2 ⊗ ℂ^3 ← ℂ^2 ⊗ ℂ^5 ⊗ ℂ^4) - B = TensorMap(randn, ℂ^2 ⊗ ℂ^4 ← ℂ^4 ⊗ ℂ^3) - C = TensorMap(randn, (ℂ^5)' ⊗ (ℂ^2)' ⊗ ℂ^2 ← (ℂ^2)' ⊗ ℂ^4) + A = randn(ℂ^2 ⊗ ℂ^3 ← ℂ^2 ⊗ ℂ^5 ⊗ ℂ^4) + B = randn(ℂ^2 ⊗ ℂ^4 ← ℂ^4 ⊗ ℂ^3) + C = randn((ℂ^5)' ⊗ (ℂ^2)' ⊗ ℂ^2 ← (ℂ^2)' ⊗ ℂ^4) A′ = force_planar(A) B′ = force_planar(B) @@ -64,7 +66,7 @@ end pB = ((2, 4), (1, 3)) pAB = ((3, 2, 1), (4, 5)) - @test force_planar(tensorcontract!(C, pAB, A, pA, :N, B, pB, :N, true, true)) ≈ + @test force_planar(tensorcontract!(C, A, pA, false, B, pB, false, pAB, true, true)) ≈ planarcontract!(C′, A′, pA, B′, pB, pAB, true, true) end end @@ -74,18 +76,18 @@ end @testset "contractcheck" begin V = ℂ^2 - A = TensorMap(rand, T, V ⊗ V ← V) - B = TensorMap(rand, T, V ⊗ V ← V') + A = rand(T, V ⊗ V ← V) + B = rand(T, V ⊗ V ← V') @tensor C1[i j; k l] := A[i j; m] * B[k l; m] @tensor contractcheck = true C2[i j; k l] := A[i j; m] * B[k l; m] @test C1 ≈ C2 - B2 = TensorMap(rand, T, V ⊗ V ← V) # wrong duality for third space + B2 = rand(T, V ⊗ V ← V) # wrong duality for third space @test_throws SpaceMismatch("incompatible spaces for m: $V ≠ $(V')") begin @tensor contractcheck = true C3[i j; k l] := A[i j; m] * B2[k l; m] end - A = TensorMap(rand, T, V ← V ⊗ V) - B = TensorMap(rand, T, V ⊗ V ← V) + A = rand(T, V ← V ⊗ V) + B = rand(T, V ⊗ V ← V) @planar C1[i; j] := A[i; k l] * τ[k l; m n] * B[m n; j] @planar contractcheck = true C2[i; j] := A[i; k l] * τ[k l; m n] * B[m n; j] @test C1 ≈ C2 @@ -101,10 +103,10 @@ end # ∂AC # ------- - x = TensorMap(randn, T, Vmps ⊗ P ← Vmps) - O = TensorMap(randn, T, Vmpo ⊗ P ← P ⊗ Vmpo) - GL = TensorMap(randn, T, Vmps ⊗ Vmpo' ← Vmps) - GR = TensorMap(randn, T, Vmps ⊗ Vmpo ← Vmps) + x = randn(T, Vmps ⊗ P ← Vmps) + O = randn(T, Vmpo ⊗ P ← P ⊗ Vmpo) + GL = randn(T, Vmps ⊗ Vmpo' ← Vmps) + GR = randn(T, Vmps ⊗ Vmpo ← Vmps) x′ = force_planar(x) O′ = force_planar(O) @@ -117,7 +119,7 @@ end # ∂AC2 # ------- - x2 = TensorMap(randn, T, Vmps ⊗ P ← Vmps ⊗ P') + x2 = randn(T, Vmps ⊗ P ← Vmps ⊗ P') x2′ = force_planar(x2) @tensor contractcheck = true y2[-1 -2; -3 -4] := GL[-1 7; 6] * x2[6 5; 1 3] * O[7 -2; 5 4] * O[4 -4; 3 2] * @@ -128,7 +130,7 @@ end # transfer matrix # ---------------- - v = TensorMap(randn, T, Vmps ← Vmps) + v = randn(T, Vmps ← Vmps) v′ = force_planar(v) @tensor ρ[-1; -2] := x[-1 2; 1] * conj(x[-2 2; 3]) * v[1; 3] @planar ρ′[-1; -2] := x′[-1 2; 1] * conj(x′[-2 2; 3]) * v′[1; 3] @@ -162,10 +164,10 @@ end @testset "MERA networks" begin Vmera = ℂ^2 - u = TensorMap(randn, T, Vmera ⊗ Vmera ← Vmera ⊗ Vmera) - w = TensorMap(randn, T, Vmera ⊗ Vmera ← Vmera) - ρ = TensorMap(randn, T, Vmera ⊗ Vmera ⊗ Vmera ← Vmera ⊗ Vmera ⊗ Vmera) - h = TensorMap(randn, T, Vmera ⊗ Vmera ⊗ Vmera ← Vmera ⊗ Vmera ⊗ Vmera) + u = randn(T, Vmera ⊗ Vmera ← Vmera ⊗ Vmera) + w = randn(T, Vmera ⊗ Vmera ← Vmera) + ρ = randn(T, Vmera ⊗ Vmera ⊗ Vmera ← Vmera ⊗ Vmera ⊗ Vmera) + h = randn(T, Vmera ⊗ Vmera ⊗ Vmera ← Vmera ⊗ Vmera ⊗ Vmera) u′ = force_planar(u) w′ = force_planar(w) @@ -193,8 +195,8 @@ end T = Float64 V1 = ℂ^2 V2 = ℂ^3 - t1 = TensorMap(rand, T, V1 ← V2) - t2 = TensorMap(rand, T, V2 ← V1) + t1 = rand(T, V1 ← V2) + t2 = rand(T, V2 ← V1) tr1 = @planar opt = true t1[a; b] * t2[b; a] / 2 tr2 = @planar opt = true t1[d; a] * t2[b; c] * 1 / 2 * τ[c b; a d] 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 diff --git a/test/tensors.jl b/test/tensors.jl index 185275ce..d9712d0c 100644 --- a/test/tensors.jl +++ b/test/tensors.jl @@ -83,7 +83,7 @@ for V in spacelist @timedtestset "Basic tensor properties" begin W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5 for T in (Int, Float32, Float64, ComplexF32, ComplexF64, BigFloat) - t = Tensor(zeros, T, W) + t = @constinferred zeros(T, W) @test @constinferred(hash(t)) == hash(deepcopy(t)) @test scalartype(t) == T @test norm(t) == 0 @@ -96,7 +96,7 @@ for V in spacelist @timedtestset "Tensor Dict conversion" begin W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5 for T in (Int, Float32, ComplexF64) - t = TensorMap(rand, T, W) + t = @constinferred rand(T, W) d = convert(Dict, t) @test t == convert(TensorMap, d) end @@ -106,9 +106,12 @@ for V in spacelist W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5 for T in (Int, Float32, ComplexF64) if T == Int - t = TensorMap(sz -> rand(-20:20, sz), W) + t = TensorMap{T}(undef, W) + for (_, b) in blocks(t) + rand!(b, -20:20) + end else - t = TensorMap(randn, T, W) + t = @constinferred randn(T, W) end a = @constinferred convert(Array, t) @test t ≈ @constinferred TensorMap(a, W) @@ -118,7 +121,7 @@ for V in spacelist @timedtestset "Basic linear algebra" begin W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5 for T in (Float32, ComplexF64) - t = TensorMap(rand, T, W) + t = @constinferred rand(T, W) @test scalartype(t) == T @test space(t) == W @test space(t') == W' @@ -136,7 +139,7 @@ for V in spacelist @test norm(t + t, p) ≈ 2 * norm(t, p) @test norm(t) ≈ norm(t') - t2 = TensorMap(rand, T, W) + t2 = @constinferred rand!(similar(t)) β = rand(T) @test @constinferred(dot(β * t2, α * t)) ≈ conj(β) * α * conj(dot(t, t2)) @test dot(t2, t) ≈ conj(dot(t, t2)) @@ -159,8 +162,8 @@ for V in spacelist @timedtestset "Basic linear algebra: test via conversion" begin W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5 for T in (Float32, ComplexF64) - t = TensorMap(rand, T, W) - t2 = TensorMap(rand, T, W) + t = rand(T, W) + t2 = @constinferred rand!(similar(t)) @test norm(t, 2) ≈ norm(convert(Array, t), 2) @test dot(t2, t) ≈ dot(convert(Array, t2), convert(Array, t)) α = rand(T) @@ -171,7 +174,7 @@ for V in spacelist @timedtestset "Real and imaginary parts" begin W = V1 ⊗ V2 for T in (Float64, ComplexF64, ComplexF32) - t = TensorMap(randn, T, W, W) + t = @constinferred randn(T, W, W) @test real(convert(Array, t)) == convert(Array, @constinferred real(t)) @test imag(convert(Array, t)) == convert(Array, @constinferred imag(t)) end @@ -179,7 +182,7 @@ for V in spacelist end @timedtestset "Tensor conversion" begin W = V1 ⊗ V2 - t = TensorMap(randn, Float64, W, W) + t = @constinferred randn(W ← W) @test typeof(convert(TensorMap, t')) == typeof(t) tc = complex(t) @test convert(typeof(tc), t) == tc @@ -190,7 +193,7 @@ for V in spacelist end @timedtestset "diag/diagm" begin W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5 - t = TensorMap(randn, ComplexF64, W) + t = randn(ComplexF64, W) d = LinearAlgebra.diag(t) D = LinearAlgebra.diagm(codomain(t), domain(t), d) @test LinearAlgebra.isdiag(D) @@ -198,8 +201,8 @@ for V in spacelist end @timedtestset "Permutations: test via inner product invariance" begin W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5 - t = Tensor(rand, ComplexF64, W) - t′ = Tensor(rand, ComplexF64, W) + t = rand(ComplexF64, W) + t′ = randn!(similar(t)) for k in 0:5 for p in permutations(1:5) p1 = ntuple(n -> p[n], k) @@ -221,7 +224,7 @@ for V in spacelist if BraidingStyle(I) isa Bosonic && hasfusiontensor(I) @timedtestset "Permutations: test via conversion" begin W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5 - t = Tensor(rand, ComplexF64, W) + t = rand(ComplexF64, W) a = convert(Array, t) for k in 0:5 for p in permutations(1:5) @@ -243,7 +246,7 @@ for V in spacelist end end @timedtestset "Full trace: test self-consistency" begin - t = Tensor(rand, ComplexF64, V1 ⊗ V2' ⊗ V2 ⊗ V1') + t = rand(ComplexF64, V1 ⊗ V2' ⊗ V2 ⊗ V1') t2 = permute(t, ((1, 2), (4, 3))) s = @constinferred tr(t2) @test conj(s) ≈ tr(t2') @@ -261,7 +264,7 @@ for V in spacelist @test ss ≈ s3 end @timedtestset "Partial trace: test self-consistency" begin - t = Tensor(rand, ComplexF64, V1 ⊗ V2' ⊗ V3 ⊗ V2 ⊗ V1' ⊗ V3') + t = rand(ComplexF64, V1 ⊗ V2' ⊗ V3 ⊗ V2 ⊗ V1' ⊗ V3') @tensor t2[a, b] := t[c, d, b, d, c, a] @tensor t4[a, b, c, d] := t[d, e, b, e, c, a] @tensor t5[a, b] := t4[a, b, c, c] @@ -269,15 +272,15 @@ for V in spacelist end if BraidingStyle(I) isa Bosonic && hasfusiontensor(I) @timedtestset "Trace: test via conversion" begin - t = Tensor(rand, ComplexF64, V1 ⊗ V2' ⊗ V3 ⊗ V2 ⊗ V1' ⊗ V3') + t = rand(ComplexF64, V1 ⊗ V2' ⊗ V3 ⊗ V2 ⊗ V1' ⊗ V3') @tensor t2[a, b] := t[c, d, b, d, c, a] @tensor t3[a, b] := convert(Array, t)[c, d, b, d, c, a] @test t3 ≈ convert(Array, t2) end end @timedtestset "Trace and contraction" begin - t1 = Tensor(rand, ComplexF64, V1 ⊗ V2 ⊗ V3) - t2 = Tensor(rand, ComplexF64, V2' ⊗ V4 ⊗ V1') + t1 = rand(ComplexF64, V1 ⊗ V2 ⊗ V3) + t2 = rand(ComplexF64, V2' ⊗ V4 ⊗ V1') t3 = t1 ⊗ t2 @tensor ta[a, b] := t1[x, y, a] * t2[y, b, x] @tensor tb[a, b] := t3[x, y, a, y, b, x] @@ -285,11 +288,11 @@ for V in spacelist end if BraidingStyle(I) isa Bosonic && hasfusiontensor(I) @timedtestset "Tensor contraction: test via conversion" begin - A1 = TensorMap(randn, ComplexF64, V1' * V2', V3') - A2 = TensorMap(randn, ComplexF64, V3 * V4, V5) - rhoL = TensorMap(randn, ComplexF64, V1, V1) - rhoR = TensorMap(randn, ComplexF64, V5, V5)' # test adjoint tensor - H = TensorMap(randn, ComplexF64, V2 * V4, V2 * V4) + A1 = randn(ComplexF64, V1' * V2', V3') + A2 = randn(ComplexF64, V3 * V4, V5) + rhoL = randn(ComplexF64, V1, V1) + rhoR = randn(ComplexF64, V5, V5)' # test adjoint tensor + H = randn(ComplexF64, V2 * V4, V2 * V4) @tensor HrA12[a, s1, s2, c] := rhoL[a, a'] * conj(A1[a', t1, b]) * A2[b, t2, c'] * rhoR[c', c] * H[s1, s2, t1, t2] @@ -307,9 +310,9 @@ for V in spacelist W1 = V1 ⊗ V2 ⊗ V3 W2 = V4 ⊗ V5 for T in (Float64, ComplexF64) - t1 = TensorMap(rand, T, W1, W1) - t2 = TensorMap(rand, T, W2, W2) - t = TensorMap(rand, T, W1, W2) + t1 = rand(T, W1, W1) + t2 = rand(T, W2, W2) + t = rand(T, W1, W2) @test t1 * (t1 \ t) ≈ t @test (t / t2) * t2 ≈ t @test t1 \ one(t1) ≈ inv(t1) @@ -326,9 +329,9 @@ for V in spacelist W1 = V1 ⊗ V2 ⊗ V3 W2 = V4 ⊗ V5 for T in (Float32, Float64, ComplexF32, ComplexF64) - t1 = TensorMap(rand, T, W1, W1) - t2 = TensorMap(rand, T, W2, W2) - t = TensorMap(rand, T, W1, W2) + t1 = rand(T, W1, W1) + t2 = rand(T, W2, W2) + t = rand(T, W1, W2) d1 = dim(W1) d2 = dim(W2) At1 = reshape(convert(Array, t1), d1, d1) @@ -362,7 +365,7 @@ for V in spacelist W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5 for T in (Float32, ComplexF64) # Test both a normal tensor and an adjoint one. - ts = (Tensor(rand, T, W), Tensor(rand, T, W)') + ts = (rand(T, W), rand(T, W)') for t in ts @testset "leftorth with $alg" for alg in (TensorKit.QR(), TensorKit.QRpos(), @@ -427,7 +430,7 @@ for V in spacelist end end @testset "empty tensor" begin - t = TensorMap(randn, T, V1 ⊗ V2, typeof(V1)()) + t = randn(T, V1 ⊗ V2, typeof(V1)()) @testset "leftorth with $alg" for alg in (TensorKit.QR(), TensorKit.QRpos(), TensorKit.QL(), TensorKit.QLpos(), @@ -467,7 +470,7 @@ for V in spacelist end end - t = Tensor(rand, T, V1 ⊗ V1' ⊗ V2 ⊗ V2') + t = rand(T, V1 ⊗ V1' ⊗ V2 ⊗ V2') @testset "eig and isposdef" begin D, V = eigen(t, ((1, 3), (2, 4))) t2 = permute(t, ((1, 3), (2, 4))) @@ -505,8 +508,8 @@ for V in spacelist for T in (Float32, ComplexF64) for p in (1, 2, 3, Inf) # Test both a normal tensor and an adjoint one. - ts = (TensorMap(randn, T, V1 ⊗ V2 ⊗ V3, V4 ⊗ V5), - TensorMap(randn, T, V4 ⊗ V5, V1 ⊗ V2 ⊗ V3)') + ts = (randn(T, V1 ⊗ V2 ⊗ V3, V4 ⊗ V5), + randn(T, V4 ⊗ V5, V1 ⊗ V2 ⊗ V3)') for t in ts U₀, S₀, V₀, = tsvd(t) t = rmul!(t, 1 / norm(S₀, p)) @@ -536,7 +539,7 @@ for V in spacelist @timedtestset "Tensor functions" begin W = V1 ⊗ V2 for T in (Float64, ComplexF64) - t = TensorMap(randn, T, W, W) + t = randn(T, W, W) s = dim(W) expt = @constinferred exp(t) @test reshape(convert(Array, expt), (s, s)) ≈ @@ -578,11 +581,11 @@ for V in spacelist end @timedtestset "Sylvester equation" begin for T in (Float32, ComplexF64) - tA = TensorMap(rand, T, V1 ⊗ V3, V1 ⊗ V3) - tB = TensorMap(rand, T, V2 ⊗ V4, V2 ⊗ V4) + tA = rand(T, V1 ⊗ V3, V1 ⊗ V3) + tB = rand(T, V2 ⊗ V4, V2 ⊗ V4) tA = 3 // 2 * leftorth(tA; alg=Polar())[1] tB = 1 // 5 * leftorth(tB; alg=Polar())[1] - tC = TensorMap(rand, T, V1 ⊗ V3, V2 ⊗ V4) + tC = rand(T, V1 ⊗ V3, V2 ⊗ V4) t = @constinferred sylvester(tA, tB, tC) @test codomain(t) == V1 ⊗ V3 @test domain(t) == V2 ⊗ V4 @@ -596,8 +599,8 @@ for V in spacelist end @timedtestset "Tensor product: test via norm preservation" begin for T in (Float32, ComplexF64) - t1 = TensorMap(rand, T, V2 ⊗ V3 ⊗ V1, V1 ⊗ V2) - t2 = TensorMap(rand, T, V2 ⊗ V1 ⊗ V3, V1 ⊗ V1) + t1 = rand(T, V2 ⊗ V3 ⊗ V1, V1 ⊗ V2) + t2 = rand(T, V2 ⊗ V1 ⊗ V3, V1 ⊗ V1) t = @constinferred (t1 ⊗ t2) @test norm(t) ≈ norm(t1) * norm(t2) end @@ -605,8 +608,8 @@ for V in spacelist if BraidingStyle(I) isa Bosonic && hasfusiontensor(I) @timedtestset "Tensor product: test via conversion" begin for T in (Float32, ComplexF64) - t1 = TensorMap(rand, T, V2 ⊗ V3 ⊗ V1, V1) - t2 = TensorMap(rand, T, V2 ⊗ V1 ⊗ V3, V2) + t1 = rand(T, V2 ⊗ V3 ⊗ V1, V1) + t2 = rand(T, V2 ⊗ V1 ⊗ V3, V2) t = @constinferred (t1 ⊗ t2) d1 = dim(codomain(t1)) d2 = dim(codomain(t2)) @@ -621,8 +624,8 @@ for V in spacelist end @timedtestset "Tensor product: test via tensor contraction" begin for T in (Float32, ComplexF64) - t1 = Tensor(rand, T, V2 ⊗ V3 ⊗ V1) - t2 = Tensor(rand, T, V2 ⊗ V1 ⊗ V3) + t1 = rand(T, V2 ⊗ V3 ⊗ V1) + t2 = rand(T, V2 ⊗ V1 ⊗ V3) t = @constinferred (t1 ⊗ t2) @tensor t′[1, 2, 3, 4, 5, 6] := t1[1, 2, 3] * t2[4, 5, 6] @test t ≈ t′ @@ -636,8 +639,8 @@ end V1, V2, V3, V4, V5 = Vlist1 W1, W2, W3, W4, W5 = Vlist2 for T in (Float32, ComplexF64) - t1 = TensorMap(rand, T, V1 ⊗ V2, V3' ⊗ V4) - t2 = TensorMap(rand, T, W2, W1 ⊗ W1') + t1 = rand(T, V1 ⊗ V2, V3' ⊗ V4) + t2 = rand(T, W2, W1 ⊗ W1') t = @constinferred (t1 ⊠ t2) d1 = dim(codomain(t1)) d2 = dim(codomain(t2))