diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 0688718f..925f3903 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -14,8 +14,145 @@ concurrency: cancel-in-progress: true jobs: + changes: + runs-on: ubuntu-latest + permissions: + pull-requests: read + outputs: + sectors: ${{ steps.filter.outputs.sectors }} + fusiontrees: ${{ steps.filter.outputs.fusiontrees }} + spaces: ${{ steps.filter.outputs.spaces }} + steps: + - uses: actions/checkout@v4 + - uses: dorny/paths-filter@v2 + id: filter + with: + filters: | + sectors: + - 'src/sectors/**' + - 'test/sectors.jl' + fusiontrees: + - 'src/fusiontrees/**' + - 'test/fusiontrees.jl' + spaces: + - 'src/spaces/**' + - 'test/spaces.jl' + test-sectors: + needs: changes + if: ${{ needs.changes.outputs.sectors == 'true' }} + name: "Sectors : Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}" + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.6' # LTS version + - '1' # automatically expands to the latest stable 1.x release of Julia + os: + - ubuntu-latest + - macOS-latest + - windows-latest + arch: + - x64 + - x86 + exclude: + - os: macOS-latest + arch: x86 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v1 + - uses: julia-actions/julia-buildpkg@latest + - uses: lkdvos/julia-runtest@main + with: + suffix: "sectors" + env: + JULIA_NUM_THREADS: 4 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v3 + with: + file: lcov.info + test-spaces: + needs: changes + if: ${{ needs.changes.outputs.spaces == 'true' }} + name: "Spaces : Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}" + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.6' # LTS version + - '1' # automatically expands to the latest stable 1.x release of Julia + os: + - ubuntu-latest + - macOS-latest + - windows-latest + arch: + - x64 + - x86 + exclude: + - os: macOS-latest + arch: x86 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v1 + - uses: julia-actions/julia-buildpkg@latest + - uses: lkdvos/julia-runtest@main + with: + suffix: "spaces" + env: + JULIA_NUM_THREADS: 4 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v3 + with: + file: lcov.info + test-fusiontrees: + needs: changes + if: ${{ needs.changes.outputs.fusiontrees == 'true' }} + name: "Fusiontrees : Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}" + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.6' # LTS version + - '1' # automatically expands to the latest stable 1.x release of Julia + os: + - ubuntu-latest + - macOS-latest + - windows-latest + arch: + - x64 + - x86 + exclude: + - os: macOS-latest + arch: x86 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v1 + - uses: julia-actions/julia-buildpkg@latest + - uses: lkdvos/julia-runtest@main + with: + suffix: "fusiontrees" + env: + JULIA_NUM_THREADS: 4 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v3 + with: + file: lcov.info test: - name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} + name: "Tests : Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}" runs-on: ${{ matrix.os }} strategy: fail-fast: false @@ -41,7 +178,9 @@ jobs: arch: ${{ matrix.arch }} - uses: julia-actions/cache@v1 - uses: julia-actions/julia-buildpkg@latest - - uses: julia-actions/julia-runtest@latest + - uses: lkdvos/julia-runtest@main + with: + suffix: '-sectors -spaces -fusiontrees' env: JULIA_NUM_THREADS: 4 - uses: julia-actions/julia-processcoverage@v1 @@ -71,6 +210,8 @@ jobs: arch: ${{ matrix.arch }} - uses: julia-actions/cache@v1 - uses: julia-actions/julia-buildpkg@latest - - uses: julia-actions/julia-runtest@latest + - uses: lkdvos/julia-runtest@main + with: + suffix: tensors env: JULIA_NUM_THREADS: 4 diff --git a/Project.toml b/Project.toml index 97b3824e..457e3120 100644 --- a/Project.toml +++ b/Project.toml @@ -32,17 +32,14 @@ WignerSymbols = "1,2" julia = "1.6" [extras] +CategoryData = "8fccf25a-f50e-468c-8fba-3cb130506274" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" -HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" 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 = ["Combinatorics", "HalfIntegers", "LinearAlgebra", "Random", "TensorOperations", "Test", "TestExtras", "WignerSymbols", "ChainRulesCore", "ChainRulesTestUtils", "FiniteDifferences"] +test = ["Combinatorics", "Random", "Test", "TestExtras", "ChainRulesCore", "ChainRulesTestUtils", "FiniteDifferences", "CategoryData"] diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 3f40d042..a81aabe3 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -11,7 +11,7 @@ module TensorKit export Sector, AbstractIrrep, Irrep export FusionStyle, UniqueFusion, MultipleFusion, MultiplicityFreeFusion, SimpleFusion, GenericFusion -export BraidingStyle, SymmetricBraiding, Bosonic, Fermionic, Anyonic +export BraidingStyle, SymmetricBraiding, Bosonic, Fermionic, Anyonic, NoBraiding export Trivial, Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep export FermionParity, FermionNumber, FermionSpin export FibonacciAnyon, IsingAnyon diff --git a/src/planar/planaroperations.jl b/src/planar/planaroperations.jl index 9f12de11..ee23bcd7 100644 --- a/src/planar/planaroperations.jl +++ b/src/planar/planaroperations.jl @@ -54,7 +54,7 @@ function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, β::Number, backend::Backend...) where {S,N₁,N₂} if BraidingStyle(sectortype(S)) == Bosonic() - return contract(C, A, pA, B, pB, pAB, α, β, backend...) + return contract!(C, A, pA, B, pB, pAB, α, β, backend...) end codA, domA = codomainind(A), domainind(A) diff --git a/src/sectors/sectors.jl b/src/sectors/sectors.jl index ed127f13..d633bda3 100644 --- a/src/sectors/sectors.jl +++ b/src/sectors/sectors.jl @@ -397,6 +397,8 @@ function pentagon_equation(a::I, b::I, c::I, d::I; kwargs...) where {I<:Sector} end function hexagon_equation(a::I, b::I, c::I; kwargs...) where {I<:Sector} + BraidingStyle(I) isa NoBraiding && + throw(ArgumentError("Hexagon equation only defined for sectors with braiding")) for e in ⊗(c, a), f in ⊗(c, b) for d in intersect(⊗(e, b), ⊗(a, f)) if FusionStyle(I) isa MultiplicityFreeFusion diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index 594e4feb..26546986 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -185,6 +185,33 @@ end 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₂}, + fusiontreetransform, + α::Number, + β::Number, + backend::Backend...) where {S,N₁,N₂} + return add_transform!(tdst, copy(tsrc), (p₁, p₂), fusiontreetransform, α, β, backend...) +end + +# VectorInterface +# --------------- +# TODO + +# Planar operations +# ----------------- +function planaradd!(C::AbstractTensorMap{S,N₁,N₂}, + A::BraidingTensor{S}, + p::Index2Tuple{N₁,N₂}, + α::Number, β::Number, + backend::Backend...) where {S,N₁,N₂} + return planaradd!(C, copy(A), p, α, β, backend...) +end + function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, A::BraidingTensor{S}, (oindA, cindA)::Index2Tuple{2,2}, @@ -278,6 +305,27 @@ function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, return C 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::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...) +end + +function planartrace!(C::AbstractTensorMap{S,N₁,N₂}, + A::BraidingTensor{S}, + p::Index2Tuple{N₁,N₂}, q::Index2Tuple{N₃,N₃}, + α::Number, β::Number, + backend::Backend...) where {S,N₁,N₂,N₃} + return planartrace!(C, copy(A), p, q, α, β, backend...) +end + # function planarcontract!(C::AbstractTensorMap{S,N₁,N₂}, # A::BraidingTensor{S}, # (oindA, cindA)::Index2Tuple{0,4}, @@ -538,5 +586,3 @@ end # end # return C # end - -has_shared_permute(t::BraidingTensor, args...) = false diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl index 9dea8d58..338b701b 100644 --- a/src/tensors/indexmanipulations.jl +++ b/src/tensors/indexmanipulations.jl @@ -28,7 +28,7 @@ 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::TensorMap{S}, (p₁, p₂)::Index2Tuple{N₁,N₂}; +function permute(t::AbstractTensorMap{S}, (p₁, p₂)::Index2Tuple{N₁,N₂}; copy::Bool=false) where {S,N₁,N₂} cod = ProductSpace{S,N₁}(map(n -> space(t, n), p₁)) dom = ProductSpace{S,N₂}(map(n -> dual(space(t, n)), p₂)) @@ -108,7 +108,7 @@ 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::TensorMap{S}, (p₁, p₂)::Index2Tuple, levels::IndexTuple; +function braid(t::AbstractTensorMap{S}, (p₁, p₂)::Index2Tuple, levels::IndexTuple; copy::Bool=false) where {S} @assert length(levels) == numind(t) if BraidingStyle(sectortype(S)) isa SymmetricBraiding @@ -161,7 +161,7 @@ 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::TensorMap{S}, +function LinearAlgebra.transpose(t::AbstractTensorMap{S}, (p₁, p₂)::Index2Tuple=_transpose_indices(t); copy::Bool=false) where {S} if sectortype(S) === Trivial diff --git a/test/ad.jl b/test/ad.jl index 3f87afca..cc446f6e 100644 --- a/test/ad.jl +++ b/test/ad.jl @@ -2,7 +2,8 @@ using ChainRulesCore using ChainRulesTestUtils using Random using FiniteDifferences -using LinearAlgebra +using TensorOperations +using TensorOperations: tensoralloc_add, tensoralloc_contract ## Test utility # ------------- @@ -120,16 +121,16 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), ℂ[Z2Irrep](0 => 3, 1 => 2)', ℂ[Z2Irrep](0 => 2, 1 => 3), ℂ[Z2Irrep](0 => 2, 1 => 2)), - (ℂ[U1Irrep](0 => 1, 1 => 2, -1 => 2), - ℂ[U1Irrep](0 => 3, 1 => 1, -1 => 1), + (ℂ[U1Irrep](0 => 1, 1 => 2, -1 => 1), + ℂ[U1Irrep](0 => 1, 1 => 1, -1 => 1), ℂ[U1Irrep](0 => 2, 1 => 2, -1 => 1)', - ℂ[U1Irrep](0 => 1, 1 => 2, -1 => 2), - ℂ[U1Irrep](0 => 1, 1 => 3, -1 => 2)'), - (ℂ[SU2Irrep](0 => 3, 1 // 2 => 1), + ℂ[U1Irrep](0 => 1, 1 => 2, -1 => 1), + ℂ[U1Irrep](0 => 1, 1 => 2, -1 => 2)'), + (ℂ[SU2Irrep](0 => 1, 1 // 2 => 1), ℂ[SU2Irrep](0 => 2, 1 => 1), ℂ[SU2Irrep](1 // 2 => 1, 1 => 1)', ℂ[SU2Irrep](0 => 2, 1 // 2 => 2), - ℂ[SU2Irrep](0 => 1, 1 // 2 => 1, 3 // 2 => 1)')) + ℂ[SU2Irrep](0 => 1, 1 // 2 => 2)')) @testset "Automatic Differentiation ($(eltype(V)))" verbose = true for V in Vlist @testset "Basic Linear Algebra ($T)" for T in (Float64, ComplexF64) @@ -171,22 +172,22 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), α = randn(T) β = randn(T) - C = _randomize!(TensorOperations.tensoralloc_add(T, pC, A, :N, false)) + C = _randomize!(tensoralloc_add(T, pC, A, :N, false)) test_rrule(tensortrace!, C, pC, A, pA, :N, α, β; atol, rtol) - C = _randomize!(TensorOperations.tensoralloc_add(T, pC, A, :C, false)) + C = _randomize!(tensoralloc_add(T, pC, A, :C, false)) test_rrule(tensortrace!, C, pC, A, pA, :C, α, β; atol, rtol) end @testset "tensoradd!" begin p = ((1, 3, 2), (5, 4)) A = TensorMap(randn, T, V[1] ⊗ V[2] ← V[3] ⊗ V[4] ⊗ V[5]) - C = _randomize!(TensorOperations.tensoralloc_add(T, p, A, :N, false)) + C = _randomize!(tensoralloc_add(T, p, A, :N, false)) α = randn(T) β = randn(T) test_rrule(tensoradd!, C, p, A, :N, α, β; atol, rtol) - C = _randomize!(TensorOperations.tensoralloc_add(T, p, A, :C, false)) + C = _randomize!(tensoralloc_add(T, p, A, :C, false)) test_rrule(tensoradd!, C, p, A, :C, α, β; atol, rtol) end @@ -199,22 +200,22 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), α = randn(T) β = randn(T) - C = _randomize!(TensorOperations.tensoralloc_contract(T, pC, A, pA, :N, - B, pB, :N, false)) + C = _randomize!(tensoralloc_contract(T, pC, A, pA, :N, + B, pB, :N, false)) test_rrule(tensorcontract!, C, pC, A, pA, :N, B, pB, :N, α, β; atol, rtol) A2 = TensorMap(randn, T, V[1]' ⊗ V[2]' ← V[3]' ⊗ V[4]' ⊗ V[5]') - C = _randomize!(TensorOperations.tensoralloc_contract(T, pC, A2, pA, :C, - B, pB, :N, false)) + C = _randomize!(tensoralloc_contract(T, pC, A2, pA, :C, + B, pB, :N, false)) test_rrule(tensorcontract!, C, pC, A2, pA, :C, B, pB, :N, α, β; atol, rtol) B2 = TensorMap(randn, T, V[3]' ⊗ V[1] ← V[2]') - C = _randomize!(TensorOperations.tensoralloc_contract(T, pC, A, pA, :N, - B2, pB, :C, false)) + C = _randomize!(tensoralloc_contract(T, pC, A, pA, :N, + B2, pB, :C, false)) test_rrule(tensorcontract!, C, pC, A, pA, :N, B2, pB, :C, α, β; atol, rtol) - C = _randomize!(TensorOperations.tensoralloc_contract(T, pC, A2, pA, :C, - B2, pB, :C, false)) + C = _randomize!(tensoralloc_contract(T, pC, A2, pA, :C, + B2, pB, :C, false)) test_rrule(tensorcontract!, C, pC, A2, pA, :C, B2, pB, :C, α, β; atol, rtol) end diff --git a/test/braidingtensor.jl b/test/braidingtensor.jl deleted file mode 100644 index b61abc61..00000000 --- a/test/braidingtensor.jl +++ /dev/null @@ -1,134 +0,0 @@ -# TODO: Make into proper tests and integrate in testset - -import TensorKit: BraidingTensor - -V1 = GradedSpace{FermionSpin}(0 => 2, 1 / 2 => 2, 1 => 1, 3 / 2 => 1) - -V2 = GradedSpace{FibonacciAnyon}(:I => 2, :τ => 2) - -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') - - ττ = copy(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')) - @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)) - @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')) - @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)) - @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)) - @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)) - @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')) - @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)) - @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)) - @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)) - @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')) - @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')) - @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')) - @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)) - @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)) - @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')) - @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) -end diff --git a/test/choosetests.jl b/test/choosetests.jl new file mode 100644 index 00000000..f86b15ec --- /dev/null +++ b/test/choosetests.jl @@ -0,0 +1,175 @@ +# This file is largely inspired by the Base Julia file test/choostests.jl + +const TESTNAMES = ["sectors", "spaces", "fusiontrees", "tensors", "planar", "ad"] +const TESTGROUPS = String[] + +const SECTORNAMES = ["Trivial", "Z2Irrep", "Z3Irrep", "Z4Irrep", "U1Irrep", "CU1Irrep", + "SU2Irrep", "NewSU2Irrep", "FibonacciAnyon", "IsingAnyon", + "FermionParity", "FermionNumber", "FermionSpin", "Z3Irrep ⊠ Z4Irrep", + "FermionNumber ⊠ SU2Irrep", "FermionSpin ⊠ SU2Irrep", + "NewSU2Irrep ⊠ NewSU2Irrep", "FibonacciAnyon", "Object{E6}", + "Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon"] +const DEFAULT_SECTORNAMES = try + if ENV["CI"] == "true" + println("Detected CI environment") + if Sys.iswindows() + ["Trivial", "Z2Irrep", "FermionParity", "Z3Irrep", "U1Irrep", + "FermionNumber", "CU1Irrep", "SU2Irrep"] + elseif Sys.isapple() + ["Trivial", "Z2Irrep", "FermionParity", "Z3Irrep", "FermionNumber", + "FermionSpin"] + else + ["Trivial", "Z2Irrep", "FermionParity", "U1Irrep", "CU1Irrep", "SU2Irrep", + "FermionSpin", "Object{E6}"] + end + else + SECTORNAMES + end +catch + SECTORNAMES +end +const SECTORGROUPS = String[] + +""" + (; tests, exit_on_error, seed) = choosetests(choices=[]) + +Selects a set of tests to be run. `choices` should be a vector of test names and/or +sectornames; if empty or set to `["all"]`, all tests are selected. + +Several options can be passed to `choosetests` by including a special token in the `choices` +argument: +- "--skip", which makes all tests coming after be skipped +- "--exit-on-error" which sets the value of `exit_on_error` +- "--revise" which loads Revise +- "--seed=SEED", which sets the value of `seed` to `SEED` (parsed as an `UInt128`); + `seed` is otherwise initialized randomly. This option can be used to reproduce failed tests. +- "--help", which prints a help message and then skips all tests. +- "--help-list", which prints the options computed without running them. + +The function returns a named tuple with the following elements: +- `tests` is a vector of fully-expanded test names +- `sectors` is a vector of fully-expanded sector names +- `exit_on_error` is true if an error in one test should cancel remaining tests to be run + (otherwise, all tests are run unconditionally) +- `seed` is a seed which will be used to initialize the global RNG for each test to be run. +""" +function choosetests(choices=[]) + tests = [] + sectors = [] + skip_tests = Set() + skip_sectors = Set() + exit_on_error = false + use_revise = false + seed = rand(RandomDevice(), UInt128) + ci_option_passed = false + dryrun = false + + # parse options + for (i, t) in enumerate(choices) + if t == "--skip" + union!(skip_tests, choices[(i + 1):end]) + break + elseif t == "--exit-on-error" + exit_on_error = true + elseif t == "--revise" + use_revise = true + elseif startswith(t, "--seed=") + seed = parse(UInt128, t[(length("--seed=") + 1):end]) + elseif t == "--ci" + ci_option_passed = true + elseif t == "--help-list" + dryrun = true + elseif t == "--help" + println(""" + USAGE: ./julia runtests.jl [options] [tests] + OPTIONS: + --exit-on-error : stop tests immediately when a test group fails + --help : prints this help message + --help-list : prints the options computed without running them + --revise : load Revise + --seed= : set the initial seed for all testgroups (parsed as a UInt128) + --skip ... : skip test or collection tagged with + TESTS: + Can be testsets, such as ($TESTNAMES), + or sector names, such as ($SECTORNAMES). + + Prefixing a name with `-` (such as `-sectors`) can be used to skip a particular test. + """) + return (; tests=[], + sectors=[], + exit_on_error=false, + use_revise=false, + seed=UInt128(0)) + elseif startswith(t, "--") + error("unknown option: $t") + elseif startswith(t, "-") + if t[2:end] in TESTNAMES || t in TESTGROUPS + push!(skip_tests, t[2:end]) + elseif t[2:end] in SECTORNAMES || t in SECTORGROUPS + push!(skip_sectors, t[2:end]) + else + error("unknown test or sector name at $i: $t") + end + else + if t in TESTNAMES || t in TESTGROUPS + push!(tests, t) + elseif t in SECTORNAMES || t in SECTORGROUPS + push!(sectors, t) + elseif t == "all-tests" + append!(tests, TESTNAMES) + elseif t == "all-sectors" + append!(sectors, SECTORNAMES) + elseif t == "default-sectors" + append!(sectors, DEFAULT_SECTORNAMES) + else + error("unknown test or sector name at $i: $t") + end + end + end + + # Filter tests + unhandled_tests = copy(skip_tests) + + if isempty(tests) + append!(tests, TESTNAMES) + end + + # Functionality to add testgroups for inclusion/exclusion: + # function filtertests!(tests, name_or_group, names=[name_or_group]) + # flt = x -> (x != name_or_group && !(x in names)) + # if name_or_group in skip_tests + # filter!(flt, tests) + # pop!(unhandled_tests, name_or_group) + # elseif name_or_group in tests + # filter!(flt, tests) + # prepend!(tests, names) + # end + # end + + filter!(!in(tests), unhandled_tests) + filter!(!in(skip_tests), tests) + + # Filter sectors + unhandled_sectors = copy(skip_sectors) + + if isempty(sectors) + append!(sectors, DEFAULT_SECTORNAMES) + end + + # Functionality to add sectorgroups for inclusion/exclusion: + # function filtersectors!(sectors, name_or_group, names=[name_or_group]) + # flt = x -> (x != name_or_group && !(x in names)) + # if name_or_group in skip_sectors + # filter!(flt, sectors) + # pop!(unhandled_sectors, name_or_group) + # elseif name_or_group in sectors + # filter!(flt, sectors) + # prepend!(sectors, names) + # end + # end + + filter!(!in(sectors), unhandled_sectors) + filter!(!in(skip_sectors), sectors) + + return (; tests, sectors, exit_on_error, use_revise, seed) +end diff --git a/test/fusiontrees.jl b/test/fusiontrees.jl index 61dc7d38..6230c96f 100644 --- a/test/fusiontrees.jl +++ b/test/fusiontrees.jl @@ -1,10 +1,13 @@ println("------------------------------------") println("Fusion Trees") println("------------------------------------") -ti = time() -@timedtestset "Fusion trees for $(TensorKit.type_repr(I))" verbose = true for I in - sectorlist + +using TensorOperations +using TensorKit: ProductSector + +@testset "$(TensorKit.type_repr(I))" verbose = true for I in sectorlist Istr = TensorKit.type_repr(I) + println("Starting tests for $Istr...") N = 5 out = ntuple(n -> randsector(I), N) isdual = ntuple(n -> rand(Bool), N) @@ -18,10 +21,11 @@ ti = time() it = @constinferred fusiontrees(out, in, isdual) @constinferred Nothing iterate(it) f = @constinferred first(it) - @testset "Fusion tree $Istr: printing" begin + + @testset "printing" begin @test eval(Meta.parse(sprint(show, f))) == f end - @testset "Fusion tree $Istr: insertat" begin + @testset "insertat" begin N = 4 out2 = ntuple(n -> randsector(I), N) in2 = rand(collect(⊗(out2...))) @@ -35,19 +39,20 @@ ti = time() isdual1 = Base.setindex(isdual1, false, i) f1 = rand(collect(fusiontrees(out1, in1, isdual1))) - trees = @constinferred TK.insertat(f1, i, f2) + trees = @constinferred TensorKit.insertat(f1, i, f2) @test norm(values(trees)) ≈ 1 - f1a, f1b = @constinferred TK.split(f1, $i) - @test length(TK.insertat(f1b, 1, f1a)) == 1 - @test first(TK.insertat(f1b, 1, f1a)) == (f1 => 1) + f1a, f1b = @constinferred TensorKit.split(f1, $i) + @test length(TensorKit.insertat(f1b, 1, f1a)) == 1 + @test first(TensorKit.insertat(f1b, 1, f1a)) == (f1 => 1) levels = ntuple(identity, N) function _reinsert_partial_tree(t, f) - (t′, c′) = first(TK.insertat(t, 1, f)) + (t′, c′) = first(TensorKit.insertat(t, 1, f)) @test c′ == one(c′) return t′ end + BraidingStyle(I) isa NoBraiding && continue braid_i_to_1 = braid(f1, levels, (i, (1:(i - 1))..., ((i + 1):N)...)) trees2 = Dict(_reinsert_partial_tree(t, f2) => c for (t, c) in braid_i_to_1) trees3 = empty(trees2) @@ -77,7 +82,7 @@ ti = time() end end end - @testset "Fusion tree $Istr: planar trace" begin + @testset "planar trace" begin if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) s = randsector(I) N = 6 @@ -89,7 +94,7 @@ ti = time() T = eltype(af) for i in 1:N - d = @constinferred TK.elementary_trace(f, i) + d = @constinferred TensorKit.elementary_trace(f, i) j = mod1(i + 1, N) inds = collect(1:(N + 1)) inds[i] = inds[j] @@ -101,7 +106,7 @@ ti = time() @test bf ≈ bf′ atol = 1e-12 end - d2 = @constinferred TK.planar_trace(f, (1, 3), (2, 4)) + d2 = @constinferred TensorKit.planar_trace(f, (1, 3), (2, 4)) oind2 = (5, 6, 7) bf2 = tensortrace(af, (:a, :a, :b, :b, :c, :d, :e)) bf2′ = zero(bf2) @@ -110,7 +115,7 @@ ti = time() end @test bf2 ≈ bf2′ atol = 1e-12 - d2 = @constinferred TK.planar_trace(f, (5, 6), (2, 1)) + d2 = @constinferred TensorKit.planar_trace(f, (5, 6), (2, 1)) oind2 = (3, 4, 7) bf2 = tensortrace(af, (:a, :b, :c, :d, :b, :a, :e)) bf2′ = zero(bf2) @@ -119,7 +124,7 @@ ti = time() end @test bf2 ≈ bf2′ atol = 1e-12 - d2 = @constinferred TK.planar_trace(f, (1, 4), (6, 3)) + d2 = @constinferred TensorKit.planar_trace(f, (1, 4), (6, 3)) bf2 = tensortrace(af, (:a, :b, :c, :c, :d, :a, :e)) bf2′ = zero(bf2) for (f2′, coeff) in d2 @@ -129,7 +134,7 @@ ti = time() q1 = (1, 3, 5) q2 = (2, 4, 6) - d3 = @constinferred TK.planar_trace(f, q1, q2) + d3 = @constinferred TensorKit.planar_trace(f, q1, q2) bf3 = tensortrace(af, (:a, :a, :b, :b, :c, :c, :d)) bf3′ = zero(bf3) for (f3′, coeff) in d3 @@ -139,7 +144,7 @@ ti = time() q1 = (1, 3, 5) q2 = (6, 2, 4) - d3 = @constinferred TK.planar_trace(f, q1, q2) + d3 = @constinferred TensorKit.planar_trace(f, q1, q2) bf3 = tensortrace(af, (:a, :b, :b, :c, :c, :a, :d)) bf3′ = zero(bf3) for (f3′, coeff) in d3 @@ -149,7 +154,7 @@ ti = time() q1 = (1, 2, 3) q2 = (6, 5, 4) - d3 = @constinferred TK.planar_trace(f, q1, q2) + d3 = @constinferred TensorKit.planar_trace(f, q1, q2) bf3 = tensortrace(af, (:a, :b, :c, :c, :b, :a, :d)) bf3′ = zero(bf3) for (f3′, coeff) in d3 @@ -159,7 +164,7 @@ ti = time() q1 = (1, 2, 4) q2 = (6, 3, 5) - d3 = @constinferred TK.planar_trace(f, q1, q2) + d3 = @constinferred TensorKit.planar_trace(f, q1, q2) bf3 = tensortrace(af, (:a, :b, :b, :c, :c, :a, :d)) bf3′ = zero(bf3) for (f3′, coeff) in d3 @@ -170,96 +175,96 @@ ti = time() end end end - @testset "Fusion tree $Istr: elementy artin braid" begin - N = length(out) - isdual = ntuple(n -> rand(Bool), N) - for in in ⊗(out...) - for i in 1:(N - 1) - for f in fusiontrees(out, in, isdual) - d1 = @constinferred TK.artin_braid(f, i) - @test norm(values(d1)) ≈ 1 - d2 = empty(d1) - for (f1, coeff1) in d1 - for (f2, coeff2) in TK.artin_braid(f1, i; inv=true) - d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 - end - end - for (f2, coeff2) in d2 - if f2 == f - @test coeff2 ≈ 1 - else - @test isapprox(coeff2, 0; atol=1e-12, rtol=1e-12) - end - end - end - end - end + BraidingStyle(I) isa NoBraiding || @testset "elementy artin braid" begin + N = length(out) + isdual = ntuple(n -> rand(Bool), N) + for in in ⊗(out...) + for i in 1:(N - 1) + for f in fusiontrees(out, in, isdual) + d1 = @constinferred TensorKit.artin_braid(f, i) + @test norm(values(d1)) ≈ 1 + d2 = empty(d1) + for (f1, coeff1) in d1 + for (f2, coeff2) in TensorKit.artin_braid(f1, i; inv=true) + d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 + end + end + for (f2, coeff2) in d2 + if f2 == f + @test coeff2 ≈ 1 + else + @test isapprox(coeff2, 0; atol=1e-12, rtol=1e-12) + end + end + end + end + end - f = rand(collect(it)) - d1 = TK.artin_braid(f, 2) - d2 = empty(d1) - for (f1, coeff1) in d1 - for (f2, coeff2) in TK.artin_braid(f1, 3) - d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 - end - end - d1 = d2 - d2 = empty(d1) - for (f1, coeff1) in d1 - for (f2, coeff2) in TK.artin_braid(f1, 3; inv=true) - d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 - end - end - d1 = d2 - d2 = empty(d1) - for (f1, coeff1) in d1 - for (f2, coeff2) in TK.artin_braid(f1, 2; inv=true) - d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 - end - end - d1 = d2 - for (f1, coeff1) in d1 - if f1 == f - @test coeff1 ≈ 1 - else - @test isapprox(coeff1, 0; atol=1e-12, rtol=1e-12) - end - end - end - @testset "Fusion tree $Istr: braiding and permuting" begin - f = rand(collect(fusiontrees(out, in, isdual))) - p = tuple(randperm(N)...) - ip = invperm(p) + f = rand(collect(it)) + d1 = TensorKit.artin_braid(f, 2) + d2 = empty(d1) + for (f1, coeff1) in d1 + for (f2, coeff2) in TensorKit.artin_braid(f1, 3) + d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 + end + end + d1 = d2 + d2 = empty(d1) + for (f1, coeff1) in d1 + for (f2, coeff2) in TensorKit.artin_braid(f1, 3; inv=true) + d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 + end + end + d1 = d2 + d2 = empty(d1) + for (f1, coeff1) in d1 + for (f2, coeff2) in TensorKit.artin_braid(f1, 2; inv=true) + d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 + end + end + d1 = d2 + for (f1, coeff1) in d1 + if f1 == f + @test coeff1 ≈ 1 + else + @test isapprox(coeff1, 0; atol=1e-12, rtol=1e-12) + end + end + end + BraidingStyle(I) isa NoBraiding || @testset "braiding and permuting" begin + f = rand(collect(fusiontrees(out, in, isdual))) + p = tuple(randperm(N)...) + ip = invperm(p) - levels = ntuple(identity, N) - d = @constinferred braid(f, levels, p) - d2 = Dict{typeof(f),valtype(d)}() - levels2 = p - for (f2, coeff) in d - for (f1, coeff2) in braid(f2, levels2, ip) - d2[f1] = get(d2, f1, zero(coeff)) + coeff2 * coeff - end - end - for (f1, coeff2) in d2 - if f1 == f - @test coeff2 ≈ 1 - else - @test isapprox(coeff2, 0; atol=1e-12, rtol=1e-12) - end - end + levels = ntuple(identity, N) + d = @constinferred braid(f, levels, p) + d2 = Dict{typeof(f),valtype(d)}() + levels2 = p + for (f2, coeff) in d + for (f1, coeff2) in braid(f2, levels2, ip) + d2[f1] = get(d2, f1, zero(coeff)) + coeff2 * coeff + end + end + for (f1, coeff2) in d2 + if f1 == f + @test coeff2 ≈ 1 + else + @test isapprox(coeff2, 0; atol=1e-12, rtol=1e-12) + end + end - if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) - Af = convert(Array, f) - Afp = permutedims(Af, (p..., N + 1)) - Afp2 = zero(Afp) - for (f1, coeff) in d - Afp2 .+= coeff .* convert(Array, f1) - end - @test Afp ≈ Afp2 - end - end + if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) + Af = convert(Array, f) + Afp = permutedims(Af, (p..., N + 1)) + Afp2 = zero(Afp) + for (f1, coeff) in d + Afp2 .+= coeff .* convert(Array, f1) + end + @test Afp ≈ Afp2 + end + end - @testset "Fusion tree $Istr: merging" begin + @testset "merging" begin N = 3 out1 = ntuple(n -> randsector(I), N) in1 = rand(collect(⊗(out1...))) @@ -268,61 +273,64 @@ ti = time() in2 = rand(collect(⊗(out2...))) f2 = rand(collect(fusiontrees(out2, in2))) - @constinferred TK.merge(f1, f2, first(in1 ⊗ in2), 1) + @constinferred TensorKit.merge(f1, f2, first(in1 ⊗ in2), 1) if !(FusionStyle(I) isa GenericFusion) - @constinferred TK.merge(f1, f2, first(in1 ⊗ in2), nothing) - @constinferred TK.merge(f1, f2, first(in1 ⊗ in2)) + @constinferred TensorKit.merge(f1, f2, first(in1 ⊗ in2), nothing) + @constinferred TensorKit.merge(f1, f2, first(in1 ⊗ in2)) end @test dim(in1) * dim(in2) ≈ sum(abs2(coeff) * dim(c) for c in in1 ⊗ in2 for μ in 1:Nsymbol(in1, in2, c) - for (f, coeff) in TK.merge(f1, f2, c, μ)) - - for c in in1 ⊗ in2 - R = Rsymbol(in1, in2, c) - for μ in 1:Nsymbol(in1, in2, c) - μ′ = FusionStyle(I) isa GenericFusion ? μ : nothing - trees1 = TK.merge(f1, f2, c, μ′) + for (f, coeff) in TensorKit.merge(f1, f2, c, μ)) + BraidingStyle(I) isa NoBraiding || for c in in1 ⊗ in2 + R = Rsymbol(in1, in2, c) + for μ in 1:Nsymbol(in1, in2, c) + μ′ = FusionStyle(I) isa GenericFusion ? μ : nothing + trees1 = TensorKit.merge(f1, f2, c, μ′) # test merge and braid interplay - trees2 = Dict{keytype(trees1),complex(valtype(trees1))}() - trees3 = Dict{keytype(trees1),complex(valtype(trees1))}() - for ν in 1:Nsymbol(in2, in1, c) - ν′ = FusionStyle(I) isa GenericFusion ? ν : nothing - for (t, coeff) in TK.merge(f2, f1, c, ν′) - trees2[t] = get(trees2, t, zero(valtype(trees2))) + coeff * R[μ, ν] - end - end - perm = ((N .+ (1:N))..., (1:N)...) - levels = ntuple(identity, 2 * N) - for (t, coeff) in trees1 - for (t′, coeff′) in braid(t, levels, perm) - trees3[t′] = get(trees3, t′, zero(valtype(trees3))) + coeff * coeff′ - end - end - for (t, coeff) in trees3 - coeff′ = get(trees2, t, zero(coeff)) - @test isapprox(coeff, coeff′; atol=1e-12, rtol=1e-12) - end + trees2 = Dict{keytype(trees1),complex(valtype(trees1))}() + trees3 = Dict{keytype(trees1),complex(valtype(trees1))}() + for ν in 1:Nsymbol(in2, in1, c) + ν′ = FusionStyle(I) isa GenericFusion ? ν : nothing + for (t, coeff) in TensorKit.merge(f2, f1, c, ν′) + trees2[t] = get(trees2, t, zero(valtype(trees2))) + + coeff * R[μ, ν] + end + end + perm = ((N .+ (1:N))..., (1:N)...) + levels = ntuple(identity, 2 * N) + for (t, coeff) in trees1 + for (t′, coeff′) in braid(t, levels, perm) + trees3[t′] = get(trees3, t′, zero(valtype(trees3))) + + coeff * coeff′ + end + end + for (t, coeff) in trees3 + coeff′ = get(trees2, t, zero(coeff)) + @test isapprox(coeff, coeff′; atol=1e-12, rtol=1e-12) + end # test via conversion - if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) - Af1 = convert(Array, f1) - Af2 = convert(Array, f2) - Af0 = convert(Array, - FusionTree((f1.coupled, f2.coupled), c, (false, false), - (), (μ,))) - _Af = TensorOperations.tensorcontract(1:(N + 2), Af1, [1:N; -1], - Af0, [-1; N + 1; N + 2]) - Af = TensorOperations.tensorcontract(1:(2N + 1), Af2, [N .+ (1:N); -1], - _Af, [1:N; -1; 2N + 1]) - Af′ = zero(Af) - for (f, coeff) in trees1 - Af′ .+= coeff .* convert(Array, f) - end - @test Af ≈ Af′ - end - end - end + if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I) + Af1 = convert(Array, f1) + Af2 = convert(Array, f2) + Af0 = convert(Array, + FusionTree((f1.coupled, f2.coupled), c, + (false, false), + (), (μ,))) + _Af = TensorOperations.tensorcontract(1:(N + 2), Af1, [1:N; -1], + Af0, [-1; N + 1; N + 2]) + Af = TensorOperations.tensorcontract(1:(2N + 1), Af2, + [N .+ (1:N); -1], + _Af, [1:N; -1; 2N + 1]) + Af′ = zero(Af) + for (f, coeff) in trees1 + Af′ .+= coeff .* convert(Array, f) + end + @test Af ≈ Af′ + end + end + end end if I <: ProductSector @@ -340,14 +348,14 @@ ti = time() f1 = rand(collect(fusiontrees(out, incoming, ntuple(n -> rand(Bool), N)))) f2 = rand(collect(fusiontrees(out[randperm(N)], incoming, ntuple(n -> rand(Bool), N)))) - @testset "Double fusion tree $Istr: repartioning" begin + @testset "repartioning" begin for n in 0:(2 * N) - d = @constinferred TK.repartition(f1, f2, $n) + d = @constinferred TensorKit.repartition(f1, f2, $n) @test dim(incoming) ≈ sum(abs2(coef) * dim(f1.coupled) for ((f1, f2), coef) in d) d2 = Dict{typeof((f1, f2)),valtype(d)}() for ((f1′, f2′), coeff) in d - for ((f1′′, f2′′), coeff2) in TK.repartition(f1′, f2′, N) + for ((f1′′, f2′′), coeff2) in TensorKit.repartition(f1′, f2′, N) d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + coeff2 * coeff end end @@ -385,7 +393,7 @@ ti = time() end end end - @testset "Double fusion tree $Istr: permutation" begin + @testset "permutation" begin if BraidingStyle(I) isa SymmetricBraiding for n in 0:(2N) p = (randperm(2 * N)...,) @@ -441,7 +449,7 @@ ti = time() end end end - @testset "Double fusion tree $Istr: transposition" begin + @testset "transposition" begin for n in 0:(2N) i0 = rand(1:(2N)) p = mod1.(i0 .+ (1:(2N)), 2N) @@ -506,15 +514,15 @@ ti = time() end end end - @testset "Double fusion tree $Istr: planar trace" begin + @testset "planar trace" begin d1 = transpose(f1, f1, (N + 1, 1:N..., ((2N):-1:(N + 3))...), (N + 2,)) - f1front, = TK.split(f1, N - 1) + f1front, = TensorKit.split(f1, N - 1) T = typeof(Fsymbol(one(I), one(I), one(I), one(I), one(I), one(I))[1, 1, 1, 1]) d2 = Dict{typeof((f1front, f1front)),T}() for ((f1′, f2′), coeff′) in d1 for ((f1′′, f2′′), coeff′′) in - TK.planar_trace(f1′, f2′, (2:N...,), (1, ((2N):-1:(N + 3))...), (N + 1,), - (N + 2,)) + TensorKit.planar_trace(f1′, f2′, (2:N...,), (1, ((2N):-1:(N + 3))...), + (N + 1,), (N + 2,)) coeff = coeff′ * coeff′′ d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + coeff end @@ -527,9 +535,6 @@ ti = time() end end end + + println("Finished tests for $Istr.") end -tf = time() -printstyled("Finished fusion tree tests in ", - string(round(tf - ti; sigdigits=3)), - " seconds."; bold=true, color=Base.info_color()) -println() diff --git a/test/planar.jl b/test/planar.jl index 02f58166..af869596 100644 --- a/test/planar.jl +++ b/test/planar.jl @@ -1,72 +1,244 @@ -using TensorKit, TensorOperations, Test -using TensorKit: planaradd!, planartrace!, planarcontract! -using TensorKit: PlanarTrivial, ℙ - -""" - force_planar(obj) - -Replace an object with a planar equivalent -- i.e. one that disallows braiding. -""" -force_planar(V::ComplexSpace) = isdual(V) ? (ℙ^dim(V))' : ℙ^dim(V) -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))) - 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))) - for (c, b) in blocks(tsrc) - copyto!(blocks(tdst)[c ⊠ PlanarTrivial()], b) - end - return tdst -end +println("------------------------------------") +println("Planar") +println("------------------------------------") -@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′ = force_planar(A) - C′ = force_planar(C) - p = ((4, 3), (5, 2, 1)) +using TensorKit: planaradd!, planartrace!, planarcontract!, BraidingTensor, + SymmetricBraiding +using TensorOperations - @test force_planar(tensoradd!(C, p, A, :N, true, true)) ≈ - planaradd!(C′, A′, p, true, true) +@testset "$(TensorKit.type_repr(I))" verbose = true for I in sectorlist + V = smallspace(I) + if isnothing(V) + "No spaces defined for $(TensorKit.type_repr(I)), skipping tests" + continue end + Istr = TensorKit.type_repr(I) + println("Starting tests for $Istr...") + V1, V2, V3, V4, V5 = V + + if BraidingStyle(I) isa SymmetricBraiding + @testset "planaradd" begin + A = TensorMap(randn, V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5) + C = TensorMap(randn, V4' ⊗ V3' ← V5 ⊗ V2' ⊗ V1') + + A′ = force_planar(A) + C′ = force_planar(C) + + p = ((4, 3), (5, 2, 1)) + + @test force_planar(tensoradd!(C, p, A, :N, true, true)) ≈ + planaradd!(C′, A′, p, true, true) + end + @testset "planartrace" begin + A = TensorMap(randn, V1 ⊗ V2 ← V1 ⊗ V4 ⊗ V5) + C = TensorMap(randn, V4' ⊗ V2 ← V5) - @testset "planartrace" begin - A = TensorMap(randn, ℂ^2 ⊗ ℂ^3 ← ℂ^2 ⊗ ℂ^5 ⊗ ℂ^4) - C = TensorMap(randn, (ℂ^5)' ⊗ ℂ^3 ← ℂ^4) - A′ = force_planar(A) - C′ = force_planar(C) - p = ((4, 2), (5,)) - q = ((1,), (3,)) + A′ = force_planar(A) + C′ = force_planar(C) - @test force_planar(tensortrace!(C, p, A, q, :N, true, true)) ≈ - planartrace!(C′, A′, p, q, true, true) + p = ((4, 2), (5,)) + q = ((1,), (3,)) + + @test force_planar(tensortrace!(C, p, A, q, :N, true, true)) ≈ + planartrace!(C′, A′, p, q, true, true) + end + + @testset "planarcontract" begin + A = TensorMap(randn, V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5) + B = TensorMap(randn, V2 ⊗ V5 ← V1 ⊗ V2) + C = TensorMap(randn, V4' ⊗ V3' ⊗ V1 ← V2' ⊗ V1) + + A′ = force_planar(A) + B′ = force_planar(B) + C′ = force_planar(C) + + pA = ((1, 3, 4), (5, 2)) + 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)) ≈ + planarcontract!(C′, A′, pA, B′, pB, pAB, true, true) + end 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) + @testset "BraidingTensor conversion" begin + for (V1, V2) in [(V1, V1), (V1', V1), (V1, V1'), (V1', V1')] + τ = BraidingTensor(V1, V2) - A′ = force_planar(A) - B′ = force_planar(B) - C′ = force_planar(C) + @test domain(τ) == V1 ⊗ V2 + @test codomain(τ) == V2 ⊗ V1 - pA = ((1, 3, 4), (5, 2)) - pB = ((2, 4), (1, 3)) - pAB = ((3, 2, 1), (4, 5)) + for (c, b) in blocks(copy(τ)) + @test b ≈ block(τ, c) + end - @test force_planar(tensorcontract!(C, pAB, A, pA, :N, B, pB, :N, true, true)) ≈ - planarcontract!(C′, A′, pA, B′, pB, pAB, true, true) + @test domain(τ') == codomain(τ) + @test codomain(τ') == domain(τ) + + for (c, b) in blocks(copy(τ')) + @test b ≈ block(τ', c) + end + end end + + t = TensorMap(randn, V1 * V1' * V1' * V1, V1 * V1') + + ττ = copy(BraidingTensor(V1, V1')) + @planar t1[-1 -2 -3 -4; -5 -6] := τ[-1 -2; 1 2] * t[1 2 -3 -4; -5 -6] + @planar t2[-1 -2 -3 -4; -5 -6] := ττ[-1 -2; 1 2] * t[1 2 -3 -4; -5 -6] + @planar t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -2 -1] * t[1 2 -3 -4; -5 -6] + @planar t4[-1 -2 -3 -4; -5 -6] := τ'[-2 2; -1 1] * t[1 2 -3 -4; -5 -6] + t5 = braid(t, ((2, 1, 3, 4), (5, 6)), (1, 2, 3, 4, 5, 6)) + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + @test t1 ≈ t5 + + ττ = copy(BraidingTensor(V1', V1')) + @planar t1[-1 -2 -3 -4; -5 -6] := τ[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6] + @planar t2[-1 -2 -3 -4; -5 -6] := ττ[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6] + @planar t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -3 -2] * t[-1 1 2 -4; -5 -6] + @planar t4[-1 -2 -3 -4; -5 -6] := τ'[-3 2; -2 1] * t[-1 1 2 -4; -5 -6] + t5 = braid(t, ((1, 3, 2, 4), (5, 6)), (1, 2, 3, 4, 5, 6)) + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + @test t1 ≈ t5 + + ττ = copy(BraidingTensor(V1', V1)) + @planar t1[-1 -2 -3 -4; -5 -6] := τ[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6] + @planar t2[-1 -2 -3 -4; -5 -6] := ττ[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6] + @planar t3[-1 -2 -3 -4; -5 -6] := τ[-3 2; -2 1] * t[-1 1 2 -4; -5 -6] + @planar t4[-1 -2 -3 -4; -5 -6] := τ'[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6] + # @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1, V1')) + @planar t1[-1 -2 -3 -4; -5 -6] := τ[-3 2; -2 1] * t[-1 1 2 -4; -5 -6] + @planar t2[-1 -2 -3 -4; -5 -6] := ττ[-3 2; -2 1] * t[-1 1 2 -4; -5 -6] + @planar t3[-1 -2 -3 -4; -5 -6] := τ[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6] + @planar t4[-1 -2 -3 -4; -5 -6] := τ'[2 1; -3 -2] * t[-1 1 2 -4; -5 -6] + # @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1, V1)) + @planar t1[-1 -2 -3 -4; -5 -6] := τ[2 1; -3 -2] * t[-1 1 2 -4; -5 -6] + @planar t2[-1 -2 -3 -4; -5 -6] := ττ[2 1; -3 -2] * t[-1 1 2 -4; -5 -6] + @planar t3[-1 -2 -3 -4; -5 -6] := τ[-2 -3; 1 2] * t[-1 1 2 -4; -5 -6] + @planar t4[-1 -2 -3 -4; -5 -6] := τ'[1 -2; 2 -3] * t[-1 1 2 -4; -5 -6] + # @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1', V1)) + @planar t1[-1 -2 -3 -4; -5 -6] := τ[-3 -4; 1 2] * t[-1 -2 1 2; -5 -6] + @planar t2[-1 -2 -3 -4; -5 -6] := ττ[-3 -4; 1 2] * t[-1 -2 1 2; -5 -6] + @planar t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -4 -3] * t[-1 -2 1 2; -5 -6] + @planar t4[-1 -2 -3 -4; -5 -6] := τ'[-4 2; -3 1] * t[-1 -2 1 2; -5 -6] + # @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1', V1)) + @planar t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ[1 2; -5 -6] + @planar t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * ττ[1 2; -5 -6] + @planar t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ[-6 -5; 2 1] + @planar t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ'[2 -6; 1 -5] + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1, V1')) + @planar t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ'[1 2; -5 -6] + @planar t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * ττ'[1 2; -5 -6] + @planar t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ'[-6 -5; 2 1] + @planar t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2] * τ[2 -6; 1 -5] + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1, V1)) + @planar t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * τ[-4 -6; 1 2] + @planar t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * ττ[-4 -6; 1 2] + @planar t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * τ[2 1; -6 -4] + @planar t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2] * τ'[-6 2; -4 1] + # @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1', V1)) + @planar t1[(); (-1, -2)] := τ[2 1; 3 4] * t[1 2 3 4; -1 -2] + @planar t2[(); (-1, -2)] := ττ[2 1; 3 4] * t[1 2 3 4; -1 -2] + @planar t3[(); (-1, -2)] := τ[4 3; 1 2] * t[1 2 3 4; -1 -2] + @planar t4[(); (-1, -2)] := τ'[1 4; 2 3] * t[1 2 3 4; -1 -2] + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1, V1)) + @planar t1[-1; -2] := τ[2 1; 3 4] * t[-1 1 2 3; -2 4] + @planar t2[-1; -2] := ττ[2 1; 3 4] * t[-1 1 2 3; -2 4] + @planar t3[-1; -2] := τ[4 3; 1 2] * t[-1 1 2 3; -2 4] + @planar t4[-1; -2] := τ'[1 4; 2 3] * t[-1 1 2 3; -2 4] + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1, V1')) + @planar t1[-1 -2] := τ[2 1; 3 4] * t[-1 -2 1 2; 4 3] + @planar t2[-1 -2] := ττ[2 1; 3 4] * t[-1 -2 1 2; 4 3] + @planar t3[-1 -2] := τ[4 3; 1 2] * t[-1 -2 1 2; 4 3] + @planar t4[-1 -2] := τ'[1 4; 2 3] * t[-1 -2 1 2; 4 3] + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1, V1')) + @planar t1[-1 -2; -3 -4] := τ[-1 3; 1 2] * t[1 2 3 -2; -3 -4] + @planar t2[-1 -2; -3 -4] := ττ[-1 3; 1 2] * t[1 2 3 -2; -3 -4] + @planar t3[-1 -2; -3 -4] := τ[2 1; 3 -1] * t[1 2 3 -2; -3 -4] + @planar t4[-1 -2; -3 -4] := τ'[3 2; -1 1] * t[1 2 3 -2; -3 -4] + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1', V1')) + @planar t1[-1 -2; -3 -4] := τ'[-2 3; 1 2] * t[-1 1 2 3; -3 -4] + @planar t2[-1 -2; -3 -4] := ττ'[-2 3; 1 2] * t[-1 1 2 3; -3 -4] + @planar t3[-1 -2; -3 -4] := τ'[2 1; 3 -2] * t[-1 1 2 3; -3 -4] + @planar t4[-1 -2; -3 -4] := τ[3 2; -2 1] * t[-1 1 2 3; -3 -4] + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1', V1)) + @planar t1[-1 -2 -3; -4] := τ[-3 3; 1 2] * t[-1 -2 1 2; -4 3] + @planar t2[-1 -2 -3; -4] := ττ[-3 3; 1 2] * t[-1 -2 1 2; -4 3] + @planar t3[-1 -2 -3; -4] := τ[2 1; 3 -3] * t[-1 -2 1 2; -4 3] + @planar t4[-1 -2 -3; -4] := τ'[3 2; -3 1] * t[-1 -2 1 2; -4 3] + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1', V1)) + @planar t1[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ[1 2; -4 3] + @planar t2[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * ττ[1 2; -4 3] + @planar t3[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ[3 -4; 2 1] + @planar t4[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ'[2 3; 1 -4] + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 + + ττ = copy(BraidingTensor(V1, V1')) + @planar t1[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ'[1 2; -4 3] + @planar t2[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * ττ'[1 2; -4 3] + @planar t3[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ'[3 -4; 2 1] + @planar t4[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2] * τ[2 3; 1 -4] + @test t1 ≈ t2 + @test t1 ≈ t3 + @test t1 ≈ t4 end @testset "@planar" verbose = true begin diff --git a/test/runtests.jl b/test/runtests.jl index 92ceb7ac..63a1ed76 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,70 +1,26 @@ -using Test -using TestExtras -using Random -using TensorKit -using Combinatorics -using TensorKit: ProductSector, fusiontensor, pentagon_equation, hexagon_equation -using TensorOperations -using Base.Iterators: take, product -# using SUNRepresentations: SUNIrrep -# const SU3Irrep = SUNIrrep{3} -using LinearAlgebra: LinearAlgebra +using Test, Random, TensorKit -include("newsectors.jl") -using .NewSectors +include("choosetests.jl") -const TK = TensorKit +choices = choosetests(ARGS) -Random.seed!(1234) +tests = choices.tests +sectors = choices.sectors +exit_on_error = choices.exit_on_error +use_revise = choices.use_revise +seed = choices.seed -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 +module Utility +include("utility.jl") 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 +@eval(Utility, const sectorlist = eval.(Meta.parse.($sectors))) + +for t in tests + modname = Symbol("Test_($t)") + m = @eval(Main, module $modname end) + @eval(m, using Test, TestExtras, Random, TensorKit, ..Utility) + !isnothing(seed) && Random.seed!(seed) + @testset "$t" verbose = true begin + Base.include(m, "$t.jl") end end - -sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewSU2Irrep, # SU3Irrep, - FibonacciAnyon, IsingAnyon, FermionParity, FermionParity ⊠ FermionParity, - Z3Irrep ⊠ Z4Irrep, FermionParity ⊠ U1Irrep ⊠ SU2Irrep, - FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, NewSU2Irrep ⊠ NewSU2Irrep, - NewSU2Irrep ⊠ SU2Irrep, FermionParity ⊠ SU2Irrep ⊠ NewSU2Irrep, - Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon) - -Ti = time() -include("sectors.jl") -include("fusiontrees.jl") -include("spaces.jl") -include("tensors.jl") -include("planar.jl") -include("ad.jl") -Tf = time() -printstyled("Finished all tests in ", - string(round((Tf - Ti) / 60; sigdigits=3)), - " minutes."; bold=true, color=Base.info_color()) -println() diff --git a/test/sectors.jl b/test/sectors.jl index fbda1e2f..a628953f 100644 --- a/test/sectors.jl +++ b/test/sectors.jl @@ -1,53 +1,117 @@ 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] - @constinferred ⊗(s..., s...) + +using TensorKit: pentagon_equation, hexagon_equation + +@testset "$(TensorKit.type_repr(I))" for I in sectorlist + @testset "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...) + BraidingStyle(I) isa NoBraiding || @constinferred Rsymbol(s...) + @constinferred Bsymbol(s...) + @constinferred Fsymbol(s..., s...) + it = @constinferred s[1] ⊗ s[2] + @constinferred ⊗(s..., s...) + end + + @testset "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 - @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 + + @testset "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 - sprev = s - i >= 10 && break + @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12) 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 + + @testset "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 + if !(BraidingStyle(I) isa NoBraiding) + @testset "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 + + if hasfusiontensor(I) + @testset "Fusion tensor and F-move" begin + 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 - if hasfusiontensor(I) - @testset "Sector $I: fusion tensor and F-move and R-move" begin + if !(BraidingStyle(I) isa NoBraiding) + @testset "Fusion tensor 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)) @@ -58,60 +122,6 @@ println("------------------------------------") @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 diff --git a/test/spaces.jl b/test/spaces.jl index 3634ef6c..1dc5dfe9 100644 --- a/test/spaces.jl +++ b/test/spaces.jl @@ -1,407 +1,406 @@ println("------------------------------------") -println("| Fields and vector spaces |") +println("Spaces") println("------------------------------------") -@timedtestset "Fields and vector spaces" verbose = true begin - @timedtestset "Fields" begin - @test isa(ℝ, Field) - @test isa(ℂ, Field) - @test eval(Meta.parse(sprint(show, ℝ))) == ℝ - @test eval(Meta.parse(sprint(show, ℂ))) == ℂ - @test ℝ ⊆ ℝ - @test ℝ ⊆ ℂ - @test ℂ ⊆ ℂ - @test !(ℂ ⊆ ℝ) - for T in (Int8, Int16, Int32, Int64, BigInt) - @test one(T) ∈ ℝ - @test one(Rational{T}) ∈ ℝ - @test !(one(Complex{T}) ∈ ℝ) - @test !(one(Complex{Rational{T}}) ∈ ℝ) - @test one(T) ∈ ℂ - @test one(Rational{T}) ∈ ℂ - @test one(Complex{T}) ∈ ℂ - @test one(Complex{Rational{T}} ∈ ℂ) +@testset "Fields" begin + @test isa(ℝ, Field) + @test isa(ℂ, Field) + @test eval(Meta.parse(sprint(show, ℝ))) == ℝ + @test eval(Meta.parse(sprint(show, ℂ))) == ℂ + @test ℝ ⊆ ℝ + @test ℝ ⊆ ℂ + @test ℂ ⊆ ℂ + @test !(ℂ ⊆ ℝ) - @test T ⊆ ℝ - @test Rational{T} ⊆ ℝ - @test !(Complex{T} ⊆ ℝ) - @test !(Complex{Rational{T}} ⊆ ℝ) - @test T ⊆ ℂ - @test Rational{T} ⊆ ℂ - @test Complex{T} ⊆ ℂ - @test Complex{Rational{T}} ⊆ ℂ - end - for T in (Float32, Float64, BigFloat) - @test one(T) ∈ ℝ - @test !(one(Complex{T}) ∈ ℝ) - @test one(T) ∈ ℂ - @test one(Complex{T} ∈ ℂ) + for T in (Int8, Int16, Int32, Int64, BigInt) + @test one(T) ∈ ℝ + @test one(Rational{T}) ∈ ℝ + @test !(one(Complex{T}) ∈ ℝ) + @test !(one(Complex{Rational{T}}) ∈ ℝ) + @test one(T) ∈ ℂ + @test one(Rational{T}) ∈ ℂ + @test one(Complex{T}) ∈ ℂ + @test one(Complex{Rational{T}} ∈ ℂ) - @test T ⊆ ℝ - @test !(Complex{T} ⊆ ℝ) - @test T ⊆ ℂ - @test Complex{T} ⊆ ℂ - end + @test T ⊆ ℝ + @test Rational{T} ⊆ ℝ + @test !(Complex{T} ⊆ ℝ) + @test !(Complex{Rational{T}} ⊆ ℝ) + @test T ⊆ ℂ + @test Rational{T} ⊆ ℂ + @test Complex{T} ⊆ ℂ + @test Complex{Rational{T}} ⊆ ℂ end + for T in (Float32, Float64, BigFloat) + @test one(T) ∈ ℝ + @test !(one(Complex{T}) ∈ ℝ) + @test one(T) ∈ ℂ + @test one(Complex{T} ∈ ℂ) - @timedtestset "ElementarySpace: CartesianSpace" begin - d = 2 - V = ℝ^d - @test eval(Meta.parse(sprint(show, V))) == V - @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V) - @test isa(V, VectorSpace) - @test isa(V, ElementarySpace) - @test isa(InnerProductStyle(V), HasInnerProduct) - @test isa(InnerProductStyle(V), EuclideanProduct) - @test isa(V, CartesianSpace) - @test !isdual(V) - @test !isdual(V') - @test V == CartesianSpace(Trivial() => d) == CartesianSpace(Dict(Trivial() => d)) - @test @constinferred(hash(V)) == hash(deepcopy(V)) - @test V == @constinferred(dual(V)) == @constinferred(conj(V)) == - @constinferred(adjoint(V)) - @test field(V) == ℝ - @test @constinferred(sectortype(V)) == Trivial - @test ((@constinferred sectors(V))...,) == (Trivial(),) - @test length(sectors(V)) == 1 - @test @constinferred(TensorKit.hassector(V, Trivial())) - @test @constinferred(dim(V)) == d == @constinferred(dim(V, Trivial())) - @test dim(@constinferred(typeof(V)())) == 0 - @test (sectors(typeof(V)())...,) == () - @test @constinferred(TensorKit.axes(V)) == Base.OneTo(d) - @test ℝ^d == ℝ[](d) == CartesianSpace(d) == typeof(V)(d) - W = @constinferred ℝ^1 - @test @constinferred(oneunit(V)) == W == oneunit(typeof(V)) - @test @constinferred(⊕(V, V)) == ℝ^(2d) - @test @constinferred(⊕(V, oneunit(V))) == ℝ^(d + 1) - @test @constinferred(⊕(V, V, V, V)) == ℝ^(4d) - @test @constinferred(fuse(V, V)) == ℝ^(d^2) - @test @constinferred(fuse(V, V', V, V')) == ℝ^(d^4) - @test @constinferred(flip(V)) == V' - @test flip(V) ≅ V - @test flip(V) ≾ V - @test flip(V) ≿ V - @test V ≺ ⊕(V, V) - @test !(V ≻ ⊕(V, V)) - @test @constinferred(infimum(V, ℝ^3)) == V - @test @constinferred(supremum(V', ℝ^3)) == ℝ^3 + @test T ⊆ ℝ + @test !(Complex{T} ⊆ ℝ) + @test T ⊆ ℂ + @test Complex{T} ⊆ ℂ end +end - @timedtestset "ElementarySpace: ComplexSpace" begin - d = 2 - V = ℂ^d - @test eval(Meta.parse(sprint(show, V))) == V - @test eval(Meta.parse(sprint(show, V'))) == V' - @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V) - @test isa(V, VectorSpace) - @test isa(V, ElementarySpace) - @test isa(InnerProductStyle(V), HasInnerProduct) - @test isa(InnerProductStyle(V), EuclideanProduct) - @test isa(V, ComplexSpace) - @test !isdual(V) - @test isdual(V') - @test V == ComplexSpace(Trivial() => d) == ComplexSpace(Dict(Trivial() => d)) - @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V') - @test @constinferred(dual(V)) == @constinferred(conj(V)) == - @constinferred(adjoint(V)) != V - @test @constinferred(field(V)) == ℂ - @test @constinferred(sectortype(V)) == Trivial - @test @constinferred(sectortype(V)) == Trivial - @test ((@constinferred sectors(V))...,) == (Trivial(),) - @test length(sectors(V)) == 1 - @test @constinferred(TensorKit.hassector(V, Trivial())) - @test @constinferred(dim(V)) == d == @constinferred(dim(V, Trivial())) - @test dim(@constinferred(typeof(V)())) == 0 - @test (sectors(typeof(V)())...,) == () - @test @constinferred(TensorKit.axes(V)) == Base.OneTo(d) - @test ℂ^d == Vect[Trivial](d) == Vect[](Trivial() => d) == ℂ[](d) == typeof(V)(d) - W = @constinferred ℂ^1 - @test @constinferred(oneunit(V)) == W == oneunit(typeof(V)) - @test @constinferred(⊕(V, V)) == ℂ^(2d) - @test_throws SpaceMismatch (⊕(V, V')) - promote_except = ErrorException("promotion of types $(typeof(ℝ^d)) and " * - "$(typeof(ℂ^d)) failed to change any arguments") - @test_throws promote_except (⊕(ℝ^d, ℂ^d)) - @test_throws promote_except (⊗(ℝ^d, ℂ^d)) - @test @constinferred(⊕(V, V)) == ℂ^(2d) - @test @constinferred(⊕(V, oneunit(V))) == ℂ^(d + 1) - @test @constinferred(⊕(V, V, V, V)) == ℂ^(4d) - @test @constinferred(fuse(V, V)) == ℂ^(d^2) - @test @constinferred(fuse(V, V', V, V')) == ℂ^(d^4) - @test @constinferred(flip(V)) == V' - @test flip(V) ≅ V - @test flip(V) ≾ V - @test flip(V) ≿ V - @test V ≺ ⊕(V, V) - @test !(V ≻ ⊕(V, V)) - @test @constinferred(infimum(V, ℂ^3)) == V - @test @constinferred(supremum(V', (ℂ^3)')) == dual(ℂ^3) == conj(ℂ^3) - end +@testset "CartesianSpace" begin + d = 2 + V = ℝ^d + @test eval(Meta.parse(sprint(show, V))) == V + @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V) + @test isa(V, VectorSpace) + @test isa(V, ElementarySpace) + @test isa(InnerProductStyle(V), HasInnerProduct) + @test isa(InnerProductStyle(V), EuclideanProduct) + @test isa(V, CartesianSpace) + @test !isdual(V) + @test !isdual(V') + @test V == CartesianSpace(Trivial() => d) == CartesianSpace(Dict(Trivial() => d)) + @test @constinferred(hash(V)) == hash(deepcopy(V)) + @test V == @constinferred(dual(V)) == @constinferred(conj(V)) == + @constinferred(adjoint(V)) + @test field(V) == ℝ + @test @constinferred(sectortype(V)) == Trivial + @test ((@constinferred sectors(V))...,) == (Trivial(),) + @test length(sectors(V)) == 1 + @test @constinferred(TensorKit.hassector(V, Trivial())) + @test @constinferred(dim(V)) == d == @constinferred(dim(V, Trivial())) + @test dim(@constinferred(typeof(V)())) == 0 + @test (sectors(typeof(V)())...,) == () + @test @constinferred(TensorKit.axes(V)) == Base.OneTo(d) + @test ℝ^d == ℝ[](d) == CartesianSpace(d) == typeof(V)(d) + W = @constinferred ℝ^1 + @test @constinferred(oneunit(V)) == W == oneunit(typeof(V)) + @test @constinferred(⊕(V, V)) == ℝ^(2d) + @test @constinferred(⊕(V, oneunit(V))) == ℝ^(d + 1) + @test @constinferred(⊕(V, V, V, V)) == ℝ^(4d) + @test @constinferred(fuse(V, V)) == ℝ^(d^2) + @test @constinferred(fuse(V, V', V, V')) == ℝ^(d^4) + @test @constinferred(flip(V)) == V' + @test flip(V) ≅ V + @test flip(V) ≾ V + @test flip(V) ≿ V + @test V ≺ ⊕(V, V) + @test !(V ≻ ⊕(V, V)) + @test @constinferred(infimum(V, ℝ^3)) == V + @test @constinferred(supremum(V', ℝ^3)) == ℝ^3 +end - @timedtestset "ElementarySpace: GeneralSpace" begin - d = 2 - V = GeneralSpace{ℂ}(d) - @test eval(Meta.parse(sprint(show, V))) == V - @test eval(Meta.parse(sprint(show, dual(V)))) == dual(V) - @test eval(Meta.parse(sprint(show, conj(V)))) == conj(V) - @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V) - @test !isdual(V) - @test isdual(V') - @test !isdual(conj(V)) - @test isdual(conj(V')) - @test !TensorKit.isconj(V) - @test !TensorKit.isconj(V') - @test TensorKit.isconj(conj(V)) - @test TensorKit.isconj(conj(V')) - @test isa(V, VectorSpace) - @test isa(V, ElementarySpace) - @test !isa(InnerProductStyle(V), HasInnerProduct) - @test !isa(InnerProductStyle(V), EuclideanProduct) - @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V') - @test @constinferred(dual(V)) != @constinferred(conj(V)) != V - @test @constinferred(field(V)) == ℂ - @test @constinferred(sectortype(V)) == Trivial - @test @constinferred(TensorKit.hassector(V, Trivial())) - @test @constinferred(dim(V)) == d == @constinferred(dim(V, Trivial())) - @test @constinferred(TensorKit.axes(V)) == Base.OneTo(d) - end +@timedtestset "ElementarySpace: ComplexSpace" begin + d = 2 + V = ℂ^d + @test eval(Meta.parse(sprint(show, V))) == V + @test eval(Meta.parse(sprint(show, V'))) == V' + @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V) + @test isa(V, VectorSpace) + @test isa(V, ElementarySpace) + @test isa(InnerProductStyle(V), HasInnerProduct) + @test isa(InnerProductStyle(V), EuclideanProduct) + @test isa(V, ComplexSpace) + @test !isdual(V) + @test isdual(V') + @test V == ComplexSpace(Trivial() => d) == ComplexSpace(Dict(Trivial() => d)) + @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V') + @test @constinferred(dual(V)) == @constinferred(conj(V)) == + @constinferred(adjoint(V)) != V + @test @constinferred(field(V)) == ℂ + @test @constinferred(sectortype(V)) == Trivial + @test @constinferred(sectortype(V)) == Trivial + @test ((@constinferred sectors(V))...,) == (Trivial(),) + @test length(sectors(V)) == 1 + @test @constinferred(TensorKit.hassector(V, Trivial())) + @test @constinferred(dim(V)) == d == @constinferred(dim(V, Trivial())) + @test dim(@constinferred(typeof(V)())) == 0 + @test (sectors(typeof(V)())...,) == () + @test @constinferred(TensorKit.axes(V)) == Base.OneTo(d) + @test ℂ^d == Vect[Trivial](d) == Vect[](Trivial() => d) == ℂ[](d) == typeof(V)(d) + W = @constinferred ℂ^1 + @test @constinferred(oneunit(V)) == W == oneunit(typeof(V)) + @test @constinferred(⊕(V, V)) == ℂ^(2d) + @test_throws SpaceMismatch (⊕(V, V')) + promote_except = ErrorException("promotion of types $(typeof(ℝ^d)) and " * + "$(typeof(ℂ^d)) failed to change any arguments") + @test_throws promote_except (⊕(ℝ^d, ℂ^d)) + @test_throws promote_except (⊗(ℝ^d, ℂ^d)) + @test @constinferred(⊕(V, V)) == ℂ^(2d) + @test @constinferred(⊕(V, oneunit(V))) == ℂ^(d + 1) + @test @constinferred(⊕(V, V, V, V)) == ℂ^(4d) + @test @constinferred(fuse(V, V)) == ℂ^(d^2) + @test @constinferred(fuse(V, V', V, V')) == ℂ^(d^4) + @test @constinferred(flip(V)) == V' + @test flip(V) ≅ V + @test flip(V) ≾ V + @test flip(V) ≿ V + @test V ≺ ⊕(V, V) + @test !(V ≻ ⊕(V, V)) + @test @constinferred(infimum(V, ℂ^3)) == V + @test @constinferred(supremum(V', (ℂ^3)')) == dual(ℂ^3) == conj(ℂ^3) +end - @timedtestset "ElementarySpace: $(TensorKit.type_repr(Vect[I]))" for I in sectorlist - if Base.IteratorSize(values(I)) === Base.IsInfinite() - set = unique(vcat(one(I), [randsector(I) for k in 1:10])) - gen = (c => 2 for c in set) - else - gen = (values(I)[k] => (k + 1) for k in 1:length(values(I))) - end - V = GradedSpace(gen) - @test eval(Meta.parse(TensorKit.type_repr(typeof(V)))) == typeof(V) - @test eval(Meta.parse(sprint(show, V))) == V - @test eval(Meta.parse(sprint(show, V'))) == V' - @test V' == GradedSpace(gen; dual=true) - @test V == @constinferred GradedSpace(gen...) - @test V' == @constinferred GradedSpace(gen...; dual=true) - @test V == @constinferred GradedSpace(tuple(gen...)) - @test V' == @constinferred GradedSpace(tuple(gen...); dual=true) - @test V == @constinferred GradedSpace(Dict(gen)) - @test V' == @constinferred GradedSpace(Dict(gen); dual=true) - @test V == @inferred Vect[I](gen) - @test V' == @constinferred Vect[I](gen; dual=true) - @test V == @constinferred Vect[I](gen...) - @test V' == @constinferred Vect[I](gen...; dual=true) - @test V == @constinferred Vect[I](Dict(gen)) - @test V' == @constinferred Vect[I](Dict(gen); dual=true) - @test V == @constinferred typeof(V)(c => dim(V, c) for c in sectors(V)) - if I isa ZNIrrep - @test V == @constinferred typeof(V)(V.dims) - @test V' == @constinferred typeof(V)(V.dims; dual=true) - end - @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V') - @test V == GradedSpace(reverse(collect(gen))...) - @test eval(Meta.parse(sprint(show, V))) == V - @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V) - # space with no sectors - @test dim(@constinferred(typeof(V)())) == 0 - # space with a single sector - W = @constinferred GradedSpace(one(I) => 1) - @test W == GradedSpace(one(I) => 1, randsector(I) => 0) - @test @constinferred(oneunit(V)) == W == oneunit(typeof(V)) - # randsector never returns trivial sector, so this cannot error - @test_throws ArgumentError GradedSpace(one(I) => 1, randsector(I) => 0, one(I) => 3) - @test eval(Meta.parse(sprint(show, W))) == W - @test isa(V, VectorSpace) - @test isa(V, ElementarySpace) - @test isa(InnerProductStyle(V), HasInnerProduct) - @test isa(InnerProductStyle(V), EuclideanProduct) - @test isa(V, GradedSpace) - @test isa(V, GradedSpace{I}) - @test @constinferred(dual(V)) == @constinferred(conj(V)) == - @constinferred(adjoint(V)) != V - @test @constinferred(field(V)) == ℂ - @test @constinferred(sectortype(V)) == I - slist = @constinferred sectors(V) - @test @constinferred(TensorKit.hassector(V, first(slist))) - @test @constinferred(dim(V)) == sum(dim(s) * dim(V, s) for s in slist) - @constinferred dim(V, first(slist)) - if hasfusiontensor(I) - @test @constinferred(TensorKit.axes(V)) == Base.OneTo(dim(V)) - end - @test @constinferred(⊕(V, V)) == Vect[I](c => 2dim(V, c) for c in sectors(V)) - @test @constinferred(⊕(V, V, V, V)) == Vect[I](c => 4dim(V, c) for c in sectors(V)) - @test @constinferred(⊕(V, oneunit(V))) == - Vect[I](c => isone(c) + dim(V, c) for c in sectors(V)) - @test @constinferred(fuse(V, oneunit(V))) == V - d = Dict{I,Int}() - for a in sectors(V), b in sectors(V) - for c in a ⊗ b - d[c] = get(d, c, 0) + dim(V, a) * dim(V, b) * Nsymbol(a, b, c) - end +@testset "GeneralSpace" begin + d = 2 + V = GeneralSpace{ℂ}(d) + @test eval(Meta.parse(sprint(show, V))) == V + @test eval(Meta.parse(sprint(show, dual(V)))) == dual(V) + @test eval(Meta.parse(sprint(show, conj(V)))) == conj(V) + @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V) + @test !isdual(V) + @test isdual(V') + @test !isdual(conj(V)) + @test isdual(conj(V')) + @test !TensorKit.isconj(V) + @test !TensorKit.isconj(V') + @test TensorKit.isconj(conj(V)) + @test TensorKit.isconj(conj(V')) + @test isa(V, VectorSpace) + @test isa(V, ElementarySpace) + @test !isa(InnerProductStyle(V), HasInnerProduct) + @test !isa(InnerProductStyle(V), EuclideanProduct) + @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V') + @test @constinferred(dual(V)) != @constinferred(conj(V)) != V + @test @constinferred(field(V)) == ℂ + @test @constinferred(sectortype(V)) == Trivial + @test @constinferred(TensorKit.hassector(V, Trivial())) + @test @constinferred(dim(V)) == d == @constinferred(dim(V, Trivial())) + @test @constinferred(TensorKit.axes(V)) == Base.OneTo(d) +end + +@testset "$(TensorKit.type_repr(Vect[I]))" for I in sectorlist + if Base.IteratorSize(values(I)) === Base.IsInfinite() + set = unique(vcat(one(I), [randsector(I) for k in 1:10])) + gen = (c => 2 for c in set) + else + gen = (values(I)[k] => (k + 1) for k in 1:length(values(I))) + end + V = GradedSpace(gen) + @test eval(Meta.parse(TensorKit.type_repr(typeof(V)))) == typeof(V) + @test eval(Meta.parse(sprint(show, V))) == V + @test eval(Meta.parse(sprint(show, V'))) == V' + @test V' == GradedSpace(gen; dual=true) + @test V == @constinferred GradedSpace(gen...) + @test V' == @constinferred GradedSpace(gen...; dual=true) + @test V == @constinferred GradedSpace(tuple(gen...)) + @test V' == @constinferred GradedSpace(tuple(gen...); dual=true) + @test V == @constinferred GradedSpace(Dict(gen)) + @test V' == @constinferred GradedSpace(Dict(gen); dual=true) + @test V == @inferred Vect[I](gen) + @test V' == @constinferred Vect[I](gen; dual=true) + @test V == @constinferred Vect[I](gen...) + @test V' == @constinferred Vect[I](gen...; dual=true) + @test V == @constinferred Vect[I](Dict(gen)) + @test V' == @constinferred Vect[I](Dict(gen); dual=true) + @test V == @constinferred typeof(V)(c => dim(V, c) for c in sectors(V)) + if I isa ZNIrrep + @test V == @constinferred typeof(V)(V.dims) + @test V' == @constinferred typeof(V)(V.dims; dual=true) + end + @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V') + @test V == GradedSpace(reverse(collect(gen))...) + @test eval(Meta.parse(sprint(show, V))) == V + @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V) + # space with no sectors + @test dim(@constinferred(typeof(V)())) == 0 + # space with a single sector + W = @constinferred GradedSpace(one(I) => 1) + @test W == GradedSpace(one(I) => 1, randsector(I) => 0) + @test @constinferred(oneunit(V)) == W == oneunit(typeof(V)) + # randsector never returns trivial sector, so this cannot error + @test_throws ArgumentError GradedSpace(one(I) => 1, randsector(I) => 0, one(I) => 3) + @test eval(Meta.parse(sprint(show, W))) == W + @test isa(V, VectorSpace) + @test isa(V, ElementarySpace) + @test isa(InnerProductStyle(V), HasInnerProduct) + @test isa(InnerProductStyle(V), EuclideanProduct) + @test isa(V, GradedSpace) + @test isa(V, GradedSpace{I}) + @test @constinferred(dual(V)) == @constinferred(conj(V)) == + @constinferred(adjoint(V)) != V + @test @constinferred(field(V)) == ℂ + @test @constinferred(sectortype(V)) == I + slist = @constinferred sectors(V) + @test @constinferred(TensorKit.hassector(V, first(slist))) + @test @constinferred(dim(V)) == sum(dim(s) * dim(V, s) for s in slist) + @constinferred dim(V, first(slist)) + if hasfusiontensor(I) + @test @constinferred(TensorKit.axes(V)) == Base.OneTo(dim(V)) + end + @test @constinferred(⊕(V, V)) == Vect[I](c => 2dim(V, c) for c in sectors(V)) + @test @constinferred(⊕(V, V, V, V)) == Vect[I](c => 4dim(V, c) for c in sectors(V)) + @test @constinferred(⊕(V, oneunit(V))) == + Vect[I](c => isone(c) + dim(V, c) for c in sectors(V)) + @test @constinferred(fuse(V, oneunit(V))) == V + d = Dict{I,Int}() + for a in sectors(V), b in sectors(V) + for c in a ⊗ b + d[c] = get(d, c, 0) + dim(V, a) * dim(V, b) * Nsymbol(a, b, c) end - @test @constinferred(fuse(V, V)) == GradedSpace(d) - @test @constinferred(flip(V)) == - Vect[I](conj(c) => dim(V, c) for c in sectors(V))' - @test flip(V) ≅ V - @test flip(V) ≾ V - @test flip(V) ≿ V - @test @constinferred(⊕(V, V)) == @constinferred supremum(V, ⊕(V, V)) - @test V == @constinferred infimum(V, ⊕(V, V)) - @test V ≺ ⊕(V, V) - @test !(V ≻ ⊕(V, V)) - @test infimum(V, GradedSpace(one(I) => 3)) == GradedSpace(one(I) => 2) - @test_throws SpaceMismatch (⊕(V, V')) end + @test @constinferred(fuse(V, V)) == GradedSpace(d) + @test @constinferred(flip(V)) == + Vect[I](conj(c) => dim(V, c) for c in sectors(V))' + @test flip(V) ≅ V + @test flip(V) ≾ V + @test flip(V) ≿ V + @test @constinferred(⊕(V, V)) == @constinferred supremum(V, ⊕(V, V)) + @test V == @constinferred infimum(V, ⊕(V, V)) + @test V ≺ ⊕(V, V) + @test !(V ≻ ⊕(V, V)) + @test infimum(V, GradedSpace(one(I) => 3)) == GradedSpace(one(I) => 2) + @test_throws SpaceMismatch (⊕(V, V')) +end - @timedtestset "ProductSpace{ℂ}" begin - V1, V2, V3, V4 = ℂ^1, ℂ^2, ℂ^3, ℂ^4 - P = @constinferred ProductSpace(V1, V2, V3, V4) - @test eval(Meta.parse(sprint(show, P))) == P - @test eval(Meta.parse(sprint(show, typeof(P)))) == typeof(P) - @test isa(P, VectorSpace) - @test isa(P, CompositeSpace) - @test spacetype(P) == ComplexSpace - @test sectortype(P) == Trivial - @test @constinferred(hash(P)) == hash(deepcopy(P)) != hash(P') - @test P == deepcopy(P) - @test P == typeof(P)(P...) - @constinferred (x -> tuple(x...))(P) - @test @constinferred(dual(P)) == P' - @test @constinferred(field(P)) == ℂ - @test @constinferred(*(V1, V2, V3, V4)) == P - @test @constinferred(⊗(V1, V2, V3, V4)) == P - @test @constinferred(⊗(V1, V2 ⊗ V3 ⊗ V4)) == P - @test @constinferred(⊗(V1 ⊗ V2, V3 ⊗ V4)) == P - @test @constinferred(⊗(V1, V2, V3 ⊗ V4)) == P - @test @constinferred(⊗(V1, V2 ⊗ V3, V4)) == P - @test @constinferred(insertunit(P, 3)) == V1 * V2 * oneunit(V1) * V3 * V4 - @test fuse(V1, V2', V3) ≅ V1 ⊗ V2' ⊗ V3 - @test fuse(V1, V2', V3) ≾ V1 ⊗ V2' ⊗ V3 - @test fuse(V1, V2', V3) ≿ V1 ⊗ V2' ⊗ V3 - @test fuse(flip(V1), V2, flip(V3)) ≅ V1 ⊗ V2 ⊗ V3 - @test @constinferred(⊗(P)) == P - @test @constinferred(⊗(V1)) == ProductSpace(V1) - @test eval(Meta.parse(sprint(show, ⊗(V1)))) == ⊗(V1) - @test @constinferred(one(V1)) == @constinferred(one(typeof(V1))) == - @constinferred(one(P)) == @constinferred(one(typeof(P))) == - ProductSpace{ComplexSpace}(()) - @test eval(Meta.parse(sprint(show, one(P)))) == one(P) - @test @constinferred(⊗(one(P), P)) == P - @test @constinferred(⊗(P, one(P))) == P - @test @constinferred(⊗(one(P), one(P))) == one(P) - @test @constinferred(adjoint(P)) == dual(P) == V4' ⊗ V3' ⊗ V2' ⊗ V1' - @test @constinferred(dims(P)) == map(dim, (V1, V2, V3, V4)) - @test @constinferred(dim(P)) == prod(dim, (V1, V2, V3, V4)) - @test @constinferred(dim(P, 2)) == dim(V2) - @test first(@constinferred(sectors(P))) == - (Trivial(), Trivial(), Trivial(), Trivial()) - cube(x) = x^3 - @test @constinferred(cube(V1)) == V1 ⊗ V1 ⊗ V1 - N = 3 - @test V1^N == V1 ⊗ V1 ⊗ V1 - @test P^2 == P ⊗ P - @test @constinferred(dims(P, first(sectors(P)))) == dims(P) - @test ((@constinferred blocksectors(P))...,) == (Trivial(),) - @test isempty(blocksectors(P ⊗ ℂ^0)) - @test isempty(@constinferred(sectors(P ⊗ ℂ^0))) - @test @constinferred(blockdim(P, first(blocksectors(P)))) == dim(P) - @test Base.IteratorEltype(P) == Base.IteratorEltype(typeof(P)) == - Base.IteratorEltype(P.spaces) - @test Base.IteratorSize(P) == Base.IteratorSize(typeof(P)) == - Base.IteratorSize(P.spaces) - @test Base.eltype(P) == Base.eltype(typeof(P)) == typeof(V1) - @test eltype(collect(P)) == typeof(V1) - @test collect(P) == [V1, V2, V3, V4] - end +@testset "ProductSpace{ℂ}" begin + V1, V2, V3, V4 = ℂ^1, ℂ^2, ℂ^3, ℂ^4 + P = @constinferred ProductSpace(V1, V2, V3, V4) + @test eval(Meta.parse(sprint(show, P))) == P + @test eval(Meta.parse(sprint(show, typeof(P)))) == typeof(P) + @test isa(P, VectorSpace) + @test isa(P, CompositeSpace) + @test spacetype(P) == ComplexSpace + @test sectortype(P) == Trivial + @test @constinferred(hash(P)) == hash(deepcopy(P)) != hash(P') + @test P == deepcopy(P) + @test P == typeof(P)(P...) + @constinferred (x -> tuple(x...))(P) + @test @constinferred(dual(P)) == P' + @test @constinferred(field(P)) == ℂ + @test @constinferred(*(V1, V2, V3, V4)) == P + @test @constinferred(⊗(V1, V2, V3, V4)) == P + @test @constinferred(⊗(V1, V2 ⊗ V3 ⊗ V4)) == P + @test @constinferred(⊗(V1 ⊗ V2, V3 ⊗ V4)) == P + @test @constinferred(⊗(V1, V2, V3 ⊗ V4)) == P + @test @constinferred(⊗(V1, V2 ⊗ V3, V4)) == P + @test @constinferred(insertunit(P, 3)) == V1 * V2 * oneunit(V1) * V3 * V4 + @test fuse(V1, V2', V3) ≅ V1 ⊗ V2' ⊗ V3 + @test fuse(V1, V2', V3) ≾ V1 ⊗ V2' ⊗ V3 + @test fuse(V1, V2', V3) ≿ V1 ⊗ V2' ⊗ V3 + @test fuse(flip(V1), V2, flip(V3)) ≅ V1 ⊗ V2 ⊗ V3 + @test @constinferred(⊗(P)) == P + @test @constinferred(⊗(V1)) == ProductSpace(V1) + @test eval(Meta.parse(sprint(show, ⊗(V1)))) == ⊗(V1) + @test @constinferred(one(V1)) == @constinferred(one(typeof(V1))) == + @constinferred(one(P)) == @constinferred(one(typeof(P))) == + ProductSpace{ComplexSpace}(()) + @test eval(Meta.parse(sprint(show, one(P)))) == one(P) + @test @constinferred(⊗(one(P), P)) == P + @test @constinferred(⊗(P, one(P))) == P + @test @constinferred(⊗(one(P), one(P))) == one(P) + @test @constinferred(adjoint(P)) == dual(P) == V4' ⊗ V3' ⊗ V2' ⊗ V1' + @test @constinferred(dims(P)) == map(dim, (V1, V2, V3, V4)) + @test @constinferred(dim(P)) == prod(dim, (V1, V2, V3, V4)) + @test @constinferred(dim(P, 2)) == dim(V2) + @test first(@constinferred(sectors(P))) == + (Trivial(), Trivial(), Trivial(), Trivial()) + cube(x) = x^3 + @test @constinferred(cube(V1)) == V1 ⊗ V1 ⊗ V1 + N = 3 + @test V1^N == V1 ⊗ V1 ⊗ V1 + @test P^2 == P ⊗ P + @test @constinferred(dims(P, first(sectors(P)))) == dims(P) + @test ((@constinferred blocksectors(P))...,) == (Trivial(),) + @test isempty(blocksectors(P ⊗ ℂ^0)) + @test isempty(@constinferred(sectors(P ⊗ ℂ^0))) + @test @constinferred(blockdim(P, first(blocksectors(P)))) == dim(P) + @test Base.IteratorEltype(P) == Base.IteratorEltype(typeof(P)) == + Base.IteratorEltype(P.spaces) + @test Base.IteratorSize(P) == Base.IteratorSize(typeof(P)) == + Base.IteratorSize(P.spaces) + @test Base.eltype(P) == Base.eltype(typeof(P)) == typeof(V1) + @test eltype(collect(P)) == typeof(V1) + @test collect(P) == [V1, V2, V3, V4] +end - @timedtestset "ProductSpace{SU₂Space}" begin - V1, V2, V3 = SU₂Space(0 => 3, 1 // 2 => 1), SU₂Space(0 => 2, 1 => 1), - SU₂Space(1 // 2 => 1, 1 => 1)' - P = @constinferred ProductSpace(V1, V2, V3) - @test eval(Meta.parse(sprint(show, P))) == P - @test eval(Meta.parse(sprint(show, typeof(P)))) == typeof(P) - @test isa(P, VectorSpace) - @test isa(P, CompositeSpace) - @test spacetype(P) == SU₂Space - @test sectortype(P) == Irrep[SU₂] == SU2Irrep - @test @constinferred(hash(P)) == hash(deepcopy(P)) != hash(P') - @test @constinferred(dual(P)) == P' - @test @constinferred(field(P)) == ℂ - @test @constinferred(*(V1, V2, V3)) == P - @test @constinferred(⊗(V1, V2, V3)) == P - @test @constinferred(adjoint(P)) == dual(P) == V3' ⊗ V2' ⊗ V1' - @test @constinferred(insertunit(P, 3; conj=true)) == V1 * V2 * oneunit(V1)' * V3 - @test fuse(V1, V2', V3) ≅ V1 ⊗ V2' ⊗ V3 - @test fuse(V1, V2', V3) ≾ V1 ⊗ V2' ⊗ V3 ≾ fuse(V1 ⊗ V2' ⊗ V3) - @test fuse(V1, V2') ⊗ V3 ≾ V1 ⊗ V2' ⊗ V3 - @test fuse(V1, V2', V3) ≿ V1 ⊗ V2' ⊗ V3 ≿ fuse(V1 ⊗ V2' ⊗ V3) - @test V1 ⊗ fuse(V2', V3) ≿ V1 ⊗ V2' ⊗ V3 - @test fuse(flip(V1) ⊗ V2) ⊗ flip(V3) ≅ V1 ⊗ V2 ⊗ V3 - @test @constinferred(⊗(V1)) == ProductSpace(V1) - @test @constinferred(one(V1)) == @constinferred(one(typeof(V1))) == - @constinferred(one(P)) == @constinferred(one(typeof(P))) == - ProductSpace{ComplexSpace}(()) - @test @constinferred(dims(P)) == map(dim, (V1, V2, V3)) - @test @constinferred(dim(P)) == prod(dim, (V1, V2, V3)) - for s in @constinferred(sectors(P)) - @test hassector(P, s) - @test @constinferred(dims(P, s)) == dim.((V1, V2, V3), s) - end - @test sum(dim(c) * blockdim(P, c) for c in @constinferred(blocksectors(P))) == - dim(P) +@testset "ProductSpace{SU₂Space}" begin + V1, V2, V3 = SU₂Space(0 => 3, 1 // 2 => 1), SU₂Space(0 => 2, 1 => 1), + SU₂Space(1 // 2 => 1, 1 => 1)' + P = @constinferred ProductSpace(V1, V2, V3) + @test eval(Meta.parse(sprint(show, P))) == P + @test eval(Meta.parse(sprint(show, typeof(P)))) == typeof(P) + @test isa(P, VectorSpace) + @test isa(P, CompositeSpace) + @test spacetype(P) == SU₂Space + @test sectortype(P) == Irrep[SU₂] == SU2Irrep + @test @constinferred(hash(P)) == hash(deepcopy(P)) != hash(P') + @test @constinferred(dual(P)) == P' + @test @constinferred(field(P)) == ℂ + @test @constinferred(*(V1, V2, V3)) == P + @test @constinferred(⊗(V1, V2, V3)) == P + @test @constinferred(adjoint(P)) == dual(P) == V3' ⊗ V2' ⊗ V1' + @test @constinferred(insertunit(P, 3; conj=true)) == V1 * V2 * oneunit(V1)' * V3 + @test fuse(V1, V2', V3) ≅ V1 ⊗ V2' ⊗ V3 + @test fuse(V1, V2', V3) ≾ V1 ⊗ V2' ⊗ V3 ≾ fuse(V1 ⊗ V2' ⊗ V3) + @test fuse(V1, V2') ⊗ V3 ≾ V1 ⊗ V2' ⊗ V3 + @test fuse(V1, V2', V3) ≿ V1 ⊗ V2' ⊗ V3 ≿ fuse(V1 ⊗ V2' ⊗ V3) + @test V1 ⊗ fuse(V2', V3) ≿ V1 ⊗ V2' ⊗ V3 + @test fuse(flip(V1) ⊗ V2) ⊗ flip(V3) ≅ V1 ⊗ V2 ⊗ V3 + @test @constinferred(⊗(V1)) == ProductSpace(V1) + @test @constinferred(one(V1)) == @constinferred(one(typeof(V1))) == + @constinferred(one(P)) == @constinferred(one(typeof(P))) == + ProductSpace{ComplexSpace}(()) + @test @constinferred(dims(P)) == map(dim, (V1, V2, V3)) + @test @constinferred(dim(P)) == prod(dim, (V1, V2, V3)) + for s in @constinferred(sectors(P)) + @test hassector(P, s) + @test @constinferred(dims(P, s)) == dim.((V1, V2, V3), s) end + @test sum(dim(c) * blockdim(P, c) for c in @constinferred(blocksectors(P))) == + dim(P) +end - @timedtestset "Deligne tensor product of spaces" begin - V1 = SU₂Space(0 => 3, 1 // 2 => 1) - V2 = SU₂Space(0 => 2, 1 => 1)' - V3 = ℤ₃Space(0 => 3, 1 => 2, 2 => 1) - V4 = ℂ^3 +@testset "Deligne tensor product of spaces" begin + V1 = SU₂Space(0 => 3, 1 // 2 => 1) + V2 = SU₂Space(0 => 2, 1 => 1)' + V3 = ℤ₃Space(0 => 3, 1 => 2, 2 => 1) + V4 = ℂ^3 - for W1 in (V1, V2, V3, V4) - for W2 in (V1, V2, V3, V4) - for W3 in (V1, V2, V3, V4) - for W4 in (V1, V2, V3, V4) - Ws = @constinferred(W1 ⊠ W2 ⊠ W3 ⊠ W4) - @test Ws == @constinferred((W1 ⊠ W2) ⊠ (W3 ⊠ W4)) == - @constinferred(((W1 ⊠ W2) ⊠ W3) ⊠ W4) == - @constinferred((W1 ⊠ (W2 ⊠ W3)) ⊠ W4) == - @constinferred(W1 ⊠ ((W2 ⊠ W3)) ⊠ W4) == - @constinferred(W1 ⊠ (W2 ⊠ (W3 ⊠ W4))) - I1, I2, I3, I4 = map(sectortype, (W1, W2, W3, W4)) - I = sectortype(Ws) - @test I == @constinferred((I1 ⊠ I2) ⊠ (I3 ⊠ I4)) == - @constinferred(((I1 ⊠ I2) ⊠ I3) ⊠ I4) == - @constinferred((I1 ⊠ (I2 ⊠ I3)) ⊠ I4) == - @constinferred(I1 ⊠ ((I2 ⊠ I3)) ⊠ I4) == - @constinferred(I1 ⊠ (I2 ⊠ (I3 ⊠ I4))) - @test dim(Ws) == dim(W1) * dim(W2) * dim(W3) * dim(W4) - end + for W1 in (V1, V2, V3, V4) + for W2 in (V1, V2, V3, V4) + for W3 in (V1, V2, V3, V4) + for W4 in (V1, V2, V3, V4) + Ws = @constinferred(W1 ⊠ W2 ⊠ W3 ⊠ W4) + @test Ws == @constinferred((W1 ⊠ W2) ⊠ (W3 ⊠ W4)) == + @constinferred(((W1 ⊠ W2) ⊠ W3) ⊠ W4) == + @constinferred((W1 ⊠ (W2 ⊠ W3)) ⊠ W4) == + @constinferred(W1 ⊠ ((W2 ⊠ W3)) ⊠ W4) == + @constinferred(W1 ⊠ (W2 ⊠ (W3 ⊠ W4))) + I1, I2, I3, I4 = map(sectortype, (W1, W2, W3, W4)) + I = sectortype(Ws) + @test I == @constinferred((I1 ⊠ I2) ⊠ (I3 ⊠ I4)) == + @constinferred(((I1 ⊠ I2) ⊠ I3) ⊠ I4) == + @constinferred((I1 ⊠ (I2 ⊠ I3)) ⊠ I4) == + @constinferred(I1 ⊠ ((I2 ⊠ I3)) ⊠ I4) == + @constinferred(I1 ⊠ (I2 ⊠ (I3 ⊠ I4))) + @test dim(Ws) == dim(W1) * dim(W2) * dim(W3) * dim(W4) end end end - @test sectortype(@constinferred((V1 ⊗ V2) ⊠ V3)) == @constinferred(Irrep[SU₂ × ℤ₃]) - @test dim((V1 ⊗ V2) ⊠ V3) == dim(V1) * dim(V2) * dim(V3) - @test sectortype((V1 ⊗ V2) ⊠ V3 ⊠ V4) == Irrep[SU₂ × ℤ₃] - @test dim((V1 ⊗ V2) ⊠ V3 ⊠ V4) == dim(V1) * dim(V2) * dim(V3) * dim(V4) - @test fuse(V2 ⊠ V4) == fuse(V4 ⊠ V2) == SU₂Space(0 => 6, 1 => 3) - @test fuse(V3 ⊠ V4) == fuse(V4 ⊠ V3) == ℤ₃Space(0 => 9, 1 => 6, 2 => 3) end + @test sectortype(@constinferred((V1 ⊗ V2) ⊠ V3)) == @constinferred(Irrep[SU₂ × ℤ₃]) + @test dim((V1 ⊗ V2) ⊠ V3) == dim(V1) * dim(V2) * dim(V3) + @test sectortype((V1 ⊗ V2) ⊠ V3 ⊠ V4) == Irrep[SU₂ × ℤ₃] + @test dim((V1 ⊗ V2) ⊠ V3 ⊠ V4) == dim(V1) * dim(V2) * dim(V3) * dim(V4) + @test fuse(V2 ⊠ V4) == fuse(V4 ⊠ V2) == SU₂Space(0 => 6, 1 => 3) + @test fuse(V3 ⊠ V4) == fuse(V4 ⊠ V3) == ℤ₃Space(0 => 9, 1 => 6, 2 => 3) +end - @timedtestset "HomSpace" begin - V1, V2, V3, V4, V5 = SU₂Space(0 => 3, 1 // 2 => 1), SU₂Space(0 => 2, 1 => 1), - SU₂Space(1 // 2 => 1, 1 => 1)', SU₂Space(0 => 2, 1 // 2 => 2), - SU₂Space(0 => 1, 1 // 2 => 1, 3 // 2 => 1)' - W = TensorKit.HomSpace(V1 ⊗ V2, V3 ⊗ V4 ⊗ V5) - @test W == (V3 ⊗ V4 ⊗ V5 → V1 ⊗ V2) - @test W == (V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5) - @test W' == (V1 ⊗ V2 → V3 ⊗ V4 ⊗ V5) - @test eval(Meta.parse(sprint(show, W))) == W - @test eval(Meta.parse(sprint(show, typeof(W)))) == typeof(W) - @test spacetype(W) == SU₂Space - @test sectortype(W) == Irrep[SU₂] - @test W[1] == V1 - @test W[2] == V2 - @test W[3] == V3' - @test W[4] == V4' - @test W[5] == V5' - @test @constinferred(hash(W)) == hash(deepcopy(W)) != hash(W') - @test W == deepcopy(W) - end +@testset "HomSpace" begin + V1, V2, V3, V4, V5 = SU₂Space(0 => 3, 1 // 2 => 1), SU₂Space(0 => 2, 1 => 1), + SU₂Space(1 // 2 => 1, 1 => 1)', SU₂Space(0 => 2, 1 // 2 => 2), + SU₂Space(0 => 1, 1 // 2 => 1, 3 // 2 => 1)' + W = TensorKit.HomSpace(V1 ⊗ V2, V3 ⊗ V4 ⊗ V5) + @test W == (V3 ⊗ V4 ⊗ V5 → V1 ⊗ V2) + @test W == (V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5) + @test W' == (V1 ⊗ V2 → V3 ⊗ V4 ⊗ V5) + @test eval(Meta.parse(sprint(show, W))) == W + @test eval(Meta.parse(sprint(show, typeof(W)))) == typeof(W) + @test spacetype(W) == SU₂Space + @test sectortype(W) == Irrep[SU₂] + @test W[1] == V1 + @test W[2] == V2 + @test W[3] == V3' + @test W[4] == V4' + @test W[5] == V5' + @test @constinferred(hash(W)) == hash(deepcopy(W)) != hash(W') + @test W == deepcopy(W) end diff --git a/test/tensors.jl b/test/tensors.jl index 3d7fd3a6..4e3a53a6 100644 --- a/test/tensors.jl +++ b/test/tensors.jl @@ -1,602 +1,541 @@ -Vtr = (ℂ^3, - (ℂ^4)', - ℂ^5, - ℂ^6, - (ℂ^7)') -Vℤ₂ = (ℂ[Z2Irrep](0 => 1, 1 => 1), - ℂ[Z2Irrep](0 => 1, 1 => 2)', - ℂ[Z2Irrep](0 => 3, 1 => 2)', - ℂ[Z2Irrep](0 => 2, 1 => 3), - ℂ[Z2Irrep](0 => 2, 1 => 5)) -Vfℤ₂ = (ℂ[FermionParity](0 => 1, 1 => 1), - ℂ[FermionParity](0 => 1, 1 => 2)', - ℂ[FermionParity](0 => 3, 1 => 2)', - ℂ[FermionParity](0 => 2, 1 => 3), - ℂ[FermionParity](0 => 2, 1 => 5)) -Vℤ₃ = (ℂ[Z3Irrep](0 => 1, 1 => 2, 2 => 2), - ℂ[Z3Irrep](0 => 3, 1 => 1, 2 => 1), - ℂ[Z3Irrep](0 => 2, 1 => 2, 2 => 1)', - ℂ[Z3Irrep](0 => 1, 1 => 2, 2 => 3), - ℂ[Z3Irrep](0 => 1, 1 => 3, 2 => 3)') -VU₁ = (ℂ[U1Irrep](0 => 1, 1 => 2, -1 => 2), - ℂ[U1Irrep](0 => 3, 1 => 1, -1 => 1), - ℂ[U1Irrep](0 => 2, 1 => 2, -1 => 1)', - ℂ[U1Irrep](0 => 1, 1 => 2, -1 => 3), - ℂ[U1Irrep](0 => 1, 1 => 3, -1 => 3)') -VfU₁ = (ℂ[FermionNumber](0 => 1, 1 => 2, -1 => 2), - ℂ[FermionNumber](0 => 3, 1 => 1, -1 => 1), - ℂ[FermionNumber](0 => 2, 1 => 2, -1 => 1)', - ℂ[FermionNumber](0 => 1, 1 => 2, -1 => 3), - ℂ[FermionNumber](0 => 1, 1 => 3, -1 => 3)') -VCU₁ = (ℂ[CU1Irrep]((0, 0) => 1, (0, 1) => 2, 1 => 1), - ℂ[CU1Irrep]((0, 0) => 3, (0, 1) => 0, 1 => 1), - ℂ[CU1Irrep]((0, 0) => 1, (0, 1) => 0, 1 => 2)', - ℂ[CU1Irrep]((0, 0) => 2, (0, 1) => 2, 1 => 1), - ℂ[CU1Irrep]((0, 0) => 2, (0, 1) => 1, 1 => 2)') -VSU₂ = (ℂ[SU2Irrep](0 => 3, 1 // 2 => 1), - ℂ[SU2Irrep](0 => 2, 1 => 1), - ℂ[SU2Irrep](1 // 2 => 1, 1 => 1)', - ℂ[SU2Irrep](0 => 2, 1 // 2 => 2), - ℂ[SU2Irrep](0 => 1, 1 // 2 => 1, 3 // 2 => 1)') -VfSU₂ = (ℂ[FermionSpin](0 => 3, 1 // 2 => 1), - ℂ[FermionSpin](0 => 2, 1 => 1), - ℂ[FermionSpin](1 // 2 => 1, 1 => 1)', - ℂ[FermionSpin](0 => 2, 1 // 2 => 2), - ℂ[FermionSpin](0 => 1, 1 // 2 => 1, 3 // 2 => 1)') -# VSU₃ = (ℂ[SU3Irrep]((0, 0, 0) => 3, (1, 0, 0) => 1), -# ℂ[SU3Irrep]((0, 0, 0) => 3, (2, 0, 0) => 1)', -# ℂ[SU3Irrep]((1, 1, 0) => 1, (2, 1, 0) => 1), -# ℂ[SU3Irrep]((1, 0, 0) => 1, (2, 0, 0) => 1), -# ℂ[SU3Irrep]((0, 0, 0) => 1, (1, 0, 0) => 1, (1, 1, 0) => 1)') - -for V in (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂)#, VSU₃) +println("------------------------------------") +println("Tensors") +println("------------------------------------") + +using Combinatorics: permutations + +@testset "$(TensorKit.type_repr(I))" verbose = true for I in sectorlist + V = smallspace(I) + if isnothing(V) + "No spaces defined for $(TensorKit.type_repr(I)), skipping tests" + continue + end + Istr = TensorKit.type_repr(I) + println("Starting tests for $Istr...") + V1, V2, V3, V4, V5 = V - @assert V3 * V4 * V2 ≿ V1' * V5' # necessary for leftorth tests - @assert V3 * V4 ≾ V1' * V2' * V5' # necessary for rightorth tests -end + @assert V3 * V4 * V2 ≿ V1' * V5' "leftorth tests assumption" + @assert V3 * V4 ≾ V1' * V2' * V5' "rightorth tests assumption" -spacelist = try - if ENV["CI"] == "true" - println("Detected running on CI") - if Sys.iswindows() - (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂) - elseif Sys.isapple() - (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VfU₁, VfSU₂)#, VSU₃) - else - (Vtr, Vℤ₂, Vfℤ₂, VU₁, VCU₁, VSU₂, VfSU₂)#, VSU₃) - end - else - (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂)#, VSU₃) + @testset "basic tensor properties ($T)" for T in (Int, Float32, Float64, ComplexF32, + ComplexF64, BigFloat) + W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5 + t = Tensor(zeros, T, W) + @test @constinferred(hash(t)) == hash(deepcopy(t)) + @test scalartype(t) == T + @test norm(t) == 0 + @test codomain(t) == W + @test space(t) == (W ← one(W)) + @test domain(t) == one(W) + @test typeof(t) == @constinferred tensormaptype(spacetype(t), 5, 0, T) end -catch - (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂)#, VSU₃) -end -for V in spacelist - I = sectortype(first(V)) - Istr = TensorKit.type_repr(I) - println("---------------------------------------") - println("Tensors with symmetry: $Istr") - println("---------------------------------------") - @timedtestset "Tensors with symmetry: $Istr" verbose = true begin - V1, V2, V3, V4, V5 = V - @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) - @test @constinferred(hash(t)) == hash(deepcopy(t)) - @test scalartype(t) == T - @test norm(t) == 0 - @test codomain(t) == W - @test space(t) == (W ← one(W)) - @test domain(t) == one(W) - @test typeof(t) == @constinferred tensormaptype(spacetype(t), 5, 0, T) - end - end - @timedtestset "Tensor Dict conversion" begin + @testset "Tensor Dict conversion ($T)" for T in (Int, Float32, ComplexF64) + W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5 + t = TensorMap(rand, T, W) + d = convert(Dict, t) + @test t == convert(TensorMap, d) + end + + if hasfusiontensor(I) + @testset "Tensor Array conversion ($T)" for T in (Int, Float32, ComplexF64) W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5 - for T in (Int, Float32, ComplexF64) - t = TensorMap(rand, T, W) - d = convert(Dict, t) - @test t == convert(TensorMap, d) - end - end - if hasfusiontensor(I) - @timedtestset "Tensor Array conversion" begin - W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5 - for T in (Int, Float32, ComplexF64) - if T == Int - t = TensorMap(sz -> rand(-20:20, sz), W) - else - t = TensorMap(randn, T, W) - end - a = @constinferred convert(Array, t) - @test t ≈ @constinferred TensorMap(a, W) - end - end + t = TensorMap(T == Int ? sz -> rand(-20:20, sz) : randn, W) + a = @constinferred convert(Array, t) + @test t ≈ @constinferred TensorMap(a, W) end - @timedtestset "Basic linear algebra" begin + end + + @testset "Basic linear algebra ($T)" for T in (Float32, ComplexF64) + W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5 + t = TensorMap(rand, T, W) + @test scalartype(t) == T + @test space(t) == W + @test space(t') == W' + @test dim(t) == dim(space(t)) + @test codomain(t) == codomain(W) + @test domain(t) == domain(W) + @test isa(@constinferred(norm(t)), real(T)) + @test norm(t)^2 ≈ dot(t, t) + α = rand(T) + @test norm(α * t) ≈ abs(α) * norm(t) + @test norm(t + t, 2) ≈ 2 * norm(t, 2) + @test norm(t + t, 1) ≈ 2 * norm(t, 1) + @test norm(t + t, Inf) ≈ 2 * norm(t, Inf) + p = 3 * rand(Float64) + @test norm(t + t, p) ≈ 2 * norm(t, p) + @test norm(t) ≈ norm(t') + + t2 = TensorMap(rand, T, W) + β = rand(T) + @test @constinferred(dot(β * t2, α * t)) ≈ conj(β) * α * conj(dot(t, t2)) + @test dot(t2, t) ≈ conj(dot(t, t2)) + @test dot(t2, t) ≈ conj(dot(t2', t')) + @test dot(t2, t) ≈ dot(t', t2') + + i1 = @constinferred(isomorphism(Matrix{T}, V1 ⊗ V2, V2 ⊗ V1)) + i2 = @constinferred(isomorphism(Matrix{T}, V2 ⊗ V1, V1 ⊗ V2)) + @test i1 * i2 == @constinferred(id(Matrix{T}, V1 ⊗ V2)) + @test i2 * i1 == @constinferred(id(Matrix{T}, V2 ⊗ V1)) + + w = @constinferred(isometry(Matrix{T}, V1 ⊗ (oneunit(V1) ⊕ oneunit(V1)), + V1)) + @test dim(w) == 2 * dim(V1 ← V1) + @test w' * w == id(Matrix{T}, V1) + @test w * w' == (w * w')^2 + end + + if hasfusiontensor(I) + @testset "Basic linear algebra: test via conversion ($T)" for T in + (Float32, ComplexF64) W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5 - for T in (Float32, ComplexF64) - t = TensorMap(rand, T, W) - @test scalartype(t) == T - @test space(t) == W - @test space(t') == W' - @test dim(t) == dim(space(t)) - @test codomain(t) == codomain(W) - @test domain(t) == domain(W) - @test isa(@constinferred(norm(t)), real(T)) - @test norm(t)^2 ≈ dot(t, t) - α = rand(T) - @test norm(α * t) ≈ abs(α) * norm(t) - @test norm(t + t, 2) ≈ 2 * norm(t, 2) - @test norm(t + t, 1) ≈ 2 * norm(t, 1) - @test norm(t + t, Inf) ≈ 2 * norm(t, Inf) - p = 3 * rand(Float64) - @test norm(t + t, p) ≈ 2 * norm(t, p) - @test norm(t) ≈ norm(t') - - t2 = TensorMap(rand, T, W) - β = rand(T) - @test @constinferred(dot(β * t2, α * t)) ≈ conj(β) * α * conj(dot(t, t2)) - @test dot(t2, t) ≈ conj(dot(t, t2)) - @test dot(t2, t) ≈ conj(dot(t2', t')) - @test dot(t2, t) ≈ dot(t', t2') - - i1 = @constinferred(isomorphism(Matrix{T}, V1 ⊗ V2, V2 ⊗ V1)) - i2 = @constinferred(isomorphism(Matrix{T}, V2 ⊗ V1, V1 ⊗ V2)) - @test i1 * i2 == @constinferred(id(Matrix{T}, V1 ⊗ V2)) - @test i2 * i1 == @constinferred(id(Matrix{T}, V2 ⊗ V1)) - - w = @constinferred(isometry(Matrix{T}, V1 ⊗ (oneunit(V1) ⊕ oneunit(V1)), - V1)) - @test dim(w) == 2 * dim(V1 ← V1) - @test w' * w == id(Matrix{T}, V1) - @test w * w' == (w * w')^2 - end - end - if hasfusiontensor(I) - @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) - @test norm(t, 2) ≈ norm(convert(Array, t), 2) - @test dot(t2, t) ≈ dot(convert(Array, t2), convert(Array, t)) - α = rand(T) - @test convert(Array, α * t) ≈ α * convert(Array, t) - @test convert(Array, t + t) ≈ 2 * convert(Array, t) - end - end - @timedtestset "Real and imaginary parts" begin - W = V1 ⊗ V2 - for T in (Float64, ComplexF64, ComplexF32) - t = TensorMap(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 - end + t = TensorMap(rand, T, W) + t2 = TensorMap(rand, T, W) + @test norm(t, 2) ≈ norm(convert(Array, t), 2) + @test dot(t2, t) ≈ dot(convert(Array, t2), convert(Array, t)) + α = rand(T) + @test convert(Array, α * t) ≈ α * convert(Array, t) + @test convert(Array, t + t) ≈ 2 * convert(Array, t) end - @timedtestset "Tensor conversion" begin + @timedtestset "Real and imaginary parts" for T in (Float64, ComplexF64, ComplexF32) W = V1 ⊗ V2 - t = TensorMap(randn, Float64, W, W) - @test typeof(convert(TensorMap, t')) == typeof(t) - tc = complex(t) - @test convert(typeof(tc), t) == tc - @test typeof(convert(typeof(tc), t)) == typeof(tc) - @test typeof(convert(typeof(tc), t')) == typeof(tc) - @test Base.promote_typeof(t, tc) == typeof(tc) - @test Base.promote_typeof(tc, t) == typeof(tc + t) + t = TensorMap(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 + end + + @timedtestset "Tensor conversion" begin + W = V1 ⊗ V2 + t = TensorMap(randn, Float64, W, W) + @test typeof(convert(TensorMap, t')) == typeof(t) + tc = complex(t) + @test convert(typeof(tc), t) == tc + @test typeof(convert(typeof(tc), t)) == typeof(tc) + @test typeof(convert(typeof(tc), t')) == typeof(tc) + @test Base.promote_typeof(t, tc) == typeof(tc) + @test Base.promote_typeof(tc, t) == typeof(tc + t) + 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) + for k in 0:5 + for p in permutations(1:5) + p1 = ntuple(n -> p[n], k) + p2 = ntuple(n -> p[k + n], 5 - k) + t2 = @constinferred permute(t, (p1, p2)) + @test norm(t2) ≈ norm(t) + t2′ = permute(t′, (p1, p2)) + @test dot(t2′, t2) ≈ dot(t′, t) ≈ dot(transpose(t2′), transpose(t2)) + end end - @timedtestset "Permutations: test via inner product invariance" begin + end + + if hasfusiontensor(I) + @timedtestset "Permutations: test via conversion" begin W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5 t = Tensor(rand, ComplexF64, W) - t′ = Tensor(rand, ComplexF64, W) for k in 0:5 for p in permutations(1:5) p1 = ntuple(n -> p[n], k) p2 = ntuple(n -> p[k + n], 5 - k) - t2 = @constinferred permute(t, (p1, p2)) - @test norm(t2) ≈ norm(t) - t2′ = permute(t′, (p1, p2)) - @test dot(t2′, t2) ≈ dot(t′, t) ≈ dot(transpose(t2′), transpose(t2)) + t2 = permute(t, (p1, p2)) + a2 = convert(Array, t2) + @test a2 ≈ permutedims(convert(Array, t), (p1..., p2...)) + @test convert(Array, transpose(t2)) ≈ + permutedims(a2, (5, 4, 3, 2, 1)) end end end - if hasfusiontensor(I) - @timedtestset "Permutations: test via conversion" begin - W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5 - t = Tensor(rand, ComplexF64, W) - for k in 0:5 - for p in permutations(1:5) - p1 = ntuple(n -> p[n], k) - p2 = ntuple(n -> p[k + n], 5 - k) - t2 = permute(t, (p1, p2)) - a2 = convert(Array, t2) - @test a2 ≈ permutedims(convert(Array, t), (p1..., p2...)) - @test convert(Array, transpose(t2)) ≈ - permutedims(a2, (5, 4, 3, 2, 1)) - end - end - end + end + + @timedtestset "Full trace: test self-consistency" begin + t = Tensor(rand, ComplexF64, V1 ⊗ V2' ⊗ V2 ⊗ V1') + t2 = permute(t, ((1, 2), (4, 3))) + s = @constinferred tr(t2) + @test conj(s) ≈ tr(t2') + if !isdual(V1) + t2 = twist!(t2, 1) end - @timedtestset "Full trace: test self-consistency" begin - t = Tensor(rand, ComplexF64, V1 ⊗ V2' ⊗ V2 ⊗ V1') - t2 = permute(t, ((1, 2), (4, 3))) - s = @constinferred tr(t2) - @test conj(s) ≈ tr(t2') - if !isdual(V1) - t2 = twist!(t2, 1) - end - if isdual(V2) - t2 = twist!(t2, 2) - end - ss = tr(t2) - @tensor s2 = t[a, b, b, a] - @tensor t3[a, b] := t[a, c, c, b] - @tensor s3 = t3[a, a] - @test ss ≈ s2 - @test ss ≈ s3 + if isdual(V2) + t2 = twist!(t2, 2) end - @timedtestset "Partial trace: test self-consistency" begin + ss = tr(t2) + @tensor s2 = t[a, b, b, a] + @tensor t3[a, b] := t[a, c, c, b] + @tensor s3 = t3[a, a] + @test ss ≈ s2 + @test ss ≈ s3 + end + + @timedtestset "Partial trace: test self-consistency" begin + t = Tensor(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] + @test t2 ≈ t5 + end + + if hasfusiontensor(I) + @timedtestset "Trace: test via conversion" begin t = Tensor(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] - @test t2 ≈ t5 + @tensor t3[a, b] := convert(Array, t)[c, d, b, d, c, a] + @test t3 ≈ convert(Array, t2) end - if hasfusiontensor(I) - @timedtestset "Trace: test via conversion" begin - t = Tensor(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') - 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] - @test ta ≈ tb + end + + @timedtestset "Trace and contraction" begin + t1 = Tensor(rand, ComplexF64, V1 ⊗ V2 ⊗ V3) + t2 = Tensor(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] + @test ta ≈ tb + end + + if 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) + @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] + + @tensor HrA12array[a, s1, s2, c] := convert(Array, rhoL)[a, a'] * + conj(convert(Array, A1)[a', t1, b]) * + convert(Array, A2)[b, t2, c'] * + convert(Array, rhoR)[c', c] * + convert(Array, H)[s1, s2, t1, t2] + + @test HrA12array ≈ convert(Array, HrA12) end - if 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) - @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] - - @tensor HrA12array[a, s1, s2, c] := convert(Array, rhoL)[a, a'] * - conj(convert(Array, A1)[a', t1, b]) * - convert(Array, A2)[b, t2, c'] * - convert(Array, rhoR)[c', c] * - convert(Array, H)[s1, s2, t1, t2] - - @test HrA12array ≈ convert(Array, HrA12) - end + end + + @timedtestset "Multiplication and inverse: test compatibility" begin + 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) + @test t1 * (t1 \ t) ≈ t + @test (t / t2) * t2 ≈ t + @test t1 \ one(t1) ≈ inv(t1) + @test one(t1) / t1 ≈ pinv(t1) + @test_throws SpaceMismatch inv(t) + @test_throws SpaceMismatch t2 \ t + @test_throws SpaceMismatch t / t1 + tp = pinv(t) * t + @test tp ≈ tp * tp end - @timedtestset "Multiplication and inverse: test compatibility" begin + end + + if hasfusiontensor(I) + @timedtestset "Multiplication and inverse: test via conversion" begin W1 = V1 ⊗ V2 ⊗ V3 W2 = V4 ⊗ V5 - for T in (Float64, ComplexF64) + 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) - @test t1 * (t1 \ t) ≈ t - @test (t / t2) * t2 ≈ t - @test t1 \ one(t1) ≈ inv(t1) - @test one(t1) / t1 ≈ pinv(t1) - @test_throws SpaceMismatch inv(t) - @test_throws SpaceMismatch t2 \ t - @test_throws SpaceMismatch t / t1 - tp = pinv(t) * t - @test tp ≈ tp * tp - end - end - if hasfusiontensor(I) - @timedtestset "Multiplication and inverse: test via conversion" begin - 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) - d1 = dim(W1) - d2 = dim(W2) - At1 = reshape(convert(Array, t1), d1, d1) - At2 = reshape(convert(Array, t2), d2, d2) - At = reshape(convert(Array, t), d1, d2) - @test reshape(convert(Array, t1 * t), d1, d2) ≈ At1 * At - @test reshape(convert(Array, t1' * t), d1, d2) ≈ At1' * At - @test reshape(convert(Array, t2 * t'), d2, d1) ≈ At2 * At' - @test reshape(convert(Array, t2' * t'), d2, d1) ≈ At2' * At' - - @test reshape(convert(Array, inv(t1)), d1, d1) ≈ inv(At1) - @test reshape(convert(Array, pinv(t)), d2, d1) ≈ pinv(At) - - if T == Float32 || T == ComplexF32 - continue - end + d1 = dim(W1) + d2 = dim(W2) + At1 = reshape(convert(Array, t1), d1, d1) + At2 = reshape(convert(Array, t2), d2, d2) + At = reshape(convert(Array, t), d1, d2) + @test reshape(convert(Array, t1 * t), d1, d2) ≈ At1 * At + @test reshape(convert(Array, t1' * t), d1, d2) ≈ At1' * At + @test reshape(convert(Array, t2 * t'), d2, d1) ≈ At2 * At' + @test reshape(convert(Array, t2' * t'), d2, d1) ≈ At2' * At' - @test reshape(convert(Array, t1 \ t), d1, d2) ≈ At1 \ At - @test reshape(convert(Array, t1' \ t), d1, d2) ≈ At1' \ At - @test reshape(convert(Array, t2 \ t'), d2, d1) ≈ At2 \ At' - @test reshape(convert(Array, t2' \ t'), d2, d1) ≈ At2' \ At' + @test reshape(convert(Array, inv(t1)), d1, d1) ≈ inv(At1) + @test reshape(convert(Array, pinv(t)), d2, d1) ≈ pinv(At) - @test reshape(convert(Array, t2 / t), d2, d1) ≈ At2 / At - @test reshape(convert(Array, t2' / t), d2, d1) ≈ At2' / At - @test reshape(convert(Array, t1 / t'), d1, d2) ≈ At1 / At' - @test reshape(convert(Array, t1' / t'), d1, d2) ≈ At1' / At' + if T == Float32 || T == ComplexF32 + continue end + + @test reshape(convert(Array, t1 \ t), d1, d2) ≈ At1 \ At + @test reshape(convert(Array, t1' \ t), d1, d2) ≈ At1' \ At + @test reshape(convert(Array, t2 \ t'), d2, d1) ≈ At2 \ At' + @test reshape(convert(Array, t2' \ t'), d2, d1) ≈ At2' \ At' + + @test reshape(convert(Array, t2 / t), d2, d1) ≈ At2 / At + @test reshape(convert(Array, t2' / t), d2, d1) ≈ At2' / At + @test reshape(convert(Array, t1 / t'), d1, d2) ≈ At1 / At' + @test reshape(convert(Array, t1' / t'), d1, d2) ≈ At1' / At' end end - @timedtestset "Factorization" begin - 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)') - for t in ts - @testset "leftorth with $alg" for alg in - (TensorKit.QR(), TensorKit.QRpos(), - TensorKit.QL(), TensorKit.QLpos(), - TensorKit.Polar(), TensorKit.SVD(), - TensorKit.SDD()) - Q, R = @constinferred leftorth(t, ((3, 4, 2), (1, 5)); alg=alg) - QdQ = Q' * Q - @test QdQ ≈ one(QdQ) - @test Q * R ≈ permute(t, ((3, 4, 2), (1, 5))) - if alg isa Polar - @test isposdef(R) - @test domain(R) == codomain(R) == space(t, 1)' ⊗ space(t, 5)' - end - end - @testset "leftnull with $alg" for alg in - (TensorKit.QR(), TensorKit.SVD(), - TensorKit.SDD()) - N = @constinferred leftnull(t, ((3, 4, 2), (1, 5)); alg=alg) - NdN = N' * N - @test NdN ≈ one(NdN) - @test norm(N' * permute(t, ((3, 4, 2), (1, 5)))) < - 100 * eps(norm(t)) - end - @testset "rightorth with $alg" for alg in - (TensorKit.RQ(), TensorKit.RQpos(), - TensorKit.LQ(), TensorKit.LQpos(), - TensorKit.Polar(), TensorKit.SVD(), - TensorKit.SDD()) - L, Q = @constinferred rightorth(t, ((3, 4), (2, 1, 5)); alg=alg) - QQd = Q * Q' - @test QQd ≈ one(QQd) - @test L * Q ≈ permute(t, ((3, 4), (2, 1, 5))) - if alg isa Polar - @test isposdef(L) - @test domain(L) == codomain(L) == space(t, 3) ⊗ space(t, 4) - end - end - @testset "rightnull with $alg" for alg in - (TensorKit.LQ(), TensorKit.SVD(), - TensorKit.SDD()) - M = @constinferred rightnull(t, ((3, 4), (2, 1, 5)); alg=alg) - MMd = M * M' - @test MMd ≈ one(MMd) - @test norm(permute(t, ((3, 4), (2, 1, 5))) * M') < - 100 * eps(norm(t)) - end - @testset "tsvd with $alg" for alg in (TensorKit.SVD(), TensorKit.SDD()) - U, S, V = @constinferred tsvd(t, ((3, 4, 2), (1, 5)); alg=alg) - UdU = U' * U - @test UdU ≈ one(UdU) - VVd = V * V' - @test VVd ≈ one(VVd) - @test U * S * V ≈ permute(t, ((3, 4, 2), (1, 5))) + end + + @timedtestset "Factorization" begin + 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)') + for t in ts + @testset "leftorth with $alg" for alg in + (TensorKit.QR(), TensorKit.QRpos(), + TensorKit.QL(), TensorKit.QLpos(), + TensorKit.Polar(), TensorKit.SVD(), + TensorKit.SDD()) + Q, R = @constinferred leftorth(t, ((3, 4, 2), (1, 5)); alg=alg) + QdQ = Q' * Q + @test QdQ ≈ one(QdQ) + @test Q * R ≈ permute(t, ((3, 4, 2), (1, 5))) + if alg isa Polar + @test isposdef(R) + @test domain(R) == codomain(R) == space(t, 1)' ⊗ space(t, 5)' end end - @testset "empty tensor" begin - t = TensorMap(randn, T, V1 ⊗ V2, typeof(V1)()) - @testset "leftorth with $alg" for alg in - (TensorKit.QR(), TensorKit.QRpos(), - TensorKit.QL(), TensorKit.QLpos(), - TensorKit.Polar(), TensorKit.SVD(), - TensorKit.SDD()) - Q, R = @constinferred leftorth(t; alg=alg) - @test Q == t - @test dim(Q) == dim(R) == 0 - end - @testset "leftnull with $alg" for alg in - (TensorKit.QR(), TensorKit.SVD(), - TensorKit.SDD()) - N = @constinferred leftnull(t; alg=alg) - @test N' * N ≈ id(domain(N)) - @test N * N' ≈ id(codomain(N)) - end - @testset "rightorth with $alg" for alg in - (TensorKit.RQ(), TensorKit.RQpos(), - TensorKit.LQ(), TensorKit.LQpos(), - TensorKit.Polar(), TensorKit.SVD(), - TensorKit.SDD()) - L, Q = @constinferred rightorth(copy(t'); alg=alg) - @test Q == t' - @test dim(Q) == dim(L) == 0 - end - @testset "rightnull with $alg" for alg in - (TensorKit.LQ(), TensorKit.SVD(), - TensorKit.SDD()) - M = @constinferred rightnull(copy(t'); alg=alg) - @test M * M' ≈ id(codomain(M)) - @test M' * M ≈ id(domain(M)) - end - @testset "tsvd with $alg" for alg in (TensorKit.SVD(), TensorKit.SDD()) - U, S, V = @constinferred tsvd(t; alg=alg) - @test U == t - @test dim(U) == dim(S) == dim(V) + @testset "leftnull with $alg" for alg in + (TensorKit.QR(), TensorKit.SVD(), + TensorKit.SDD()) + N = @constinferred leftnull(t, ((3, 4, 2), (1, 5)); alg=alg) + NdN = N' * N + @test NdN ≈ one(NdN) + @test norm(N' * permute(t, ((3, 4, 2), (1, 5)))) < + 100 * eps(norm(t)) + end + @testset "rightorth with $alg" for alg in + (TensorKit.RQ(), TensorKit.RQpos(), + TensorKit.LQ(), TensorKit.LQpos(), + TensorKit.Polar(), TensorKit.SVD(), + TensorKit.SDD()) + L, Q = @constinferred rightorth(t, ((3, 4), (2, 1, 5)); alg=alg) + QQd = Q * Q' + @test QQd ≈ one(QQd) + @test L * Q ≈ permute(t, ((3, 4), (2, 1, 5))) + if alg isa Polar + @test isposdef(L) + @test domain(L) == codomain(L) == space(t, 3) ⊗ space(t, 4) end end - - t = Tensor(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))) - @test t2 * V ≈ V * D - - # Somehow moving these test before the previous one gives rise to errors - # with T=Float32 on x86 platforms. Is this an OpenBLAS issue? - VdV = V' * V; - VdV = (VdV + VdV') / 2 - @test isposdef(VdV) - - @test !isposdef(t2) # unlikely for non-hermitian map - t2 = (t2 + t2') - D, V = eigen(t2) - VdV = V' * V - @test VdV ≈ one(VdV) - D̃, Ṽ = @constinferred eigh(t2) - @test D ≈ D̃ - @test V ≈ Ṽ - λ = minimum(minimum(real(LinearAlgebra.diag(b))) - for (c, b) in blocks(D)) - @test isposdef(t2) == isposdef(λ) - @test isposdef(t2 - λ * one(t2) + 0.1 * one(t2)) - @test !isposdef(t2 - λ * one(t2) - 0.1 * one(t2)) + @testset "rightnull with $alg" for alg in + (TensorKit.LQ(), TensorKit.SVD(), + TensorKit.SDD()) + M = @constinferred rightnull(t, ((3, 4), (2, 1, 5)); alg=alg) + MMd = M * M' + @test MMd ≈ one(MMd) + @test norm(permute(t, ((3, 4), (2, 1, 5))) * M') < + 100 * eps(norm(t)) end - end - end - @timedtestset "Tensor truncation" begin - 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)') - for t in ts - U₀, S₀, V₀, = tsvd(t) - t = rmul!(t, 1 / norm(S₀, p)) - U, S, V, ϵ = @constinferred tsvd(t; trunc=truncerr(5e-1), p=p) - # @show p, ϵ - # @show domain(S) - # @test min(space(S,1), space(S₀,1)) != space(S₀,1) - U′, S′, V′, ϵ′ = tsvd(t; trunc=truncerr(nextfloat(ϵ)), p=p) - @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) - U′, S′, V′, ϵ′ = tsvd(t; trunc=truncdim(ceil(Int, dim(domain(S)))), - p=p) - @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) - U′, S′, V′, ϵ′ = tsvd(t; trunc=truncspace(space(S, 1)), p=p) - @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) - # results with truncationcutoff cannot be compared because they don't take degeneracy into account, and thus truncate differently - U, S, V, ϵ = tsvd(t; trunc=truncbelow(1 / dim(domain(S₀))), p=p) - # @show p, ϵ - # @show domain(S) - # @test min(space(S,1), space(S₀,1)) != space(S₀,1) - U′, S′, V′, ϵ′ = tsvd(t; trunc=truncspace(space(S, 1)), p=p) - @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) - end + @testset "tsvd with $alg" for alg in (TensorKit.SVD(), TensorKit.SDD()) + U, S, V = @constinferred tsvd(t, ((3, 4, 2), (1, 5)); alg=alg) + UdU = U' * U + @test UdU ≈ one(UdU) + VVd = V * V' + @test VVd ≈ one(VVd) + @test U * S * V ≈ permute(t, ((3, 4, 2), (1, 5))) end end - end - if hasfusiontensor(I) - @timedtestset "Tensor functions" begin - W = V1 ⊗ V2 - for T in (Float64, ComplexF64) - t = TensorMap(randn, T, W, W) - s = dim(W) - expt = @constinferred exp(t) - @test reshape(convert(Array, expt), (s, s)) ≈ - exp(reshape(convert(Array, t), (s, s))) - - @test (@constinferred sqrt(t))^2 ≈ t - @test reshape(convert(Array, sqrt(t^2)), (s, s)) ≈ - sqrt(reshape(convert(Array, t^2), (s, s))) - - @test exp(@constinferred log(expt)) ≈ expt - @test reshape(convert(Array, log(expt)), (s, s)) ≈ - log(reshape(convert(Array, expt), (s, s))) - - @test (@constinferred cos(t))^2 + (@constinferred sin(t))^2 ≈ id(W) - @test (@constinferred tan(t)) ≈ sin(t) / cos(t) - @test (@constinferred cot(t)) ≈ cos(t) / sin(t) - @test (@constinferred cosh(t))^2 - (@constinferred sinh(t))^2 ≈ id(W) - @test (@constinferred tanh(t)) ≈ sinh(t) / cosh(t) - @test (@constinferred coth(t)) ≈ cosh(t) / sinh(t) - - t1 = sin(t) - @test sin(@constinferred asin(t1)) ≈ t1 - t2 = cos(t) - @test cos(@constinferred acos(t2)) ≈ t2 - t3 = sinh(t) - @test sinh(@constinferred asinh(t3)) ≈ t3 - t4 = cosh(t) - @test cosh(@constinferred acosh(t4)) ≈ t4 - t5 = tan(t) - @test tan(@constinferred atan(t5)) ≈ t5 - t6 = cot(t) - @test cot(@constinferred acot(t6)) ≈ t6 - t7 = tanh(t) - @test tanh(@constinferred atanh(t7)) ≈ t7 - t8 = coth(t) - @test coth(@constinferred acoth(t8)) ≈ t8 + @testset "empty tensor" begin + t = TensorMap(randn, T, V1 ⊗ V2, typeof(V1)()) + @testset "leftorth with $alg" for alg in + (TensorKit.QR(), TensorKit.QRpos(), + TensorKit.QL(), TensorKit.QLpos(), + TensorKit.Polar(), TensorKit.SVD(), + TensorKit.SDD()) + Q, R = @constinferred leftorth(t; alg=alg) + @test Q == t + @test dim(Q) == dim(R) == 0 + end + @testset "leftnull with $alg" for alg in + (TensorKit.QR(), TensorKit.SVD(), + TensorKit.SDD()) + N = @constinferred leftnull(t; alg=alg) + @test N' * N ≈ id(domain(N)) + @test N * N' ≈ id(codomain(N)) + end + @testset "rightorth with $alg" for alg in + (TensorKit.RQ(), TensorKit.RQpos(), + TensorKit.LQ(), TensorKit.LQpos(), + TensorKit.Polar(), TensorKit.SVD(), + TensorKit.SDD()) + L, Q = @constinferred rightorth(copy(t'); alg=alg) + @test Q == t' + @test dim(Q) == dim(L) == 0 + end + @testset "rightnull with $alg" for alg in + (TensorKit.LQ(), TensorKit.SVD(), + TensorKit.SDD()) + M = @constinferred rightnull(copy(t'); alg=alg) + @test M * M' ≈ id(codomain(M)) + @test M' * M ≈ id(domain(M)) end + @testset "tsvd with $alg" for alg in (TensorKit.SVD(), TensorKit.SDD()) + U, S, V = @constinferred tsvd(t; alg=alg) + @test U == t + @test dim(U) == dim(S) == dim(V) + end + end + + t = Tensor(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))) + @test t2 * V ≈ V * D + @test !isposdef(t2) # unlikely for non-hermitian map + t2 = (t2 + t2') + D, V = eigen(t2) + VdV = V' * V + @test VdV ≈ one(VdV) + D̃, Ṽ = @constinferred eigh(t2) + @test D ≈ D̃ + @test V ≈ Ṽ + λ = minimum(minimum(real(LinearAlgebra.diag(b))) + for (c, b) in blocks(D)) + @test isposdef(t2) == isposdef(λ) + @test isposdef(t2 - λ * one(t2) + 0.1 * one(t2)) + @test !isposdef(t2 - λ * one(t2) - 0.1 * one(t2)) end 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 = 3 // 2 * leftorth(tA; alg=Polar())[1] - tB = 1 // 5 * leftorth(tB; alg=Polar())[1] - tC = TensorMap(rand, T, V1 ⊗ V3, V2 ⊗ V4) - t = @constinferred sylvester(tA, tB, tC) - @test codomain(t) == V1 ⊗ V3 - @test domain(t) == V2 ⊗ V4 - @test norm(tA * t + t * tB + tC) < - (norm(tA) + norm(tB) + norm(tC)) * eps(real(T))^(2 / 3) - if hasfusiontensor(I) - matrix(x) = reshape(convert(Array, x), dim(codomain(x)), dim(domain(x))) - @test matrix(t) ≈ sylvester(matrix(tA), matrix(tB), matrix(tC)) + end + @timedtestset "Tensor truncation" begin + 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)') + for t in ts + U₀, S₀, V₀, = tsvd(t) + t = rmul!(t, 1 / norm(S₀, p)) + U, S, V, ϵ = @constinferred tsvd(t; trunc=truncerr(5e-1), p=p) + # @show p, ϵ + # @show domain(S) + # @test min(space(S,1), space(S₀,1)) != space(S₀,1) + U′, S′, V′, ϵ′ = tsvd(t; trunc=truncerr(nextfloat(ϵ)), p=p) + @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) + U′, S′, V′, ϵ′ = tsvd(t; trunc=truncdim(ceil(Int, dim(domain(S)))), + p=p) + @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) + U′, S′, V′, ϵ′ = tsvd(t; trunc=truncspace(space(S, 1)), p=p) + @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) + # results with truncationcutoff cannot be compared because they don't take degeneracy into account, and thus truncate differently + U, S, V, ϵ = tsvd(t; trunc=truncbelow(1 / dim(domain(S₀))), p=p) + # @show p, ϵ + # @show domain(S) + # @test min(space(S,1), space(S₀,1)) != space(S₀,1) + U′, S′, V′, ϵ′ = tsvd(t; trunc=truncspace(space(S, 1)), p=p) + @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) end end 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) - t = @constinferred (t1 ⊗ t2) - @test norm(t) ≈ norm(t1) * norm(t2) + end + if hasfusiontensor(I) + @timedtestset "Tensor functions" begin + W = V1 ⊗ V2 + for T in (Float64, ComplexF64) + t = TensorMap(randn, T, W, W) + s = dim(W) + expt = @constinferred exp(t) + @test reshape(convert(Array, expt), (s, s)) ≈ + exp(reshape(convert(Array, t), (s, s))) + + @test (@constinferred sqrt(t))^2 ≈ t + @test reshape(convert(Array, sqrt(t^2)), (s, s)) ≈ + sqrt(reshape(convert(Array, t^2), (s, s))) + + @test exp(@constinferred log(expt)) ≈ expt + @test reshape(convert(Array, log(expt)), (s, s)) ≈ + log(reshape(convert(Array, expt), (s, s))) + + @test (@constinferred cos(t))^2 + (@constinferred sin(t))^2 ≈ id(W) + @test (@constinferred tan(t)) ≈ sin(t) / cos(t) + @test (@constinferred cot(t)) ≈ cos(t) / sin(t) + @test (@constinferred cosh(t))^2 - (@constinferred sinh(t))^2 ≈ id(W) + @test (@constinferred tanh(t)) ≈ sinh(t) / cosh(t) + @test (@constinferred coth(t)) ≈ cosh(t) / sinh(t) + + t1 = sin(t) + @test sin(@constinferred asin(t1)) ≈ t1 + t2 = cos(t) + @test cos(@constinferred acos(t2)) ≈ t2 + t3 = sinh(t) + @test sinh(@constinferred asinh(t3)) ≈ t3 + t4 = cosh(t) + @test cosh(@constinferred acosh(t4)) ≈ t4 + t5 = tan(t) + @test tan(@constinferred atan(t5)) ≈ t5 + t6 = cot(t) + @test cot(@constinferred acot(t6)) ≈ t6 + t7 = tanh(t) + @test tanh(@constinferred atanh(t7)) ≈ t7 + t8 = coth(t) + @test coth(@constinferred acoth(t8)) ≈ t8 end end - if 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) - t = @constinferred (t1 ⊗ t2) - d1 = dim(codomain(t1)) - d2 = dim(codomain(t2)) - d3 = dim(domain(t1)) - d4 = dim(domain(t2)) - At = convert(Array, t) - @test reshape(At, (d1, d2, d3, d4)) ≈ - reshape(convert(Array, t1), (d1, 1, d3, 1)) .* - reshape(convert(Array, t2), (1, d2, 1, d4)) - end + 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 = 3 // 2 * leftorth(tA; alg=Polar())[1] + tB = 1 // 5 * leftorth(tB; alg=Polar())[1] + tC = TensorMap(rand, T, V1 ⊗ V3, V2 ⊗ V4) + t = @constinferred sylvester(tA, tB, tC) + @test codomain(t) == V1 ⊗ V3 + @test domain(t) == V2 ⊗ V4 + @test norm(tA * t + t * tB + tC) < + (norm(tA) + norm(tB) + norm(tC)) * eps(real(T))^(2 / 3) + if hasfusiontensor(I) + matrix(x) = reshape(convert(Array, x), dim(codomain(x)), dim(domain(x))) + @test matrix(t) ≈ sylvester(matrix(tA), matrix(tB), matrix(tC)) end end - @timedtestset "Tensor product: test via tensor contraction" begin + 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) + t = @constinferred (t1 ⊗ t2) + @test norm(t) ≈ norm(t1) * norm(t2) + end + end + + if hasfusiontensor(I) + @timedtestset "Tensor product: test via conversion" begin for T in (Float32, ComplexF64) - t1 = Tensor(rand, T, V2 ⊗ V3 ⊗ V1) - t2 = Tensor(rand, T, V2 ⊗ V1 ⊗ V3) + t1 = TensorMap(rand, T, V2 ⊗ V3 ⊗ V1, V1) + t2 = TensorMap(rand, T, V2 ⊗ V1 ⊗ V3, V2) t = @constinferred (t1 ⊗ t2) - @tensor t′[1, 2, 3, 4, 5, 6] := t1[1, 2, 3] * t2[4, 5, 6] - @test t ≈ t′ + d1 = dim(codomain(t1)) + d2 = dim(codomain(t2)) + d3 = dim(domain(t1)) + d4 = dim(domain(t2)) + At = convert(Array, t) + @test reshape(At, (d1, d2, d3, d4)) ≈ + reshape(convert(Array, t1), (d1, 1, d3, 1)) .* + reshape(convert(Array, t2), (1, d2, 1, d4)) end end 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) + t = @constinferred (t1 ⊗ t2) + @tensor t′[1, 2, 3, 4, 5, 6] := t1[1, 2, 3] * t2[4, 5, 6] + @test t ≈ t′ + end + end + + println("Finished tests for $Istr.") end @timedtestset "Deligne tensor product: test via conversion" begin + Vtr = smallspace(Trivial) + VSU₂ = smallspace(SU2Irrep) + Vℤ₂ = smallspace(Z2Irrep) @testset for Vlist1 in (Vtr, VSU₂), Vlist2 in (Vtr, Vℤ₂) V1, V2, V3, V4, V5 = Vlist1 W1, W2, W3, W4, W5 = Vlist2 diff --git a/test/utility.jl b/test/utility.jl new file mode 100644 index 00000000..005b9c9b --- /dev/null +++ b/test/utility.jl @@ -0,0 +1,148 @@ +export smallset, randsector, hasfusiontensor, NewSU2Irrep +export smallspace, sectorlist, fusiontensor, force_planar +export LinearAlgebra + +using Base.Iterators: take, product +using TensorKit +using TensorKit: ProductSector, Trivial, fusiontensor +using TensorKit: ℙ, PlanarTrivial +using Random +using LinearAlgebra: LinearAlgebra +using Combinatorics +using CategoryData + +include("newsectors.jl") +using .NewSectors + +# TODO: add some sector with multiplicity +# sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewSU2Irrep, +# FibonacciAnyon, IsingAnyon, FermionParity, FermionNumber, FermionSpin, +# Z3Irrep ⊠ Z4Irrep, FermionParity ⊠ U1Irrep ⊠ SU2Irrep, +# FermionParity ⊠ SU2Irrep ⊠ SU2Irrep, NewSU2Irrep ⊠ NewSU2Irrep, +# Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon, Object{E6}) + +const spacelist = Dict(Trivial => (ℂ^3, (ℂ^4)', ℂ^5, ℂ^6, (ℂ^7)'), + Z2Irrep => (ℂ[Z2Irrep](0 => 1, 1 => 1), ℂ[Z2Irrep](0 => 1, 1 => 2)', + ℂ[Z2Irrep](0 => 3, 1 => 2)', + ℂ[Z2Irrep](0 => 2, 1 => 3), ℂ[Z2Irrep](0 => 2, 1 => 5)), + FermionParity => (ℂ[FermionParity](0 => 1, 1 => 1), + ℂ[FermionParity](0 => 1, 1 => 2)', + ℂ[FermionParity](0 => 3, 1 => 2)', + ℂ[FermionParity](0 => 2, 1 => 3), + ℂ[FermionParity](0 => 2, 1 => 5)), + Z3Irrep => (ℂ[Z3Irrep](0 => 1, 1 => 2, 2 => 2), + ℂ[Z3Irrep](0 => 3, 1 => 1, 2 => 1), + ℂ[Z3Irrep](0 => 2, 1 => 2, 2 => 1)', + ℂ[Z3Irrep](0 => 1, 1 => 2, 2 => 3), + ℂ[Z3Irrep](0 => 1, 1 => 3, 2 => 3)'), + U1Irrep => (ℂ[U1Irrep](0 => 1, 1 => 2, -1 => 2), + ℂ[U1Irrep](0 => 3, 1 => 1, -1 => 1), + ℂ[U1Irrep](0 => 2, 1 => 2, -1 => 1)', + ℂ[U1Irrep](0 => 1, 1 => 2, -1 => 3), + ℂ[U1Irrep](0 => 1, 1 => 3, -1 => 3)'), + FermionNumber => (ℂ[FermionNumber](0 => 1, 1 => 2, -1 => 2), + ℂ[FermionNumber](0 => 3, 1 => 1, -1 => 1), + ℂ[FermionNumber](0 => 2, 1 => 2, -1 => 1)', + ℂ[FermionNumber](0 => 1, 1 => 2, -1 => 3), + ℂ[FermionNumber](0 => 1, 1 => 3, -1 => 3)'), + CU1Irrep => (ℂ[CU1Irrep]((0, 0) => 1, (0, 1) => 2, 1 => 1), + ℂ[CU1Irrep]((0, 0) => 3, (0, 1) => 0, 1 => 1), + ℂ[CU1Irrep]((0, 0) => 1, (0, 1) => 0, 1 => 2)', + ℂ[CU1Irrep]((0, 0) => 2, (0, 1) => 2, 1 => 1), + ℂ[CU1Irrep]((0, 0) => 2, (0, 1) => 1, 1 => 2)'), + SU2Irrep => (ℂ[SU2Irrep](0 => 3, 1 // 2 => 1), + ℂ[SU2Irrep](0 => 2, 1 => 1), + ℂ[SU2Irrep](1 // 2 => 1, 1 => 1)', + ℂ[SU2Irrep](0 => 2, 1 // 2 => 2), + ℂ[SU2Irrep](0 => 1, 1 // 2 => 1, 3 // 2 => 1)'), + FermionSpin => (ℂ[FermionSpin](0 => 3, 1 // 2 => 1), + ℂ[FermionSpin](0 => 2, 1 => 1), + ℂ[FermionSpin](1 // 2 => 1, 1 => 1)', + ℂ[FermionSpin](0 => 2, 1 // 2 => 2), + ℂ[FermionSpin](0 => 1, 1 // 2 => 1, 3 // 2 => 1)'), + Object{E6} => (ℂ[Object{E6}](1 => 1, 3 => 1), + ℂ[Object{E6}](2 => 1, 3 => 1)', + ℂ[Object{E6}](1 => 1, 3 => 1), + ℂ[Object{E6}](1 => 1, 3 => 1), + ℂ[Object{E6}](1 => 1, 2 => 1, 3 => 1)'), + FibonacciAnyon => (GradedSpace{FibonacciAnyon}(:I => 2, :τ => 2), + GradedSpace{FibonacciAnyon}(:I => 1, :τ => 2)', + GradedSpace{FibonacciAnyon}(:I => 2, :τ => 2), + GradedSpace{FibonacciAnyon}(:I => 2, :τ => 1), + GradedSpace{FibonacciAnyon}(:I => 2, :τ => 2)'), + IsingAnyon => (GradedSpace{IsingAnyon}(:I => 2, :psi => 1, + :sigma => 1), + GradedSpace{IsingAnyon}(:I => 2, + :sigma => 1)', + GradedSpace{IsingAnyon}(:I => 2, :psi => 1, + :sigma => 1), + GradedSpace{IsingAnyon}(:I => 2, :psi => 1), + GradedSpace{IsingAnyon}(:I => 2, :psi => 1, + :sigma => 1)')) + +function smallspace(::Type{I}) where {I<:Sector} + return get(spacelist, I, nothing) +end + +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 + +""" + force_planar(obj) + +Replace an object with a planar equivalent -- i.e. one that disallows braiding. +""" +force_planar(V::ComplexSpace) = isdual(V) ? (ℙ^dim(V))' : ℙ^dim(V) +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))) + copyto!(blocks(tdst)[PlanarTrivial()], blocks(tsrc)[Trivial()]) + return tdst +end +function force_planar(tsrc::TensorMap{<:GradedSpace}) + if BraidingStyle(sectortype(tsrc)) isa TensorKit.SymmetricBraiding + return tsrc + end + tdst = TensorMap(undef, scalartype(tsrc), + force_planar(codomain(tsrc)) ← force_planar(domain(tsrc))) + for (c, b) in blocks(tsrc) + copyto!(blocks(tdst)[c ⊠ PlanarTrivial()], b) + end + return tdst +end