Update sector tests

Project.toml

@@ -35,7 +35,7 @@ Serialization = "1"
 SparseArrayKit = "0.4"
 SparseArrays = "1"
 TensorKit = "0.16.4, 0.17"
-TensorKitSectors = "0.3.11"
+TensorKitSectors = "0.3.12"
 TensorOperations = "5"
 Test = "1"
 TestExtras = "0.3"

src/clebschgordan.jl

@@ -20,10 +20,10 @@ CGC(s1::I, s2::I, s3::I) where {I <: SUNIrrep} = CGC(sectorscalartype(I), s1, s2
 CGC(::Type{T}, s1::I, s2::I, s3::I) where {T, I <: SUNIrrep} = SUN_CGC(T, s1, s2, s3)
 
 @cached function SUN_CGC(::Type{T}, s1::I, s2::I, s3::I)::SparseArray{T, 4} where {T <: Real, I <: SUNIrrep}
-    if isone(s1)
+    if isunit(s1)
         @assert s2 == s3
         CGC = trivial_CGC(T, s2, true)
-    elseif isone(s2)
+    elseif isunit(s2)
         @assert s1 == s3
         CGC = trivial_CGC(T, s1, false)
     else

src/sector.jl

@@ -134,13 +134,14 @@ TensorKitSectors.fusionscalartype(::Type{<:SUNIrrep}) = Float64
 TensorKitSectors.Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I <: SUNIrrep} =
     SUN_Fsymbol(a, b, c, d, e, f)
 @cached function SUN_Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I)::Array{Float64, 4} where {I <: SUNIrrep}
-    N1 = Nsymbol(a, b, e)
-    N2 = Nsymbol(e, c, d)
-    N3 = Nsymbol(b, c, f)
-    N4 = Nsymbol(a, f, d)
-
-    (N1 == 0 || N2 == 0 || N3 == 0 || N4 == 0) &&
-        return fill(zero(fusionscalartype(I)), N1, N2, N3, N4)
+    return TensorKitSectors.Fsymbol_from_fusiontensor(a, b, c, d, e, f)
+end
+MemoizationKit.default_cache_settings(::typeof(SUN_Fsymbol)) = (; maxsize = 10^5)
+function TensorKitSectors.Fsymbol_from_fusiontensor(
+        a::I, b::I, c::I, d::I, e::I, f::I
+    ) where {I <: SUNIrrep}
+    Nabe, Necd, Nbcd, Nafd = Nsymbol(a, b, e), Nsymbol(e, c, d), Nsymbol(b, c, f), Nsymbol(a, f, d)
+    iszero(Nabe * Necd * Nbcd * Nafd) && return zeros(fusionscalartype(I), Nabe, Necd, Nbcd, Nafd)
 
     # computing first diagonal element
     A = fusiontensor(a, b, e)
@@ -152,20 +153,42 @@ TensorKitSectors.Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I <: SUNIrre
         A[1, 2, 3, -1] * B[3, 4, -2]
     return Array(F)
 end
-MemoizationKit.default_cache_settings(::typeof(SUN_Fsymbol)) = (; maxsize = 10^5)
+
+# remove the view calls
+function TensorKitSectors.Asymbol_from_fusiontensor(a::I, b::I, c::I) where {I <: SUNIrrep}
+    Nabc = Nsymbol(a, b, c)
+    iszero(Nabc) && return zeros(fusionscalartype(typeof(a)), Nabc, Nabc)
+
+    C1 = fusiontensor(a, b, c)[:, 1, :, :]
+    C2 = fusiontensor(dual(a), c, b)[:, :, 1, :]
+    Za = sqrtdim(a) * fusiontensor(a, dual(a), leftunit(a))[:, :, 1, 1]
+    @tensor A[-1, -2] := sqrtdim(b) / sqrtdim(c) * conj(Za[1, 2]) * C1[1, 3, -1] * C2[2, 3, -2]
+    return Array(A)
+end
+function TensorKitSectors.Bsymbol_from_fusiontensor(a::I, b::I, c::I) where {I <: SUNIrrep}
+    Nabc = Nsymbol(a, b, c)
+    iszero(Nabc) && return zeros(fusionscalartype(typeof(a)), Nabc, Nabc)
+
+    C1 = fusiontensor(a, b, c)[1, :, :, :]
+    C2 = fusiontensor(c, dual(b), a)[:, :, 1, :]
+    Zb = sqrtdim(b) * fusiontensor(b, dual(b), leftunit(b))[:, :, 1, 1]
+    @tensor B[-1, -2] := sqrtdim(a) / sqrtdim(c) * conj(Zb[1, 2]) * C1[1, 3, -1] * C2[3, 2, -2]
+    return Array(B)
+end
 
 TensorKitSectors.braidingscalartype(::Type{<:SUNIrrep}) = Float64
 TensorKitSectors.Rsymbol(a::I, b::I, c::I) where {I <: SUNIrrep} = SUN_Rsymbol(a, b, c)
 @cached function SUN_Rsymbol(a::I, b::I, c::I)::Array{Float64, 2} where {I <: SUNIrrep}
