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/fusiontrees.jl b/test/fusiontrees.jl index 61dc7d38..34840f0c 100644 --- a/test/fusiontrees.jl +++ b/test/fusiontrees.jl @@ -1,10 +1,14 @@ println("------------------------------------") println("Fusion Trees") println("------------------------------------") -ti = time() -@timedtestset "Fusion trees for $(TensorKit.type_repr(I))" verbose = true for I in - sectorlist + +module FusionTreeTests + +include("utility.jl") + +@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 +22,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,16 +40,16 @@ 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 @@ -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,17 +175,17 @@ ti = time() end end end - @testset "Fusion tree $Istr: elementy artin braid" begin + @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 TK.artin_braid(f, i) + d1 = @constinferred TensorKit.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) + for (f2, coeff2) in TensorKit.artin_braid(f1, i; inv=true) d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 end end @@ -196,24 +201,24 @@ ti = time() end f = rand(collect(it)) - d1 = TK.artin_braid(f, 2) + d1 = TensorKit.artin_braid(f, 2) d2 = empty(d1) for (f1, coeff1) in d1 - for (f2, coeff2) in TK.artin_braid(f1, 3) + 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 TK.artin_braid(f1, 3; inv=true) + 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 TK.artin_braid(f1, 2; inv=true) + for (f2, coeff2) in TensorKit.artin_braid(f1, 2; inv=true) d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 end end @@ -226,7 +231,7 @@ ti = time() end end end - @testset "Fusion tree $Istr: braiding and permuting" begin + @testset "braiding and permuting" begin f = rand(collect(fusiontrees(out, in, isdual))) p = tuple(randperm(N)...) ip = invperm(p) @@ -259,7 +264,7 @@ ti = time() end end - @testset "Fusion tree $Istr: merging" begin + @testset "merging" begin N = 3 out1 = ntuple(n -> randsector(I), N) in1 = rand(collect(⊗(out1...))) @@ -268,27 +273,27 @@ 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 (f, coeff) in TensorKit.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, μ′) + 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, ν′) + for (t, coeff) in TensorKit.merge(f2, f1, c, ν′) trees2[t] = get(trees2, t, zero(valtype(trees2))) + coeff * R[μ, ν] end end @@ -340,14 +345,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 +390,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 +446,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 +511,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 +532,8 @@ ti = time() end end end + + println("Finished tests for $Istr.") +end + 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..a623ffce 100644 --- a/test/planar.jl +++ b/test/planar.jl @@ -1,6 +1,15 @@ -using TensorKit, TensorOperations, Test +println("------------------------------------") +println("Planar") +println("------------------------------------") + +module PlanarTests + +using Test, TestExtras +using TensorKit using TensorKit: planaradd!, planartrace!, planarcontract! using TensorKit: PlanarTrivial, ℙ +using TensorOperations +import TensorKit: BraidingTensor """ force_planar(obj) @@ -27,125 +36,317 @@ function force_planar(tsrc::TensorMap{<:GradedSpace}) return tdst end -@testset "planar methods" verbose = true begin - @testset "planaradd" begin - A = TensorMap(randn, ℂ^2 ⊗ ℂ^3 ← ℂ^6 ⊗ ℂ^5 ⊗ ℂ^4) - C = TensorMap(randn, (ℂ^5)' ⊗ (ℂ^6)' ← ℂ^4 ⊗ (ℂ^3)' ⊗ (ℂ^2)') - A′ = force_planar(A) - C′ = force_planar(C) - p = ((4, 3), (5, 2, 1)) +# @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)) - @test force_planar(tensoradd!(C, p, A, :N, true, true)) ≈ - planaradd!(C′, A′, p, true, true) - end +# @test force_planar(tensoradd!(C, p, A, :N, true, true)) ≈ +# planaradd!(C′, A′, p, true, true) +# end - @testset "planartrace" begin - A = TensorMap(randn, ℂ^2 ⊗ ℂ^3 ← ℂ^2 ⊗ ℂ^5 ⊗ ℂ^4) - C = TensorMap(randn, (ℂ^5)' ⊗ ℂ^3 ← ℂ^4) - A′ = force_planar(A) - C′ = force_planar(C) - p = ((4, 2), (5,)) - q = ((1,), (3,)) +# @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,)) - @test force_planar(tensortrace!(C, p, A, q, :N, true, true)) ≈ - planartrace!(C′, A′, p, q, true, true) - end +# @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, ℂ^2 ⊗ ℂ^3 ← ℂ^2 ⊗ ℂ^5 ⊗ ℂ^4) - B = TensorMap(randn, ℂ^2 ⊗ ℂ^4 ← ℂ^4 ⊗ ℂ^3) - C = TensorMap(randn, (ℂ^5)' ⊗ (ℂ^2)' ⊗ ℂ^2 ← (ℂ^2)' ⊗ ℂ^4) +# @testset "planarcontract" begin +# A = TensorMap(randn, ℂ^2 ⊗ ℂ^3 ← ℂ^2 ⊗ ℂ^5 ⊗ ℂ^4) +# B = TensorMap(randn, ℂ^2 ⊗ ℂ^4 ← ℂ^4 ⊗ ℂ^3) +# C = TensorMap(randn, (ℂ^5)' ⊗ (ℂ^2)' ⊗ ℂ^2 ← (ℂ^2)' ⊗ ℂ^4) - A′ = force_planar(A) - B′ = force_planar(B) - C′ = force_planar(C) +# 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)) +# 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 +# @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 "@planar" verbose = true begin - T = ComplexF64 - @testset "MPS networks" begin - P = ℂ^2 - Vmps = ℂ^12 - Vmpo = ℂ^4 - - # ∂AC - # ------- - x = TensorMap(randn, T, Vmps ⊗ P ← Vmps) - O = TensorMap(randn, T, Vmpo ⊗ P ← P ⊗ Vmpo) - GL = TensorMap(randn, T, Vmps ⊗ Vmpo' ← Vmps) - GR = TensorMap(randn, T, Vmps ⊗ Vmpo ← Vmps) - - x′ = force_planar(x) - O′ = force_planar(O) - GL′ = force_planar(GL) - GR′ = force_planar(GR) - - @tensor y[-1 -2; -3] := GL[-1 2; 1] * x[1 3; 4] * O[2 -2; 3 5] * GR[4 5; -3] - @planar y′[-1 -2; -3] := GL′[-1 2; 1] * x′[1 3; 4] * O′[2 -2; 3 5] * GR′[4 5; -3] - @test force_planar(y) ≈ y′ - - # ∂AC2 - # ------- - x2 = TensorMap(randn, T, Vmps ⊗ P ← Vmps ⊗ P') - x2′ = force_planar(x2) - @tensor contractcheck = true y2[-1 -2; -3 -4] := GL[-1 7; 6] * x2[6 5; 1 3] * - O[7 -2; 5 4] * O[4 -4; 3 2] * - GR[1 2; -3] - @planar y2′[-1 -2; -3 -4] := GL′[-1 7; 6] * x2′[6 5; 1 3] * O′[7 -2; 5 4] * - O′[4 -4; 3 2] * GR′[1 2; -3] - @test force_planar(y2) ≈ y2′ - - # transfer matrix - # ---------------- - v = TensorMap(randn, T, Vmps ← Vmps) - v′ = force_planar(v) - @tensor ρ[-1; -2] := x[-1 2; 1] * conj(x[-2 2; 3]) * v[1; 3] - @planar ρ′[-1; -2] := x′[-1 2; 1] * conj(x′[-2 2; 3]) * v′[1; 3] - @test force_planar(ρ) ≈ ρ′ - - @tensor ρ2[-1 -2; -3] := GL[1 -2; 3] * x[3 2; -3] * conj(x[1 2; -1]) - @plansor ρ3[-1 -2; -3] := GL[1 2; 4] * x[4 5; -3] * τ[2 3; 5 -2] * conj(x[1 3; -1]) - @planar ρ2′[-1 -2; -3] := GL′[1 2; 4] * x′[4 5; -3] * τ[2 3; 5 -2] * - conj(x′[1 3; -1]) - @test force_planar(ρ2) ≈ ρ2′ - @test ρ2 ≈ ρ3 - end +# @testset "@planar" verbose = true begin +# T = ComplexF64 +# @testset "MPS networks" begin +# P = ℂ^2 +# Vmps = ℂ^12 +# Vmpo = ℂ^4 - @testset "MERA networks" begin - Vmera = ℂ^2 - - u = TensorMap(randn, T, Vmera ⊗ Vmera ← Vmera ⊗ Vmera) - w = TensorMap(randn, T, Vmera ⊗ Vmera ← Vmera) - ρ = TensorMap(randn, T, Vmera ⊗ Vmera ⊗ Vmera ← Vmera ⊗ Vmera ⊗ Vmera) - h = TensorMap(randn, T, Vmera ⊗ Vmera ⊗ Vmera ← Vmera ⊗ Vmera ⊗ Vmera) - - u′ = force_planar(u) - w′ = force_planar(w) - ρ′ = force_planar(ρ) - h′ = force_planar(h) - - @tensor begin - C = (((((((h[9 3 4; 5 1 2] * u[1 2; 7 12]) * conj(u[3 4; 11 13])) * - (u[8 5; 15 6] * w[6 7; 19])) * - (conj(u[8 9; 17 10]) * conj(w[10 11; 22]))) * - ((w[12 14; 20] * conj(w[13 14; 23])) * ρ[18 19 20; 21 22 23])) * - w[16 15; 18]) * conj(w[16 17; 21])) - end - @planar begin - C′ = (((((((h′[9 3 4; 5 1 2] * u′[1 2; 7 12]) * conj(u′[3 4; 11 13])) * - (u′[8 5; 15 6] * w′[6 7; 19])) * - (conj(u′[8 9; 17 10]) * conj(w′[10 11; 22]))) * - ((w′[12 14; 20] * conj(w′[13 14; 23])) * ρ′[18 19 20; 21 22 23])) * - w′[16 15; 18]) * conj(w′[16 17; 21])) +# # ∂AC +# # ------- +# x = TensorMap(randn, T, Vmps ⊗ P ← Vmps) +# O = TensorMap(randn, T, Vmpo ⊗ P ← P ⊗ Vmpo) +# GL = TensorMap(randn, T, Vmps ⊗ Vmpo' ← Vmps) +# GR = TensorMap(randn, T, Vmps ⊗ Vmpo ← Vmps) + +# x′ = force_planar(x) +# O′ = force_planar(O) +# GL′ = force_planar(GL) +# GR′ = force_planar(GR) + +# @tensor y[-1 -2; -3] := GL[-1 2; 1] * x[1 3; 4] * O[2 -2; 3 5] * GR[4 5; -3] +# @planar y′[-1 -2; -3] := GL′[-1 2; 1] * x′[1 3; 4] * O′[2 -2; 3 5] * GR′[4 5; -3] +# @test force_planar(y) ≈ y′ + +# # ∂AC2 +# # ------- +# x2 = TensorMap(randn, T, Vmps ⊗ P ← Vmps ⊗ P') +# x2′ = force_planar(x2) +# @tensor contractcheck = true y2[-1 -2; -3 -4] := GL[-1 7; 6] * x2[6 5; 1 3] * +# O[7 -2; 5 4] * O[4 -4; 3 2] * +# GR[1 2; -3] +# @planar y2′[-1 -2; -3 -4] := GL′[-1 7; 6] * x2′[6 5; 1 3] * O′[7 -2; 5 4] * +# O′[4 -4; 3 2] * GR′[1 2; -3] +# @test force_planar(y2) ≈ y2′ + +# # transfer matrix +# # ---------------- +# v = TensorMap(randn, T, Vmps ← Vmps) +# v′ = force_planar(v) +# @tensor ρ[-1; -2] := x[-1 2; 1] * conj(x[-2 2; 3]) * v[1; 3] +# @planar ρ′[-1; -2] := x′[-1 2; 1] * conj(x′[-2 2; 3]) * v′[1; 3] +# @test force_planar(ρ) ≈ ρ′ + +# @tensor ρ2[-1 -2; -3] := GL[1 -2; 3] * x[3 2; -3] * conj(x[1 2; -1]) +# @plansor ρ3[-1 -2; -3] := GL[1 2; 4] * x[4 5; -3] * τ[2 3; 5 -2] * conj(x[1 3; -1]) +# @planar ρ2′[-1 -2; -3] := GL′[1 2; 4] * x′[4 5; -3] * τ[2 3; 5 -2] * +# conj(x′[1 3; -1]) +# @test force_planar(ρ2) ≈ ρ2′ +# @test ρ2 ≈ ρ3 +# end + +# @testset "MERA networks" begin +# Vmera = ℂ^2 + +# u = TensorMap(randn, T, Vmera ⊗ Vmera ← Vmera ⊗ Vmera) +# w = TensorMap(randn, T, Vmera ⊗ Vmera ← Vmera) +# ρ = TensorMap(randn, T, Vmera ⊗ Vmera ⊗ Vmera ← Vmera ⊗ Vmera ⊗ Vmera) +# h = TensorMap(randn, T, Vmera ⊗ Vmera ⊗ Vmera ← Vmera ⊗ Vmera ⊗ Vmera) + +# u′ = force_planar(u) +# w′ = force_planar(w) +# ρ′ = force_planar(ρ) +# h′ = force_planar(h) + +# @tensor begin +# C = (((((((h[9 3 4; 5 1 2] * u[1 2; 7 12]) * conj(u[3 4; 11 13])) * +# (u[8 5; 15 6] * w[6 7; 19])) * +# (conj(u[8 9; 17 10]) * conj(w[10 11; 22]))) * +# ((w[12 14; 20] * conj(w[13 14; 23])) * ρ[18 19 20; 21 22 23])) * +# w[16 15; 18]) * conj(w[16 17; 21])) +# end +# @planar begin +# C′ = (((((((h′[9 3 4; 5 1 2] * u′[1 2; 7 12]) * conj(u′[3 4; 11 13])) * +# (u′[8 5; 15 6] * w′[6 7; 19])) * +# (conj(u′[8 9; 17 10]) * conj(w′[10 11; 22]))) * +# ((w′[12 14; 20] * conj(w′[13 14; 23])) * ρ′[18 19 20; 21 22 23])) * +# w′[16 15; 18]) * conj(w′[16 17; 21])) +# end +# @test C ≈ C′ +# end +# end + +spacelist = Dict(FermionParity => GradedSpace{FermionParity}(0 => 1, 1 => 1), + FermionSpin => GradedSpace{FermionSpin}(0 => 2, 1 / 2 => 2, 1 => 1, + 3 / 2 => 1), + FibonacciAnyon => GradedSpace{FibonacciAnyon}(:I => 2, :τ => 2), + IsingAnyon => GradedSpace{IsingAnyon}(:I => 2, :psi => 1, :sigma => 1), +) + +@testset "BraidingTensor ($(TensorKit.type_repr(I)))" for I in keys(spacelist) + V = spacelist[I] + + @testset "BraidingTensor conversion" begin + for (V1, V2) in ((V, V), (V', V), (V, V'), (V', V')) + τ = BraidingTensor(V1, V2) + + @test domain(τ) == V1 ⊗ V2 + @test codomain(τ) == V2 ⊗ V1 + + for (c, b) in blocks(copy(τ)) + @test b ≈ block(τ, c) + end + + @test domain(τ') == codomain(τ) + @test codomain(τ') == domain(τ) + + for (c, b) in blocks(copy(τ')) + @test b ≈ block(τ', c) + end end - @test C ≈ C′ end + + t = TensorMap(randn, V * V' * V' * V, V * V') + + ττ = copy(BraidingTensor(V, V')) + @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(V', V')) + @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(V', V)) + @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(V, V')) + @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(V, V)) + @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(V', V)) + @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(V', V)) + @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(V, V')) + @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(V, V)) + @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(V', V)) + @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(V, V)) + @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(V, V')) + @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(V, V')) + @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(V', V')) + @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(V', V)) + @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(V', V)) + @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(V, V')) + @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 + end diff --git a/test/runtests.jl b/test/runtests.jl index 4738f065..082c69a6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,69 +1,35 @@ 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 -include("newsectors.jl") -using .NewSectors +const GROUP = uppercase(get(ENV, "GROUP", "ALL")) -const TK = TensorKit +@testset verbose=true begin + if GROUP == "ALL" || GROUP == "SECTORS" + @testset "Sectors" verbose = true begin + include("sectors.jl") + end + end -Random.seed!(1234) + if GROUP == "ALL" || GROUP == "SPACES" + @testset "Spaces" verbose = true begin + include("spaces.jl") + end + 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) + if GROUP == "ALL" || GROUP == "FUSIONTREES" + @testset "Fusiontrees" verbose = true begin + include("fusiontrees.jl") + end 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) + + if GROUP == "ALL" || GROUP == "TENSORS" + @testset "Tensors" verbose = true begin + include("tensors.jl") + end + end + + if GROUP == "ALL" || GROUP == "PLANAR" + @testset "Planar" verbose = true begin + include("planar.jl") end 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") -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..a2e08ae0 100644 --- a/test/sectors.jl +++ b/test/sectors.jl @@ -1,118 +1,129 @@ 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...) - end - @testset "Sector $Istr: Value iterator" begin - @test eltype(values(I)) == I - sprev = one(I) - for (i, s) in enumerate(values(I)) - @test !isless(s, sprev) # confirm compatibility with sort order - if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector - @test_throws ArgumentError values(I)[i] - @test_throws ArgumentError TensorKit.findindex(values(I), s) - elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int}) - @test s == @constinferred (values(I)[i]) - @test TensorKit.findindex(values(I), s) == i - end - sprev = s - i >= 10 && break - end - @test one(I) == first(values(I)) + +module SectorTests + +include("utility.jl") + +@testset "$(TensorKit.type_repr(I))" verbose = true 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...) + @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 TensorKit.findindex(values(I), one(I)) + @test_throws ArgumentError values(I)[i] + @test_throws ArgumentError TensorKit.findindex(values(I), s) 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 + @test s == @constinferred (values(I)[i]) + @test TensorKit.findindex(values(I), s) == i end + sprev = s + i >= 10 && break end - if hasfusiontensor(I) - @testset "Sector $I: fusion tensor and F-move and R-move" begin - for a in smallset(I), b in smallset(I) - for c in ⊗(a, b) - X1 = permutedims(fusiontensor(a, b, c), (2, 1, 3, 4)) - X2 = fusiontensor(b, a, c) - l = dim(a) * dim(b) * dim(c) - R = LinearAlgebra.transpose(Rsymbol(a, b, c)) - sz = (l, convert(Int, Nsymbol(a, b, c))) - @test reshape(X1, sz) ≈ reshape(X2, sz) * R - end - end - for a in smallset(I), b in smallset(I), c in smallset(I) - for e in ⊗(a, b), f in ⊗(b, c) - for d in intersect(⊗(e, c), ⊗(a, f)) - X1 = fusiontensor(a, b, e) - X2 = fusiontensor(e, c, d) - Y1 = fusiontensor(b, c, f) - Y2 = fusiontensor(a, f, d) - @tensor f1[-1, -2, -3, -4] := conj(Y2[a, f, d, -4]) * - conj(Y1[b, c, f, -3]) * - X1[a, b, e, -1] * X2[e, c, d, -2] - if FusionStyle(I) isa MultiplicityFreeFusion - f2 = fill(Fsymbol(a, b, c, d, e, f) * dim(d), (1, 1, 1, 1)) - else - f2 = Fsymbol(a, b, c, d, e, f) * dim(d) - end - @test isapprox(f1, f2; atol=1e-12, rtol=1e-12) + @test one(I) == first(values(I)) + if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector + @test_throws ArgumentError TensorKit.findindex(values(I), one(I)) + elseif hasmethod(Base.getindex, Tuple{typeof(values(I)),Int}) + @test (@constinferred TensorKit.findindex(values(I), one(I))) == 1 + for s in smallset(I) + @test (@constinferred values(I)[TensorKit.findindex(values(I), s)]) == s + end + end + end + + @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 + @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12) end end - @testset "Sector $Istr: Unitarity of F-move" begin + 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 + @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 + + if hasfusiontensor(I) + @testset "Fusion tensor and 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 + 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 - F = hvcat(length(fs), Fblocks...) + @test isapprox(f1, f2; atol=1e-12, rtol=1e-12) 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) + + @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)) + X2 = fusiontensor(b, a, c) + l = dim(a) * dim(b) * dim(c) + R = LinearAlgebra.transpose(Rsymbol(a, b, c)) + sz = (l, convert(Int, Nsymbol(a, b, c))) + @test reshape(X1, sz) ≈ reshape(X2, sz) * R + end end end end end + +end diff --git a/test/spaces.jl b/test/spaces.jl index 8d02c97f..4003c590 100644 --- a/test/spaces.jl +++ b/test/spaces.jl @@ -1,405 +1,410 @@ 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}} ∈ ℂ) +module SpaceTests - @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} ∈ ℂ) +include("utility.jl") - @test T ⊆ ℝ - @test !(Complex{T} ⊆ ℝ) - @test T ⊆ ℂ - @test Complex{T} ⊆ ℂ - end - end +@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 !(ℂ ⊆ ℝ) - @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 - end + 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}} ∈ ℂ) - @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')) - @test_throws MethodError (⊕(ℝ^d, ℂ^d)) - @test_throws MethodError (⊗(ℝ^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) + @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: 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) + @test T ⊆ ℝ + @test !(Complex{T} ⊆ ℝ) + @test T ⊆ ℂ + @test Complex{T} ⊆ ℂ end +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 "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 + +@testset "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')) + @test_throws MethodError (⊕(ℝ^d, ℂ^d)) + @test_throws MethodError (⊗(ℝ^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 "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 + +@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 - @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 end diff --git a/test/tensors.jl b/test/tensors.jl index 910424d1..4597cd13 100644 --- a/test/tensors.jl +++ b/test/tensors.jl @@ -1,598 +1,562 @@ -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₃) - 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 +println("------------------------------------") +println("Tensors") +println("------------------------------------") + +module TensorTests -spacelist = try +include("utility.jl") + +sectortypes = try if ENV["CI"] == "true" println("Detected running on CI") if Sys.iswindows() - (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂) + (Trivial, Z2Irrep, FermionParity, Z3Irrep, U1Irrep, + FermionNumber, CU1Irrep, SU2Irrep) elseif Sys.isapple() - (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VfU₁, VfSU₂)#, VSU₃) + (Trivial, Z2Irrep, FermionParity, Z3Irrep, FermionNumber, FermionSpin) else - (Vtr, Vℤ₂, Vfℤ₂, VU₁, VCU₁, VSU₂, VfSU₂)#, VSU₃) + (Trivial, Z2Irrep, FermionParity, U1Irrep, CU1Irrep, SU2Irrep, FermionSpin) end else - (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂)#, VSU₃) + keys(spacelist) end catch - (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂)#, VSU₃) + keys(spacelist) +end + +for V in getindex.(Ref(spacelist), sectortypes) + 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 -for V in spacelist - I = sectortype(first(V)) +@testset "$(TensorKit.type_repr(I))" verbose = true for I in sectortypes 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 + println("Starting tests for $Istr...") + + V1, V2, V3, V4, V5 = spacelist[I] + + @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 + + @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 + t = TensorMap(T == Int ? sz -> rand(-20:20, sz) : randn, W) + a = @constinferred convert(Array, t) + @test t ≈ @constinferred TensorMap(a, W) 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 - 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 + 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 - 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 - 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 - @timedtestset "Permutations: test via inner product invariance" begin + 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 + 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 - 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 + @tensor t3[a, b] := convert(Array, t)[c, d, b, d, c, a] + @test t3 ≈ convert(Array, t2) 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))) - D̃, Ṽ = @constinferred eig(t, ((1, 3), (2, 4))) - @test D ≈ D̃ - @test V ≈ Ṽ - VdV = V' * V - VdV = (VdV + VdV') / 2 - @test isposdef(VdV) - 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)) + @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))) + D̃, Ṽ = @constinferred eig(t, ((1, 3), (2, 4))) + @test D ≈ D̃ + @test V ≈ Ṽ + VdV = V' * V + VdV = (VdV + VdV') / 2 + @test isposdef(VdV) + 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 @@ -614,3 +578,5 @@ end end end end + +end diff --git a/test/utility.jl b/test/utility.jl new file mode 100644 index 00000000..ce532728 --- /dev/null +++ b/test/utility.jl @@ -0,0 +1,95 @@ +using Base.Iterators: take, product +using TensorKit +using TensorKit: ProductSector, Trivial, fusiontensor +using TensorKit: pentagon_equation, hexagon_equation +using TensorOperations +using Test, TestExtras +using Random +using LinearAlgebra: LinearAlgebra +using Combinatorics + +Random.seed!(1234) + +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) + +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)')) + +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