-    N1 = Nsymbol(a, b, c)
-    N2 = Nsymbol(b, a, c)
-
-    (N1 == 0 || N2 == 0) && return fill(zero(braidingscalartype(I)), N1, N2)
+    return TensorKitSectors.Rsymbol_from_fusiontensor(a, b, c)
+end
+MemoizationKit.default_cache_settings(::typeof(SUN_Rsymbol)) = (; maxsize = 10^5)
+function TensorKitSectors.Rsymbol_from_fusiontensor(a::I, b::I, c::I) where {I <: SUNIrrep}
+    Nabc = Nsymbol(a, b, c)
+    iszero(Nabc) && return fill(zero(braidingscalartype(I)), Nabc, Nabc)
 
     A = fusiontensor(a, b, c)[:, :, 1, :]
     B = fusiontensor(b, a, c)[:, :, 1, :]
 
     @tensor R[-1; -2] := conj(B[1, 2, -2]) * A[2, 1, -1]
     return Array(R)
 end
-MemoizationKit.default_cache_settings(::typeof(SUN_Rsymbol)) = (; maxsize = 10^5)

test/caching.jl

@@ -12,6 +12,8 @@ disk_hits() = only(p.second.hits for p in MemoizationKit.disk_cache_stats() if p
 @testset "CGC disk cache" begin
     # only remove cache if running on CI
     get(ENV, "CI", "false") == "true" && empty_disk_caches!(SUN_CGC)
+    # CGCs still in RAM from earlier tests would not be written to disk
+    empty_caches!(SUN_CGC)
 
     precompute_disk_cache(3, 1)
     @test only(disk_cache_info(SUN_CGC)).second.entries > 0

test/fusiontrees.jl

@@ -1,15 +1,19 @@
+# the trees grow quickly with the dimensions, so use the smallest nontrivial irreps instead of
+# `randsector`, which draws from the larger `smallset` of the sector tests
+fusiontree_sector(::Type{I}) where {I <: SUNIrrep} = rand(filter(!isunit, collect(take(values(I), 5))))
+
 for I in sectorlist
     println("------------------------------------")
     println("Fusion Trees $I")
     println("------------------------------------")
     ti = time()
     N = 5
-    out = ntuple(n -> randsector(I), N)
+    out = ntuple(n -> fusiontree_sector(I), N)
     isdual = ntuple(n -> rand(Bool), N)
     in = rand(collect(⊗(out...)))
     numtrees = count(n -> true, fusiontrees(out, in, isdual))
     while !(0 < numtrees < 30)
-        out = ntuple(n -> randsector(I), N)
+        out = ntuple(n -> fusiontree_sector(I), N)
         in = rand(collect(⊗(out...)))
         numtrees = count(n -> true, fusiontrees(out, in, isdual))
     end
@@ -33,12 +37,12 @@ for I in sectorlist
 
     @timedtestset "Fusion tree $I: insertat" begin
         N = 3
-        out2 = ntuple(n -> randsector(I), N)
+        out2 = ntuple(n -> fusiontree_sector(I), N)
         in2 = rand(collect(⊗(out2...)))
         isdual2 = ntuple(n -> rand(Bool), N)
         f2 = rand(collect(fusiontrees(out2, in2, isdual2)))
         for i in 1:N
-            out1 = ntuple(n -> randsector(I), N)
+            out1 = ntuple(n -> fusiontree_sector(I), N)
             out1 = Base.setindex(out1, in2, i)
             in1 = rand(collect(⊗(out1...)))
             isdual1 = ntuple(n -> rand(Bool), N)
@@ -69,10 +73,10 @@ for I in sectorlist
 
     @timedtestset "Fusion tree $I: merging" begin
         N = 3
-        out1 = ntuple(n -> randsector(I), N)
+        out1 = ntuple(n -> fusiontree_sector(I), N)
         in1 = rand(collect(⊗(out1...)))
         f1 = rand(collect(fusiontrees(out1, in1)))
-        out2 = ntuple(n -> randsector(I), N)
+        out2 = ntuple(n -> fusiontree_sector(I), N)
         in2 = rand(collect(⊗(out2...)))
         f2 = rand(collect(fusiontrees(out2, in2)))
 
@@ -109,10 +113,10 @@ for I in sectorlist
     end
 
     N = 3
-    out = ntuple(n -> randsector(I), N)
+    out = ntuple(n -> fusiontree_sector(I), N)
     numtrees = count(n -> true, fusiontrees((out..., map(dual, out)...)))
     while !(0 < numtrees < 100)
-        out = ntuple(n -> randsector(I), N)
+        out = ntuple(n -> fusiontree_sector(I), N)
         numtrees = count(n -> true, fusiontrees((out..., map(dual, out)...)))
     end
     incoming = rand(collect(⊗(out...)))

test/generic.jl

@@ -36,7 +36,7 @@ using SUNRepresentations.TensorKitSectors: ⊗, Nsymbol, dual
     @test I1.I == @constinferred weight(I1)
 
     @test inv(cartanmatrix(I1)) ≈ inverse_cartanmatrix(I1)
-    @test SUNIrrep{N}("1") === one(SUNIrrep{N})
+    @test SUNIrrep{N}("1") === unit(SUNIrrep{N})
 end
 
 @timedtestset "Direct product of SUNIrrep{$N}:" for N in 2:6

test/runtests.jl

@@ -1,15 +1,16 @@
 using Test
 using TestExtras
-using Aqua
 using Random
+using TensorKitSectors
+using MemoizationKit: MemoizationKit, set_cache_size!
 using TensorKit
 using SUNRepresentations
 using Combinatorics
 using TensorKit
 using TensorKit: ProductSector, fusiontensor, pentagon_equation, hexagon_equation,
     FusionTreeBlock
 using TensorOperations
-using Base.Iterators: take, product
+using Base.Iterators: take
 using SparseArrayKit: SparseArray
 using LinearAlgebra: LinearAlgebra
 
@@ -18,45 +19,59 @@ const TK = TensorKit
 _isone(x; kwargs...) = isapprox(x, one(x); kwargs...)
 
 Random.seed!(1234)
+Ti = time()
+
+# sector tests
+testsuite_path = joinpath(
+    dirname(dirname(pathof(TensorKitSectors))), # TensorKitSectors root
+    "test", "testsuite.jl"
+)
+include(testsuite_path)
+using .SectorTestSuite: randsector, smallset
 
-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
+# the pentagon test scales steeply with the dimensions, so keep these small enough for CI
+function SectorTestSuite.smallset(::Type{I}) where {I <: SUNIrrep}
+    N = SUNRepresentations.rank(I)
+    size = N == 3 ? 5 : N == 4 ? 4 : 3
+    dim_limit = N == 3 ? 11 : 16
+    return smallset(I, size, dim_limit)
 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
+function SectorTestSuite.smallset(::Type{ProductSector{Tuple{I1, I2}}}) where {I1 <: SUNIrrep, I2 <: SUNIrrep}
+    products = [a ⊠ b for a in smallset(I1) for b in smallset(I2) if dim(a) * dim(b) <= 30]
+    return first(shuffle!(products), 5)
 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)
+
+sectorlist = (SU3Irrep, SU4Irrep, SU5Irrep, SU3Irrep ⊠ SU3Irrep)
+# the sector tests visit many intermediate irreps, so keep the RAM caches small (in bytes)
+set_cache_size!(SUNRepresentations.SUN_CGC, 2^30; by = MemoizationKit.cachesize)
+set_cache_size!(SUNRepresentations.SUN_Fsymbol, 2^28; by = MemoizationKit.cachesize)
+set_cache_size!(SUNRepresentations.SUN_Rsymbol, 2^26; by = MemoizationKit.cachesize)
+@testset "Sector tests" begin
+    for sector in sectorlist
+        SectorTestSuite.test_sector(sector)
     end
-    return a
+    include("sectors.jl")
 end
 
-Ti = time()
 module GenericTests
     using Test
     using TestExtras
     using Random
     using SUNRepresentations
+    using TensorKitSectors
     include("generic.jl")
 end
 
-include("caching.jl")
+@testset "Caching tests" begin
+    include("caching.jl")
+end
 include("casimir.jl")
 include("clebschgordan.jl")
-sectorlist = (SU3Irrep, SU4Irrep, SU5Irrep, SU3Irrep ⊠ SU3Irrep)
-include("sectors.jl")
 sectorlist = (SU3Irrep, SU4Irrep, SU5Irrep)
 include("fusiontrees.jl")
 
 @testset "Aqua" verbose = true begin
+    using Aqua
     # RationalRoots has ambiguities with Base/Core, so only test SUNRepresentations ambiguities
     # Intentional piracy of Rep[SU{N}] etc
     Aqua.test_all(SUNRepresentations; ambiguities = false, piracies = (; treat_as_own = [SU]))

test/sectors.jl

@@ -2,110 +2,19 @@ for I in sectorlist
     println("------------------------------------")
     println("Sector $I")
     println("------------------------------------")
-    ti = time()
-    @testset "Sector $I: Basic properties" begin
+    @testset "Sector $I: Additional basic properties" begin
         s = (randsector(I), randsector(I), randsector(I))
 
         mode_old = SUNRepresentations.display_mode("dimension")
         for mode in ["dimension", "dynkin", "weight"]
             SUNRepresentations.display_mode(mode)
-            @test eval(Meta.parse(sprint(show, s[1]))) == s[1]
         end
         SUNRepresentations.display_mode(mode_old)
-
-        @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 frobenius_schur_phase(s[1])
-        @constinferred frobenius_schur_indicator(s[1])
-        @constinferred Nsymbol(s...)
-        @constinferred Rsymbol(s...)
-        @constinferred Bsymbol(s...)
-        @constinferred Fsymbol(s..., s...)
-        it = @constinferred s[1] ⊗ s[2]
-        @constinferred ⊗(s..., s...)
         for i in 1:3
-            @test 1 == @constinferred Nsymbol(s[i], conj(s[i]), one(s[i]))
-        end
-    end
-    @timedtestset "Sector $I: Value iterator" begin
-        @test eltype(values(I)) == I
-        sprev = one(I)
-        for (i, s) in enumerate(values(I))
-            @test !isless(s, sprev) # confirm compatibility with sort order
-            @test s == @constinferred (values(I)[i])
-            @test TensorKit.findindex(values(I), s) == i
-            sprev = s
-            i >= 10 && break
-        end
-        @test one(I) == first(values(I))
-        @test (@constinferred TensorKit.findindex(values(I), one(I))) == 1
-        for s in smallset(I)
-            @test (@constinferred values(I)[TensorKit.findindex(values(I), s)]) == s
+            @test 1 == @constinferred Nsymbol(s[i], dual(s[i]), unit(s[i]))
         end
         @test all(i -> TensorKit.findindex(values(I), values(I)[i]) == i, 1:1000)
     end
-    @timedtestset "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]
-                    f2 = Fsymbol(a, b, c, d, e, f) * dim(d)
-                    @test isapprox(f1, f2; atol = 1000 * eps(), rtol = 1000 * eps())
-                end
-            end
-        end
-    end
-    @timedtestset "Sector $I: Unitarity of F-move" begin
-        for a in smallset(I), b in smallset(I), c in smallset(I), 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))))
-            Fblocks = Vector{Any}()
-            for e in es, 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
-            F = hvcat(length(fs), Fblocks...)
-            @test F' * F ≈ one(F)
-        end
-    end
-    @testset "Sector $I: 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 = 1.0e-12, rtol = 1.0e-12)
-        end
-    end
-    @testset "Sector $I: Hexagon equation" begin
-        for a in smallset(I), b in smallset(I), c in smallset(I)
-            @test hexagon_equation(a, b, c; atol = 1.0e-12, rtol = 1.0e-12)
-        end
-    end
-    tf = time()
-    printstyled(
-        "Finished sector $I tests in ",
-        string(round(tf - ti; sigdigits = 3)),
-        " seconds."; bold = true, color = Base.info_color()
-    )
-    println()
 
     try
         s = sprint(SUNRepresentations.cache_info)