diff --git a/NDTensors/src/imports.jl b/NDTensors/src/imports.jl index 04ea3d8d6b..7453b1f20e 100644 --- a/NDTensors/src/imports.jl +++ b/NDTensors/src/imports.jl @@ -33,9 +33,9 @@ for lib in [ :MetalExtensions, :BroadcastMapConversion, :RankFactorization, - :Sectors, :LabelledNumbers, :GradedAxes, + :Sectors, :TensorAlgebra, :SparseArrayInterface, :SparseArrayDOKs, diff --git a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/Project.toml b/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/Project.toml deleted file mode 100644 index 9b1d5ccd25..0000000000 --- a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/Project.toml +++ /dev/null @@ -1,2 +0,0 @@ -[deps] -NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf" diff --git a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/src/GradedAxesSectorsExt.jl b/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/src/GradedAxesSectorsExt.jl deleted file mode 100644 index aa3056438e..0000000000 --- a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/src/GradedAxesSectorsExt.jl +++ /dev/null @@ -1,9 +0,0 @@ -module GradedAxesSectorsExt -using ..GradedAxes: GradedAxes -using ...Sectors: Sectors, AbstractCategory, ⊗ # , dual - -GradedAxes.fuse_labels(c1::AbstractCategory, c2::AbstractCategory) = only(c1 ⊗ c2) - -# TODO: Decide the fate of `dual`. -## GradedAxes.dual(c::AbstractCategory) = dual(c) -end diff --git a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/Project.toml b/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/Project.toml deleted file mode 100644 index ef491a529c..0000000000 --- a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/Project.toml +++ /dev/null @@ -1,3 +0,0 @@ -[deps] -NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/runtests.jl b/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/runtests.jl deleted file mode 100644 index 371e7e57cd..0000000000 --- a/NDTensors/src/lib/GradedAxes/ext/GradedAxesSectorsExt/test/runtests.jl +++ /dev/null @@ -1,15 +0,0 @@ -@eval module $(gensym()) -using NDTensors.GradedAxes: dual, fuse_labels -using NDTensors.Sectors: U1, Z -using Test: @test, @testset - -@testset "GradedAxesSectorsExt" begin - @test fuse_labels(U1(1), U1(2)) == U1(3) - @test dual(U1(2)) == U1(-2) - - @test fuse_labels(Z{2}(1), Z{2}(1)) == Z{2}(0) - @test fuse_labels(Z{2}(0), Z{2}(1)) == Z{2}(1) - @test dual(Z{2}(1)) == Z{2}(1) - @test dual(Z{2}(0)) == Z{2}(0) -end -end diff --git a/NDTensors/src/lib/GradedAxes/src/GradedAxes.jl b/NDTensors/src/lib/GradedAxes/src/GradedAxes.jl index 7756b71575..d605289f81 100644 --- a/NDTensors/src/lib/GradedAxes/src/GradedAxes.jl +++ b/NDTensors/src/lib/GradedAxes/src/GradedAxes.jl @@ -1,8 +1,7 @@ module GradedAxes include("blockedunitrange.jl") include("gradedunitrange.jl") -include("fusion.jl") include("dual.jl") include("unitrangedual.jl") -include("../ext/GradedAxesSectorsExt/src/GradedAxesSectorsExt.jl") +include("fusion.jl") end diff --git a/NDTensors/src/lib/GradedAxes/src/dual.jl b/NDTensors/src/lib/GradedAxes/src/dual.jl index 985ebe33cc..ff11b1103c 100644 --- a/NDTensors/src/lib/GradedAxes/src/dual.jl +++ b/NDTensors/src/lib/GradedAxes/src/dual.jl @@ -5,3 +5,5 @@ using NDTensors.LabelledNumbers: label_dual(x) = label_dual(LabelledStyle(x), x) label_dual(::NotLabelled, x) = x label_dual(::IsLabelled, x) = labelled(unlabel(x), dual(label(x))) + +flip(g::AbstractGradedUnitRange) = dual(gradedrange(label_dual.(blocklengths(g)))) diff --git a/NDTensors/src/lib/GradedAxes/src/fusion.jl b/NDTensors/src/lib/GradedAxes/src/fusion.jl index 5193449f2f..2e5513c7ee 100644 --- a/NDTensors/src/lib/GradedAxes/src/fusion.jl +++ b/NDTensors/src/lib/GradedAxes/src/fusion.jl @@ -6,6 +6,10 @@ OneToOne() = OneToOne{Bool}() Base.first(a::OneToOne) = one(eltype(a)) Base.last(a::OneToOne) = one(eltype(a)) +gradedisequal(::AbstractUnitRange, ::OneToOne) = false +gradedisequal(::OneToOne, ::AbstractUnitRange) = false +gradedisequal(::OneToOne, ::OneToOne) = true + # https://github.com/ITensor/ITensors.jl/blob/v0.3.57/NDTensors/src/lib/GradedAxes/src/tensor_product.jl # https://en.wikipedia.org/wiki/Tensor_product # https://github.com/KeitaNakamura/Tensorial.jl @@ -18,7 +22,7 @@ function tensor_product( return foldl(tensor_product, (a1, a2, a3, a_rest...)) end -function tensor_product(a1::AbstractUnitRange, a2::AbstractUnitRange) +function tensor_product(::AbstractUnitRange, ::AbstractUnitRange) return error("Not implemented yet.") end @@ -26,18 +30,31 @@ function tensor_product(a1::Base.OneTo, a2::Base.OneTo) return Base.OneTo(length(a1) * length(a2)) end -function tensor_product(a1::OneToOne, a2::AbstractUnitRange) +function tensor_product(::OneToOne, a2::AbstractUnitRange) return a2 end -function tensor_product(a1::AbstractUnitRange, a2::OneToOne) +function tensor_product(a1::AbstractUnitRange, ::OneToOne) return a1 end -function tensor_product(a1::OneToOne, a2::OneToOne) +function tensor_product(::OneToOne, ::OneToOne) return OneToOne() end +# Handle dual. Always return a non-dual GradedUnitRange. +function tensor_product(a1::AbstractUnitRange, a2::UnitRangeDual) + return tensor_product(a1, flip(a2)) +end + +function tensor_product(a1::UnitRangeDual, a2::AbstractUnitRange) + return tensor_product(flip(a1), a2) +end + +function tensor_product(a1::UnitRangeDual, a2::UnitRangeDual) + return tensor_product(flip(a1), flip(a2)) +end + function fuse_labels(x, y) return error( "`fuse_labels` not implemented for object of type `$(typeof(x))` and `$(typeof(y))`." @@ -53,21 +70,28 @@ function fuse_blocklengths(x::LabelledInteger, y::LabelledInteger) return labelled(unlabel(x) * unlabel(y), fuse_labels(label(x), label(y))) end +flatten_maybe_nested(v::Vector{<:Integer}) = v +flatten_maybe_nested(v::Vector{<:AbstractGradedUnitRange}) = reduce(vcat, blocklengths.(v)) + using BlockArrays: blockedrange, blocks function tensor_product(a1::AbstractBlockedUnitRange, a2::AbstractBlockedUnitRange) - blocklengths = map(vec(collect(Iterators.product(blocks(a1), blocks(a2))))) do x - return mapreduce(length, fuse_blocklengths, x) - end + maybe_nested = map( + it -> mapreduce(length, fuse_blocklengths, it), + Iterators.flatten((Iterators.product(blocks(a1), blocks(a2)),)), + ) + blocklengths = flatten_maybe_nested(maybe_nested) return blockedrange(blocklengths) end function blocksortperm(a::AbstractBlockedUnitRange) - # TODO: Figure out how to deal with dual sectors. - # TODO: `rev=isdual(a)` may not be correct for symmetries beyond `U(1)`. - ## return Block.(sortperm(nondual_sectors(a); rev=isdual(a))) return Block.(sortperm(blocklabels(a))) end +# convention: sort UnitRangeDual according to nondual blocks +function blocksortperm(a::UnitRangeDual) + return Block.(sortperm(blocklabels(nondual(a)))) +end + using BlockArrays: Block, BlockVector using SplitApplyCombine: groupcount # Get the permutation for sorting, then group by common elements. @@ -83,24 +107,38 @@ end # Get the permutation for sorting, then group by common elements. # groupsortperm([2, 1, 2, 3]) == [[2], [1, 3], [4]] function blockmergesortperm(a::AbstractBlockedUnitRange) - # If it is dual, reverse the sorting so the sectors - # end up sorted in the same way whether or not the space - # is dual. - # TODO: Figure out how to deal with dual sectors. - # TODO: `rev=isdual(a)` may not be correct for symmetries beyond `U(1)`. - ## return Block.(groupsortperm(nondual_sectors(a); rev=isdual(a))) return Block.(groupsortperm(blocklabels(a))) end # Used by `TensorAlgebra.splitdims` in `BlockSparseArraysGradedAxesExt`. invblockperm(a::Vector{<:Block{1}}) = Block.(invperm(Int.(a))) -# Used by `TensorAlgebra.fusedims` in `BlockSparseArraysGradedAxesExt`. -function blockmergesortperm(a::GradedUnitRange) - # If it is dual, reverse the sorting so the sectors - # end up sorted in the same way whether or not the space - # is dual. - # TODO: Figure out how to deal with dual sectors. - # TODO: `rev=isdual(a)` may not be correct for symmetries beyond `U(1)`. - return Block.(groupsortperm(blocklabels(a))) +function blockmergesortperm(a::UnitRangeDual) + return Block.(groupsortperm(blocklabels(nondual(a)))) +end + +function blockmergesort(g::AbstractGradedUnitRange) + glabels = blocklabels(g) + gblocklengths = blocklengths(g) + new_blocklengths = map( + la -> labelled(sum(gblocklengths[findall(==(la), glabels)]; init=0), la), + sort(unique(glabels)), + ) + return GradedAxes.gradedrange(new_blocklengths) +end + +blockmergesort(g::UnitRangeDual) = dual(blockmergesort(flip(g))) +blockmergesort(g::OneToOne) = g + +# fusion_product produces a sorted, non-dual GradedUnitRange +function fusion_product(g1, g2) + return blockmergesort(tensor_product(g1, g2)) +end + +fusion_product(g::AbstractUnitRange) = blockmergesort(g) +fusion_product(g::UnitRangeDual) = fusion_product(flip(g)) + +# recursive fusion_product. Simpler than reduce + fix type stability issues with reduce +function fusion_product(g1, g2, g3...) + return fusion_product(fusion_product(g1, g2), g3...) end diff --git a/NDTensors/src/lib/GradedAxes/src/gradedunitrange.jl b/NDTensors/src/lib/GradedAxes/src/gradedunitrange.jl index 57e9420d88..5961d4098d 100644 --- a/NDTensors/src/lib/GradedAxes/src/gradedunitrange.jl +++ b/NDTensors/src/lib/GradedAxes/src/gradedunitrange.jl @@ -12,7 +12,7 @@ using BlockArrays: blockedrange, BlockIndexRange, blockfirsts, - blocklasts, + blockisequal, blocklength, blocklengths, findblock, @@ -37,6 +37,11 @@ function Base.OrdinalRange{T,T}(a::GradedOneTo{<:LabelledInteger{T}}) where {T} return unlabel_blocks(a) end +# == is just a range comparison that ignores labels. Need dedicated function to check equality. +function gradedisequal(a1::AbstractUnitRange, a2::AbstractUnitRange) + return blockisequal(a1, a2) && (blocklabels(a1) == blocklabels(a2)) +end + # This is only needed in certain Julia versions below 1.10 # (for example Julia 1.6). # TODO: Delete this once we drop Julia 1.6 support. diff --git a/NDTensors/src/lib/GradedAxes/src/unitrangedual.jl b/NDTensors/src/lib/GradedAxes/src/unitrangedual.jl index aa04cc1600..3f5249728c 100644 --- a/NDTensors/src/lib/GradedAxes/src/unitrangedual.jl +++ b/NDTensors/src/lib/GradedAxes/src/unitrangedual.jl @@ -6,7 +6,10 @@ UnitRangeDual(a::AbstractUnitRange) = UnitRangeDual{eltype(a),typeof(a)}(a) dual(a::AbstractUnitRange) = UnitRangeDual(a) nondual(a::UnitRangeDual) = a.nondual_unitrange dual(a::UnitRangeDual) = nondual(a) +flip(a::UnitRangeDual) = dual(flip(nondual(a))) nondual(a::AbstractUnitRange) = a +isdual(::AbstractGradedUnitRange) = false +isdual(::UnitRangeDual) = true ## TODO: Define this to instantiate a dual unit range. ## materialize_dual(a::UnitRangeDual) = materialize_dual(nondual(a)) @@ -16,6 +19,16 @@ Base.step(a::UnitRangeDual) = label_dual(step(nondual(a))) Base.view(a::UnitRangeDual, index::Block{1}) = a[index] +function Base.show(io::IO, a::UnitRangeDual) + return print(io, UnitRangeDual, "(", blocklasts(a), ")") +end + +function Base.show(io::IO, mimetype::MIME"text/plain", a::UnitRangeDual) + return Base.invoke( + show, Tuple{typeof(io),MIME"text/plain",AbstractArray}, io, mimetype, a + ) +end + function Base.getindex(a::UnitRangeDual, indices::AbstractUnitRange{<:Integer}) return dual(getindex(nondual(a), indices)) end @@ -92,6 +105,14 @@ BlockArrays.blockaxes(a::UnitRangeDual) = blockaxes(nondual(a)) BlockArrays.blockfirsts(a::UnitRangeDual) = label_dual.(blockfirsts(nondual(a))) BlockArrays.blocklasts(a::UnitRangeDual) = label_dual.(blocklasts(nondual(a))) BlockArrays.findblock(a::UnitRangeDual, index::Integer) = findblock(nondual(a), index) + +blocklabels(a::UnitRangeDual) = dual.(blocklabels(nondual(a))) + +gradedisequal(::UnitRangeDual, ::AbstractGradedUnitRange) = false +gradedisequal(::AbstractGradedUnitRange, ::UnitRangeDual) = false +function gradedisequal(a1::UnitRangeDual, a2::UnitRangeDual) + return gradedisequal(nondual(a1), nondual(a2)) +end function BlockArrays.combine_blockaxes(a1::UnitRangeDual, a2::UnitRangeDual) return dual(combine_blockaxes(dual(a1), dual(a2))) end diff --git a/NDTensors/src/lib/GradedAxes/test/runtests.jl b/NDTensors/src/lib/GradedAxes/test/runtests.jl index 09335af5e8..c0fdca21be 100644 --- a/NDTensors/src/lib/GradedAxes/test/runtests.jl +++ b/NDTensors/src/lib/GradedAxes/test/runtests.jl @@ -2,7 +2,7 @@ using Test: @testset @testset "GradedAxes" begin include("test_basics.jl") - include("test_tensor_product.jl") include("test_dual.jl") + include("test_tensor_product.jl") end end diff --git a/NDTensors/src/lib/GradedAxes/test/test_basics.jl b/NDTensors/src/lib/GradedAxes/test/test_basics.jl index e1dcc67174..2f3cd03bc7 100644 --- a/NDTensors/src/lib/GradedAxes/test/test_basics.jl +++ b/NDTensors/src/lib/GradedAxes/test/test_basics.jl @@ -9,7 +9,8 @@ using BlockArrays: blocklength, blocklengths, blocks -using NDTensors.GradedAxes: GradedOneTo, GradedUnitRange, blocklabels, gradedrange +using NDTensors.GradedAxes: + GradedOneTo, GradedUnitRange, blocklabels, gradedisequal, gradedrange using NDTensors.LabelledNumbers: LabelledUnitRange, islabelled, label, labelled, unlabel using Test: @test, @test_broken, @testset @testset "GradedAxes basics" begin @@ -40,6 +41,7 @@ using Test: @test, @test_broken, @testset @test label(x) == "y" end @test isnothing(iterate(a, labelled(5, "y"))) + @test gradedisequal(a, a) @test length(a) == 5 @test step(a) == 1 @test !islabelled(step(a)) diff --git a/NDTensors/src/lib/GradedAxes/test/test_dual.jl b/NDTensors/src/lib/GradedAxes/test/test_dual.jl index 0fb15d31bf..0cce88a6dd 100644 --- a/NDTensors/src/lib/GradedAxes/test/test_dual.jl +++ b/NDTensors/src/lib/GradedAxes/test/test_dual.jl @@ -1,19 +1,40 @@ @eval module $(gensym()) -using BlockArrays: Block, blockaxes, blockfirsts, blocklasts, blocks, findblock -using NDTensors.GradedAxes: GradedAxes, UnitRangeDual, dual, gradedrange, nondual +using BlockArrays: + Block, blockaxes, blockfirsts, blocklasts, blocklength, blocklengths, blocks, findblock +using NDTensors.GradedAxes: + GradedAxes, + UnitRangeDual, + blocklabels, + blockmergesortperm, + blocksortperm, + dual, + flip, + gradedisequal, + gradedrange, + isdual, + nondual using NDTensors.LabelledNumbers: LabelledInteger, label, labelled using Test: @test, @test_broken, @testset struct U1 n::Int end GradedAxes.dual(c::U1) = U1(-c.n) +Base.isless(c1::U1, c2::U1) = c1.n < c2.n @testset "dual" begin a = gradedrange([U1(0) => 2, U1(1) => 3]) ad = dual(a) @test eltype(ad) == LabelledInteger{Int,U1} - @test dual(ad) == a - @test nondual(ad) == a - @test nondual(a) == a + + @test gradedisequal(dual(ad), a) + @test gradedisequal(nondual(ad), a) + @test gradedisequal(nondual(a), a) + @test gradedisequal(ad, ad) + @test !gradedisequal(a, ad) + @test !gradedisequal(ad, a) + + @test isdual(ad) + @test !isdual(a) + @test blockfirsts(ad) == [labelled(1, U1(0)), labelled(3, U1(-1))] @test blocklasts(ad) == [labelled(2, U1(0)), labelled(5, U1(-1))] @test findblock(ad, 4) == Block(2) @@ -34,5 +55,36 @@ GradedAxes.dual(c::U1) = U1(-c.n) @test label(ad[[Block(2), Block(1)]][Block(1)]) == U1(-1) @test ad[[Block(2)[1:2], Block(1)[1:2]]][Block(1)] == 3:4 @test label(ad[[Block(2)[1:2], Block(1)[1:2]]][Block(1)]) == U1(-1) + @test blocksortperm(a) == [Block(1), Block(2)] + @test blocksortperm(ad) == [Block(1), Block(2)] + @test blocklength(blockmergesortperm(a)) == 2 + @test blocklength(blockmergesortperm(ad)) == 2 + @test blockmergesortperm(a) == [Block(1), Block(2)] + @test blockmergesortperm(ad) == [Block(1), Block(2)] +end + +@testset "flip" begin + a = gradedrange([U1(0) => 2, U1(1) => 3]) + ad = dual(a) + @test gradedisequal(flip(a), dual(gradedrange([U1(0) => 2, U1(-1) => 3]))) + @test gradedisequal(flip(ad), gradedrange([U1(0) => 2, U1(-1) => 3])) + + @test blocklabels(a) == [U1(0), U1(1)] + @test blocklabels(dual(a)) == [U1(0), U1(-1)] + @test blocklabels(flip(a)) == [U1(0), U1(1)] + @test blocklabels(flip(dual(a))) == [U1(0), U1(-1)] + @test blocklabels(dual(flip(a))) == [U1(0), U1(-1)] + + @test blocklengths(a) == [2, 3] + @test blocklengths(dual(a)) == [2, 3] + @test blocklengths(flip(a)) == [2, 3] + @test blocklengths(flip(dual(a))) == [2, 3] + @test blocklengths(dual(flip(a))) == [2, 3] + + @test !isdual(a) + @test isdual(dual(a)) + @test isdual(flip(a)) + @test !isdual(flip(dual(a))) + @test !isdual(dual(flip(a))) end end diff --git a/NDTensors/src/lib/GradedAxes/test/test_tensor_product.jl b/NDTensors/src/lib/GradedAxes/test/test_tensor_product.jl index c10ae4bf95..7b533f79c5 100644 --- a/NDTensors/src/lib/GradedAxes/test/test_tensor_product.jl +++ b/NDTensors/src/lib/GradedAxes/test/test_tensor_product.jl @@ -1,11 +1,83 @@ @eval module $(gensym()) -using NDTensors.GradedAxes: GradedAxes, GradedOneTo, gradedrange, tensor_product using Test: @test, @testset + +using BlockArrays: blocklength, blocklengths + +using NDTensors.GradedAxes: + GradedAxes, + GradedOneTo, + OneToOne, + dual, + fusion_product, + flip, + gradedrange, + gradedisequal, + isdual, + tensor_product + +struct U1 + n::Int +end +GradedAxes.dual(c::U1) = U1(-c.n) +Base.isless(c1::U1, c2::U1) = c1.n < c2.n +GradedAxes.fuse_labels(x::U1, y::U1) = U1(x.n + y.n) + @testset "GradedAxes.tensor_product" begin GradedAxes.fuse_labels(x::String, y::String) = x * y + + g0 = OneToOne() + @test gradedisequal(tensor_product(g0, g0), g0) + a = gradedrange(["x" => 2, "y" => 3]) b = tensor_product(a, a) + @test b isa GradedOneTo @test length(b) == 25 + @test blocklength(b) == 4 + @test blocklengths(b) == [4, 6, 6, 9] + @test gradedisequal(b, gradedrange(["xx" => 4, "yx" => 6, "xy" => 6, "yy" => 9])) + + c = tensor_product(a, a, a) + @test c isa GradedOneTo + @test length(c) == 125 + @test blocklength(c) == 8 +end + +@testset "GradedAxes.fusion_product" begin + g0 = OneToOne() + @test gradedisequal(fusion_product(g0, g0), g0) + + a = gradedrange([U1(1) => 1, U1(2) => 3, U1(1) => 1]) + + b = fusion_product(a) + @test gradedisequal(b, gradedrange([U1(1) => 2, U1(2) => 3])) + + c = fusion_product(a, a) + @test gradedisequal(c, gradedrange([U1(2) => 4, U1(3) => 12, U1(4) => 9])) + + d = fusion_product(a, a, a) + @test gradedisequal(d, gradedrange([U1(3) => 8, U1(4) => 36, U1(5) => 54, U1(6) => 27])) +end + +@testset "dual and tensor_product" begin + a = gradedrange([U1(1) => 1, U1(2) => 3, U1(1) => 1]) + ad = dual(a) + + b = fusion_product(ad) @test b isa GradedOneTo + @test !isdual(b) + @test gradedisequal(b, gradedrange([U1(-2) => 3, U1(-1) => 2])) + + c = fusion_product(ad, ad) + @test c isa GradedOneTo + @test !isdual(c) + @test gradedisequal(c, gradedrange([U1(-4) => 9, U1(-3) => 12, U1(-2) => 4])) + + d = fusion_product(ad, a) + @test !isdual(d) + @test gradedisequal(d, gradedrange([U1(-1) => 6, U1(0) => 13, U1(1) => 6])) + + e = fusion_product(a, ad) + @test !isdual(d) + @test gradedisequal(e, d) end end diff --git a/NDTensors/src/lib/Sectors/src/Sectors.jl b/NDTensors/src/lib/Sectors/src/Sectors.jl index fc813e602e..f851dc4002 100644 --- a/NDTensors/src/lib/Sectors/src/Sectors.jl +++ b/NDTensors/src/lib/Sectors/src/Sectors.jl @@ -1,13 +1,13 @@ module Sectors +include("symmetry_style.jl") include("abstractcategory.jl") -include("category_definitions/u1.jl") -include("category_definitions/zn.jl") -include("category_definitions/su.jl") -include("category_definitions/su2.jl") -include("category_definitions/su2k.jl") include("category_definitions/fib.jl") include("category_definitions/ising.jl") +include("category_definitions/su.jl") +include("category_definitions/su2k.jl") +include("category_definitions/u1.jl") +include("category_definitions/zn.jl") include("namedtuple_operations.jl") include("category_product.jl") diff --git a/NDTensors/src/lib/Sectors/src/abstractcategory.jl b/NDTensors/src/lib/Sectors/src/abstractcategory.jl index 4aaf56c6c5..51311ab201 100644 --- a/NDTensors/src/lib/Sectors/src/abstractcategory.jl +++ b/NDTensors/src/lib/Sectors/src/abstractcategory.jl @@ -1,45 +1,121 @@ +# This file defines the abstract type AbstractCategory +# all fusion categories (Z{2}, SU2, Ising...) are subtypes of AbstractCategory + abstract type AbstractCategory end -label(c::AbstractCategory) = error("method `label` not defined for type $(typeof(c))") +# =================================== Base interface ===================================== +function Base.isless(c1::C, c2::C) where {C<:AbstractCategory} + return isless(category_label(c1), category_label(c2)) +end -function dimension(c::AbstractCategory) - return error("method `dimension` not defined for type $(typeof(c))") +# ================================= Sectors interface ==================================== +trivial(x) = trivial(typeof(x)) +function trivial(axis_type::Type{<:AbstractUnitRange}) + return GradedAxes.gradedrange([trivial(eltype(axis_type))]) # always returns nondual +end +function trivial(la_type::Type{<:LabelledNumbers.LabelledInteger}) + return la_type(1, trivial(LabelledNumbers.label_type(la_type))) +end +function trivial(type::Type) + return error("`trivial` not defined for type $(type).") end -function label_fusion_rule(category_type::Type{<:AbstractCategory}, l1, l2) - return error("`label_fusion_rule` not defined for type $(category_type).") +istrivial(c::AbstractCategory) = (c == trivial(c)) + +function category_label(c::AbstractCategory) + return error("method `category_label` not defined for type $(typeof(c))") +end + +block_dimensions(g::AbstractUnitRange) = block_dimensions(SymmetryStyle(g), g) +block_dimensions(::AbelianGroup, g) = GradedAxes.unlabel.(BlockArrays.blocklengths(g)) +function block_dimensions(::SymmetryStyle, g) + return Sectors.quantum_dimension.(GradedAxes.blocklabels(g)) .* + BlockArrays.blocklengths(g) end -function fusion_rule(c1::AbstractCategory, c2::AbstractCategory) - category_type = typeof(c1) - return [category_type(l) for l in label_fusion_rule(category_type, label(c1), label(c2))] +quantum_dimension(x) = quantum_dimension(SymmetryStyle(x), x) + +function quantum_dimension(::SymmetryStyle, c::AbstractCategory) + return error("method `quantum_dimension` not defined for type $(typeof(c))") end +quantum_dimension(::AbelianGroup, ::AbstractCategory) = 1 +quantum_dimension(::EmptyCategory, ::AbstractCategory) = 1 +quantum_dimension(::SymmetryStyle, g::AbstractUnitRange) = sum(block_dimensions(g)) +quantum_dimension(::AbelianGroup, g::AbstractUnitRange) = length(g) + +# =============================== Fusion rule interface ================================== ⊗(c1::AbstractCategory, c2::AbstractCategory) = fusion_rule(c1, c2) -⊕(c1::AbstractCategory, c2::AbstractCategory) = [c1, c2] -⊕(cs::Vector{<:AbstractCategory}, c::AbstractCategory) = [cs; c] -⊕(c::AbstractCategory, cs::Vector{<:AbstractCategory}) = [c; cs] +function fusion_rule(c1, c2) + return fusion_rule(combine_styles(SymmetryStyle(c1), SymmetryStyle(c2)), c1, c2) +end -function Base.show(io::IO, cs::Vector{<:AbstractCategory}) - (length(cs) <= 1) && print(io, "[") - symbol = "" - for c in cs - print(io, symbol, c) - symbol = " ⊕ " - end - (length(cs) <= 1) && print(io, "]") - return nothing +function fusion_rule(::SymmetryStyle, c1::C, c2::C) where {C<:AbstractCategory} + degen, labels = label_fusion_rule(C, category_label(c1), category_label(c2)) + return GradedAxes.gradedrange(LabelledNumbers.labelled.(degen, C.(labels))) end -function trivial(category_type::Type{<:AbstractCategory}) - return error("`trivial` not defined for type $(category_type).") +# abelian case: return Category +function fusion_rule(::AbelianGroup, c1::C, c2::C) where {C<:AbstractCategory} + return C(label_fusion_rule(C, category_label(c1), category_label(c2))) end -istrivial(c::AbstractCategory) = (c == trivial(typeof(c))) +function fusion_rule( + ::SymmetryStyle, l1::LabelledNumbers.LabelledInteger, l2::LabelledNumbers.LabelledInteger +) + fused = LabelledNumbers.label(l1) ⊗ LabelledNumbers.label(l2) + v = + LabelledNumbers.labelled.( + l1 * l2 .* BlockArrays.blocklengths(fused), GradedAxes.blocklabels(fused) + ) + return GradedAxes.gradedrange(v) +end + +function fusion_rule( + ::AbelianGroup, l1::LabelledNumbers.LabelledInteger, l2::LabelledNumbers.LabelledInteger +) + fused = LabelledNumbers.label(l1) ⊗ LabelledNumbers.label(l2) + return LabelledNumbers.labelled(l1 * l2, fused) +end -function dual(category_type::Type{<:AbstractCategory}) - return error("`dual` not defined for type $(category_type).") +function fusion_rule( + ::EmptyCategory, l1::LabelledNumbers.LabelledInteger, l2::LabelledNumbers.LabelledInteger +) + return LabelledNumbers.labelled(l1 * l2, sector()) end -Base.isless(c1::AbstractCategory, c2::AbstractCategory) = isless(label(c1), label(c2)) +function label_fusion_rule(category_type::Type{<:AbstractCategory}, ::Any, ::Any) + return error("`label_fusion_rule` not defined for type $(category_type).") +end + +# ================================ GradedAxes interface ================================== +# tensor_product interface +function GradedAxes.fuse_blocklengths( + l1::LabelledNumbers.LabelledInteger{<:Integer,<:Sectors.AbstractCategory}, + l2::LabelledNumbers.LabelledInteger{<:Integer,<:Sectors.AbstractCategory}, +) + return fusion_rule(l1, l2) +end + +# cast to range +to_graded_axis(c::AbstractCategory) = to_graded_axis(LabelledNumbers.labelled(1, c)) +to_graded_axis(l::LabelledNumbers.LabelledInteger) = GradedAxes.gradedrange([l]) +to_graded_axis(g::AbstractUnitRange) = g + +# allow to fuse a category with a GradedUnitRange +function GradedAxes.tensor_product(c::AbstractCategory, g::AbstractUnitRange) + return GradedAxes.tensor_product(to_graded_axis(c), g) +end + +function GradedAxes.tensor_product(g::AbstractUnitRange, c::AbstractCategory) + return GradedAxes.tensor_product(c, g) +end + +function GradedAxes.tensor_product(c1::AbstractCategory, c2::AbstractCategory) + return to_graded_axis(fusion_rule(c1, c2)) +end + +function GradedAxes.fusion_product(c::AbstractCategory) + return GradedAxes.fusion_product(to_graded_axis(c)) +end diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/fib.jl b/NDTensors/src/lib/Sectors/src/category_definitions/fib.jl index b7b8a89d49..ff7f69795e 100644 --- a/NDTensors/src/lib/Sectors/src/category_definitions/fib.jl +++ b/NDTensors/src/lib/Sectors/src/category_definitions/fib.jl @@ -18,13 +18,15 @@ function Fib(s::AbstractString) return error("Unrecognized input \"$s\" to Fib constructor") end -dual(f::Fib) = f +SymmetryStyle(::Fib) = NonGroupCategory() -label(f::Fib) = f.l +GradedAxes.dual(f::Fib) = f + +category_label(f::Fib) = f.l trivial(::Type{Fib}) = Fib(0) -dimension(f::Fib) = istrivial(f) ? 1 : ((1 + √5) / 2) +quantum_dimension(::NonGroupCategory, f::Fib) = istrivial(f) ? 1.0 : ((1 + √5) / 2) # Fusion rules identical to su2₃ label_fusion_rule(::Type{Fib}, l1, l2) = label_fusion_rule(su2{3}, l1, l2) diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/ising.jl b/NDTensors/src/lib/Sectors/src/category_definitions/ising.jl index 1d1089046b..d390f524c9 100644 --- a/NDTensors/src/lib/Sectors/src/category_definitions/ising.jl +++ b/NDTensors/src/lib/Sectors/src/category_definitions/ising.jl @@ -1,4 +1,4 @@ -using HalfIntegers: Half, twice +using HalfIntegers: HalfIntegers # # Ising category @@ -7,7 +7,7 @@ using HalfIntegers: Half, twice # struct Ising <: AbstractCategory - l::Half{Int} + l::HalfIntegers.Half{Int} end # TODO: Use `Val` dispatch here? @@ -18,19 +18,21 @@ function Ising(s::AbstractString) return error("Unrecognized input \"$s\" to Ising constructor") end -dual(i::Ising) = i +SymmetryStyle(::Ising) = NonGroupCategory() -label(i::Ising) = i.l +GradedAxes.dual(i::Ising) = i + +category_label(i::Ising) = i.l trivial(::Type{Ising}) = Ising(0) -dimension(i::Ising) = (label(i) == 1//2) ? √2 : 1 +quantum_dimension(::NonGroupCategory, i::Ising) = (category_label(i) == 1//2) ? √2 : 1.0 # Fusion rules identical to su2₂ label_fusion_rule(::Type{Ising}, l1, l2) = label_fusion_rule(su2{2}, l1, l2) # TODO: Use `Val` dispatch here? -label_to_str(i::Ising) = ("1", "σ", "ψ")[twice(label(i)) + 1] +label_to_str(i::Ising) = ("1", "σ", "ψ")[HalfIntegers.twice(category_label(i)) + 1] function Base.show(io::IO, f::Ising) return print(io, "Ising(", label_to_str(f), ")") diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/su.jl b/NDTensors/src/lib/Sectors/src/category_definitions/su.jl index 8fc82df329..9937e1c362 100644 --- a/NDTensors/src/lib/Sectors/src/category_definitions/su.jl +++ b/NDTensors/src/lib/Sectors/src/category_definitions/su.jl @@ -1,28 +1,40 @@ # -# Special unitary group SU{N} +# Special unitary group SU(N) # -struct SU{N} <: AbstractCategory +struct SU{N,M} <: AbstractCategory # l is the first row of the # Gelfand-Tsetlin (GT) pattern describing # an SU(N) irrep - #TODO: any way this could be NTuple{N-1,Int} ? - l::NTuple{N,Int} + # this first row is identical to the Young tableau of the irrep + l::NTuple{M,Int} + + # M is there to avoid storing a N-Tuple with an extra zero. + # inner constructor enforces M = N - 1 + # It does NOT check for Young Tableau validity (non-increasing positive integers) + function SU{N,M}(t::NTuple{M,Integer}) where {N,M} + return N == M + 1 && M > 0 ? new{N,M}(t) : error("Invalid tuple length") + end end -label(s::SU) = s.l +SU{N}(t::Tuple) where {N} = SU{N,length(t)}(t) +SU(t::Tuple) = SU{length(t) + 1}(t) # infer N from tuple length + +SymmetryStyle(::SU) = NonAbelianGroup() + +category_label(s::SU) = s.l groupdim(::SU{N}) where {N} = N -trivial(::Type{SU{N}}) where {N} = SU{N}(ntuple(_ -> 0, Val(N))) +trivial(::Type{<:SU{N}}) where {N} = SU{N}(ntuple(_ -> 0, Val(N - 1))) -fundamental(::Type{SU{N}}) where {N} = SU{N}(ntuple(i -> Int(i == 1), Val(N))) +fundamental(::Type{<:SU{N}}) where {N} = SU{N}(ntuple(i -> i == 1, Val(N - 1))) -adjoint(::Type{SU{N}}) where {N} = SU{N}((ntuple(i -> Int(i == 1) + Int(i < N), Val(N)))) +adjoint(::Type{<:SU{N}}) where {N} = SU{N}((ntuple(i -> 1 + (i == 1), Val(N - 1)))) -function dimension(s::SU) +function quantum_dimension(::NonAbelianGroup, s::SU) N = groupdim(s) - l = label(s) + l = (category_label(s)..., 0) d = 1 for k1 in 1:N, k2 in (k1 + 1):N d *= ((k2 - k1) + (l[k1] - l[k2]))//(k2 - k1) @@ -30,23 +42,28 @@ function dimension(s::SU) return Int(d) end -function dual(s::SU) - l = label(s) - nl = ((reverse(cumsum(l[begin:(end - 1)] .- l[(begin + 1):end]))..., 0)) +function GradedAxes.dual(s::SU) + l = category_label(s) + nl = reverse(cumsum((l[begin:(end - 1)] .- l[(begin + 1):end]..., l[end]))) return typeof(s)(nl) end +function Base.show(io::IO, s::SU) + disp = join([string(l) for l in category_label(s)], ", ") + return print(io, "SU(", groupdim(s), ")[", disp, "]") +end + # display SU(N) irrep as a Young tableau with utf8 box char function Base.show(io::IO, ::MIME"text/plain", s::SU) - l = label(s) - if l[1] == 0 # singlet = no box - println(io, "●") - return nothing + if istrivial(s) # singlet = no box + return print(io, "●") end - println("┌─" * "┬─"^(l[1] - 1) * "┐") + N = groupdim(s) + l = category_label(s) + println(io, "┌─" * "┬─"^(l[1] - 1) * "┐") i = 1 - while l[i + 1] != 0 + while i < N - 1 && l[i + 1] != 0 println( io, "├─", @@ -58,26 +75,100 @@ function Base.show(io::IO, ::MIME"text/plain", s::SU) i += 1 end - println(io, "└─", "┴─"^max(0, l[i] - 1), "┘") + print(io, "└─", "┴─"^max(0, l[i] - 1), "┘") return nothing end # # Specializations for the case SU{2} -# Where irreps specified by dimension "d" # -dimension(s::SU{2}) = 1 + label(s)[1] +# optimize implementation +quantum_dimension(s::SU{2}) = category_label(s)[1] + 1 -SU{2}(d::Integer) = SU{2}((d - 1, 0)) +GradedAxes.dual(s::SU{2}) = s -dual(s::SU{2}) = s - -function fusion_rule(s1::SU{2}, s2::SU{2}) - d1, d2 = dimension(s1), dimension(s2) - return [SU{2}(d) for d in (abs(d1 - d2) + 1):2:(d1 + d2 - 1)] +function label_fusion_rule(::Type{<:SU{2}}, s1, s2) + labels = collect((i,) for i in (abs(s1[1] - s2[1])):2:(s1[1] + s2[1])) + degen = ones(Int, length(labels)) + return degen, labels end +# define angular momentum-like interface using half-integers +SU2(h::Number) = SU{2,1}((HalfIntegers.twice(HalfIntegers.HalfInteger(h)),)) + +# display SU2 using half-integers function Base.show(io::IO, s::SU{2}) - return print(io, "SU{2}(", dimension(s), ")") + return print(io, "SU(2)[S=", HalfIntegers.half(quantum_dimension(s) - 1), "]") +end + +function Base.show(io::IO, ::MIME"text/plain", s::SU{2}) + return print(io, "S = ", HalfIntegers.half(quantum_dimension(s) - 1)) +end + +# +# Specializations for the case SU{3} +# aimed for testing non-abelian non self-conjugate representations +# TODO replace with generic implementation +# + +function label_fusion_rule(::Type{<:SU{3}}, left, right) + # Compute SU(3) fusion rules using Littlewood-Richardson rule for Young tableaus. + # See e.g. Di Francesco, Mathieu and Sénéchal, section 13.5.3. + if sum(right) > sum(left) # impose more boxes in left Young tableau + return label_fusion_rule(SU{3}, right, left) + end + + if right[1] == 0 # avoid issues with singlet + return [1], [left] + end + + left_row1 = left[1] + left_row2 = left[2] + right_row1 = right[1] + right_row2 = right[2] + + irreps = [] + + # put a23 boxes on 2nd or 3rd line + a23max1 = 2 * left_row1 # row2a <= row1a + a23max2 = right_row1 # a2 + a3 <= total number of a + a23max = min(a23max1, a23max2) + for a23 in 0:a23max + a3min1 = left_row2 + 2 * a23 - left_row1 - right_row1 + a3min2 = left_row2 - left_row1 + a23 # no a below a: row2a <= row1 + a3min = max(0, a3min1, a3min2) + a3max1 = left_row2 # row3a <= row2a + a3max2 = a23 # a3 <= a2 + a3 + a3max3 = right_row1 - right_row2 # more a than b, right to left: b2 + b3 <= a1 + a2 + a3max = min(a3max1, a3max2, a3max3) + for a3 in a3min:a3max + a2 = a23 - a3 + row1a = left_row1 + right_row1 - a23 + row2a = left_row2 + a23 - a3 + + # cannot put any b on 1st line: row1ab = row1a + b3min1 = row2a + right_row2 - row1a # row2ab <= row1ab = row1a + b3min2 = right_row2 + a23 - right_row1 + b3min = max(0, b3min1, b3min2) + b3max1 = right_row2 # only other.row2 b boxes + b3max2 = (row2a + right_row2 - a3) ÷ 2 # row3ab >= row2ab + b3max3 = right_row1 - a3 # more a than b, right to left: b2 <= a1 + b3max4 = row2a - a3 # no b below b: row2a >= row3ab + b3max = min(b3max1, b3max2, b3max3, b3max4) + for b3 in b3min:b3max + b2 = right_row2 - b3 + row2ab = row2a + b2 + row3ab = a3 + b3 + yt = (row1a - row3ab, row2ab - row3ab) + + push!(irreps, yt) + end + end + end + + unique_labels = sort(unique(irreps)) + degen = [count(==(irr), irreps) for irr in unique_labels] + + return degen, unique_labels end diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/su2.jl b/NDTensors/src/lib/Sectors/src/category_definitions/su2.jl deleted file mode 100644 index 2996b63f08..0000000000 --- a/NDTensors/src/lib/Sectors/src/category_definitions/su2.jl +++ /dev/null @@ -1,22 +0,0 @@ -using HalfIntegers: Half, half, twice - -# -# Conventional SU2 group -# using "J" labels -# - -struct SU2 <: AbstractCategory - j::Half{Int} -end - -dual(s::SU2) = s - -label(s::SU2) = s.j - -trivial(::Type{SU2}) = SU2(0) -fundamental(::Type{SU2}) = SU2(half(1)) -adjoint(::Type{SU2}) = SU2(1) - -dimension(s::SU2) = twice(label(s)) + 1 - -label_fusion_rule(::Type{SU2}, j1, j2) = abs(j1 - j2):(j1 + j2) diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/su2k.jl b/NDTensors/src/lib/Sectors/src/category_definitions/su2k.jl index 86cf5f5abf..736e9987b8 100644 --- a/NDTensors/src/lib/Sectors/src/category_definitions/su2k.jl +++ b/NDTensors/src/lib/Sectors/src/category_definitions/su2k.jl @@ -1,21 +1,23 @@ -using HalfIntegers: Half - # # Quantum 'group' su2ₖ # struct su2{k} <: AbstractCategory - j::Half{Int} + j::HalfIntegers.Half{Int} end +SymmetryStyle(::su2) = NonGroupCategory() + dual(s::su2) = s -label(s::su2) = s.j +category_label(s::su2) = s.j -level(s::su2{k}) where {k} = k +level(::su2{k}) where {k} = k trivial(::Type{su2{k}}) where {k} = su2{k}(0) function label_fusion_rule(::Type{su2{k}}, j1, j2) where {k} - return abs(j1 - j2):min(k - j1 - j2, j1 + j2) + labels = collect(abs(j1 - j2):min(k - j1 - j2, j1 + j2)) + degen = ones(Int, length(labels)) + return degen, labels end diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/u1.jl b/NDTensors/src/lib/Sectors/src/category_definitions/u1.jl index 7af8d22b9f..055b08f630 100644 --- a/NDTensors/src/lib/Sectors/src/category_definitions/u1.jl +++ b/NDTensors/src/lib/Sectors/src/category_definitions/u1.jl @@ -1,19 +1,25 @@ -using HalfIntegers: Half - # # U₁ group (circle group, or particle number, total Sz etc.) # -struct U1 <: AbstractCategory - n::Half{Int} +# Parametric type to allow both integer label as well as +# HalfInteger for easy conversion to/from SU(2) +struct U1{T} <: AbstractCategory + n::T end -dual(u::U1) = U1(-u.n) +SymmetryStyle(::U1) = AbelianGroup() + +GradedAxes.dual(u::U1) = U1(-u.n) -label(u::U1) = u.n +category_label(u::U1) = u.n -dimension(::U1) = 1 +trivial(::Type{U1}) = trivial(U1{Int}) +trivial(::Type{U1{T}}) where {T} = U1(T(0)) -trivial(::Type{U1}) = U1(0) +label_fusion_rule(::Type{<:U1}, n1, n2) = n1 + n2 -label_fusion_rule(::Type{U1}, n1, n2) = (n1 + n2,) +# hide label type in printing +function Base.show(io::IO, u::U1) + return print(io, "U(1)[", category_label(u), "]") +end diff --git a/NDTensors/src/lib/Sectors/src/category_definitions/zn.jl b/NDTensors/src/lib/Sectors/src/category_definitions/zn.jl index 1348052d97..665080b73d 100644 --- a/NDTensors/src/lib/Sectors/src/category_definitions/zn.jl +++ b/NDTensors/src/lib/Sectors/src/category_definitions/zn.jl @@ -1,23 +1,23 @@ -using HalfIntegers: Half - # # Cyclic group Zₙ # struct Z{N} <: AbstractCategory - m::Half{Int} + m::Int Z{N}(m) where {N} = new{N}(m % N) end -label(c::Z) = c.m +SymmetryStyle(::Z) = AbelianGroup() + +category_label(c::Z) = c.m modulus(::Type{Z{N}}) where {N} = N modulus(c::Z) = modulus(typeof(c)) -dimension(::Z) = 1 - trivial(category_type::Type{<:Z}) = category_type(0) -label_fusion_rule(category_type::Type{<:Z}, n1, n2) = ((n1 + n2) % modulus(category_type),) +function label_fusion_rule(category_type::Type{<:Z}, n1, n2) + return (n1 + n2) % modulus(category_type) +end -dual(c::Z) = typeof(c)(mod(-label(c), modulus(c))) +GradedAxes.dual(c::Z) = typeof(c)(mod(-category_label(c), modulus(c))) diff --git a/NDTensors/src/lib/Sectors/src/category_product.jl b/NDTensors/src/lib/Sectors/src/category_product.jl index 12baf35be1..c80c454547 100644 --- a/NDTensors/src/lib/Sectors/src/category_product.jl +++ b/NDTensors/src/lib/Sectors/src/category_product.jl @@ -1,4 +1,7 @@ +# This files defines a structure for Cartesian product of 2 or more fusion categories +# e.g. U(1)×U(1), U(1)×SU2(2)×SU(3) +# ===================================== Definition ======================================= struct CategoryProduct{Categories} <: AbstractCategory cats::Categories global _CategoryProduct(l) = new{typeof(l)}(l) @@ -8,22 +11,30 @@ CategoryProduct(c::CategoryProduct) = _CategoryProduct(categories(c)) categories(s::CategoryProduct) = s.cats -Base.isempty(S::CategoryProduct) = isempty(categories(S)) -Base.length(S::CategoryProduct) = length(categories(S)) -Base.getindex(S::CategoryProduct, args...) = getindex(categories(S), args...) +# ================================= Sectors interface ==================================== +function SymmetryStyle(c::CategoryProduct) + return reduce(combine_styles, map(SymmetryStyle, categories(c)); init=EmptyCategory()) +end + +function quantum_dimension(::NonAbelianGroup, s::CategoryProduct) + return prod(map(quantum_dimension, categories(s))) +end -function fusion_rule(s1::CategoryProduct, s2::CategoryProduct) - return [ - CategoryProduct(l) for l in categories_fusion_rule(categories(s1), categories(s2)) - ] +function quantum_dimension(::NonGroupCategory, s::CategoryProduct) + return prod(map(quantum_dimension, categories(s))) end +GradedAxes.dual(s::CategoryProduct) = CategoryProduct(map(GradedAxes.dual, categories(s))) + +trivial(type::Type{<:CategoryProduct}) = sector(categories_trivial(categories_type(type))) + +# =================================== Base interface ===================================== function Base.:(==)(A::CategoryProduct, B::CategoryProduct) - return categories_equal(categories(A), categories(B)) + return categories_isequal(categories(A), categories(B)) end function Base.show(io::IO, s::CategoryProduct) - (length(s) < 2) && print(io, "sector") + (length(categories(s)) < 2) && print(io, "sector") print(io, "(") symbol = "" for p in pairs(categories(s)) @@ -34,27 +45,145 @@ function Base.show(io::IO, s::CategoryProduct) return print(io, ")") end -category_show(io::IO, k, v) = print(io, v) - +category_show(io::IO, ::Int, v) = print(io, v) category_show(io::IO, k::Symbol, v) = print(io, "($k=$v,)") +function Base.isless(s1::C, s2::C) where {C<:CategoryProduct} + return categories_isless(categories(s1), categories(s2)) +end +function Base.isless(s1::CategoryProduct, s2::CategoryProduct) + return categories_isless(categories(s1), categories(s2)) +end + +# ======================================= shared ========================================= +# there are 2 implementations for CategoryProduct +# - ordered-like with a Tuple +# - dictionary-like with a NamedTuple + +function categories_fusion_rule(cats1, cats2) + diff_cat = CategoryProduct(find_diff(cats1, cats2)) + shared1, shared2 = find_common(cats1, cats2) + fused = map(fusion_rule, values(shared1), values(shared2)) + return recover_key(typeof(shared1), fused) × diff_cat +end + +# get clean results when mixing implementations +categories_isequal(::Tuple, ::NamedTuple) = false +categories_isequal(::NamedTuple, ::Tuple) = false + +categories_isless(::NamedTuple, ::Tuple) = throw(ArgumentError("Not implemented")) +categories_isless(::Tuple, ::NamedTuple) = throw(ArgumentError("Not implemented")) + +categories_type(::Type{<:CategoryProduct{T}}) where {T} = T + +recover_key(T::Type, t::Tuple{Vararg{AbstractCategory}}) = sector(T, t) +recover_key(T::Type, c::AbstractCategory) = recover_key(T, (c,)) +recover_key(T::Type, c::CategoryProduct) = recover_key(T, categories(c)) +function recover_key(T::Type, fused::Tuple) + # here fused contains at leat one GradedUnitRange + g0 = reduce(×, fused) + # convention: keep unsorted blocklabels as produced by F order loops in × + new_labels = recover_key.(T, GradedAxes.blocklabels(g0)) + new_blocklengths = + LabelledNumbers.labelled.(GradedAxes.unlabel.(BlockArrays.blocklengths(g0)), new_labels) + return GradedAxes.gradedrange(new_blocklengths) +end + +sector(T::Type{<:CategoryProduct}, cats::Tuple) = sector(categories_type(T), cats) +sector(T::Type, cats::Tuple) = sector(T(cats)) # recover NamedTuple +sector(args...; kws...) = CategoryProduct(args...; kws...) + +# ================================= Cartesian Product ==================================== ×(c1::AbstractCategory, c2::AbstractCategory) = ×(CategoryProduct(c1), CategoryProduct(c2)) function ×(p1::CategoryProduct, p2::CategoryProduct) return CategoryProduct(categories_product(categories(p1), categories(p2))) end -categories_product(l1::NamedTuple, l2::NamedTuple) = union_keys(l1, l2) - -categories_product(l1::Tuple, l2::Tuple) = (l1..., l2...) - +×(a, g::AbstractUnitRange) = ×(to_graded_axis(a), g) +×(g::AbstractUnitRange, b) = ×(g, to_graded_axis(b)) ×(nt1::NamedTuple, nt2::NamedTuple) = ×(CategoryProduct(nt1), CategoryProduct(nt2)) ×(c1::NamedTuple, c2::AbstractCategory) = ×(CategoryProduct(c1), CategoryProduct(c2)) ×(c1::AbstractCategory, c2::NamedTuple) = ×(CategoryProduct(c1), CategoryProduct(c2)) -# -# Dictionary-like implementation -# +function ×(l1::LabelledNumbers.LabelledInteger, l2::LabelledNumbers.LabelledInteger) + c3 = LabelledNumbers.label(l1) × LabelledNumbers.label(l2) + m3 = LabelledNumbers.unlabel(l1) * LabelledNumbers.unlabel(l2) + return LabelledNumbers.labelled(m3, c3) +end + +function ×(g1::AbstractUnitRange, g2::AbstractUnitRange) + v = map( + ((l1, l2),) -> l1 × l2, + Iterators.flatten(( + Iterators.product(BlockArrays.blocklengths(g1), BlockArrays.blocklengths(g2)), + ),), + ) + return GradedAxes.gradedrange(v) +end + +# ==================================== Fusion rules ====================================== +# generic case: fusion returns a GradedAxes, even for fusion with Empty +function fusion_rule(::SymmetryStyle, s1::CategoryProduct, s2::CategoryProduct) + return to_graded_axis(categories_fusion_rule(categories(s1), categories(s2))) +end + +# Abelian case: fusion returns CategoryProduct +function fusion_rule(::AbelianGroup, s1::CategoryProduct, s2::CategoryProduct) + return categories_fusion_rule(categories(s1), categories(s2)) +end + +# Empty case +function fusion_rule( + ::EmptyCategory, ::CategoryProduct{Tuple{}}, ::CategoryProduct{Tuple{}} +) + return sector() +end + +# EmptyCategory acts as trivial on any AbstractCategory, not just CategoryProduct +function fusion_rule(::SymmetryStyle, ::CategoryProduct{Tuple{}}, c::AbstractCategory) + return to_graded_axis(c) +end +function fusion_rule(::SymmetryStyle, ::CategoryProduct{Tuple{}}, c::CategoryProduct) + return to_graded_axis(c) +end +function fusion_rule(::SymmetryStyle, c::AbstractCategory, ::CategoryProduct{Tuple{}}) + return to_graded_axis(c) +end +function fusion_rule(::SymmetryStyle, c::CategoryProduct, ::CategoryProduct{Tuple{}}) + return to_graded_axis(c) +end + +# abelian case: return Category +fusion_rule(::AbelianGroup, ::CategoryProduct{Tuple{}}, c::AbstractCategory) = c +fusion_rule(::AbelianGroup, ::CategoryProduct{Tuple{}}, c::CategoryProduct) = c +fusion_rule(::AbelianGroup, c::AbstractCategory, ::CategoryProduct{Tuple{}}) = c +fusion_rule(::AbelianGroup, c::CategoryProduct, ::CategoryProduct{Tuple{}}) = c + +# =============================== Ordered implementation ================================= +CategoryProduct(t::Tuple) = _CategoryProduct(t) +CategoryProduct(cats::AbstractCategory...) = CategoryProduct((cats...,)) + +categories_isequal(o1::Tuple, o2::Tuple) = (o1 == o2) + +categories_isless(::Tuple, ::Tuple) = throw(ArgumentError("Not implemented")) +categories_isless(t1::T, t2::T) where {T<:Tuple} = t1 < t2 +categories_product(l1::Tuple, l2::Tuple) = (l1..., l2...) + +categories_trivial(type::Type{<:Tuple}) = trivial.(fieldtypes(type)) + +function find_common(t1::Tuple, t2::Tuple) + n = min(length(t1), length(t2)) + return t1[begin:n], t2[begin:n] +end + +function find_diff(t1::Tuple, t2::Tuple) + n1 = length(t1) + n2 = length(t2) + return n1 < n2 ? t2[(n1 + 1):end] : t1[(n2 + 1):end] +end + +# =========================== Dictionary-like implementation ============================= function CategoryProduct(nt::NamedTuple) categories = sort_keys(nt) return _CategoryProduct(categories) @@ -62,50 +191,60 @@ end CategoryProduct(; kws...) = CategoryProduct((; kws...)) +# avoid having 2 different kinds of EmptyCategory: cast empty NamedTuple to Tuple{} +CategoryProduct(::NamedTuple{()}) = CategoryProduct(()) + function CategoryProduct(pairs::Pair...) keys = ntuple(n -> Symbol(pairs[n][1]), length(pairs)) vals = ntuple(n -> pairs[n][2], length(pairs)) return CategoryProduct(NamedTuple{keys}(vals)) end -function categories_fusion_rule(A::NamedTuple, B::NamedTuple) - qs = [A] - for (la, lb) in zip(pairs(intersect_keys(A, B)), pairs(intersect_keys(B, A))) - @assert la[1] == lb[1] - fused_vals = ⊗(la[2], lb[2]) - qs = [union_keys((; la[1] => v), q) for v in fused_vals for q in qs] - end - # Include sectors of B not in A - qs = [union_keys(q, B) for q in qs] - return qs +# sector() acts as empty NamedTuple +function categories_isequal(nt::NamedTuple, ::Tuple{}) + return categories_isequal(nt, categories_trivial(typeof(nt))) end - -function categories_equal(A::NamedTuple, B::NamedTuple) - common_categories = zip(pairs(intersect_keys(A, B)), pairs(intersect_keys(B, A))) - common_categories_match = all(nl -> (nl[1] == nl[2]), common_categories) - unique_categories_zero = all(l -> istrivial(l), symdiff_keys(A, B)) - return common_categories_match && unique_categories_zero +function categories_isequal(::Tuple{}, nt::NamedTuple) + return categories_isequal(categories_trivial(typeof(nt)), nt) +end +function categories_isequal(nt1::NamedTuple, nt2::NamedTuple) + return ==(sym_categories_insert_unspecified(nt1, nt2)...) end -# -# Ordered implementation -# +function categories_isless(nt::NamedTuple, ::Tuple{}) + return categories_isless(nt, categories_trivial(typeof(nt))) +end +function categories_isless(::Tuple{}, nt::NamedTuple) + return categories_isless(categories_trivial(typeof(nt)), nt) +end +function categories_isless(nt1::NamedTuple, nt2::NamedTuple) + return isless(sym_categories_insert_unspecified(nt1, nt2)...) +end -CategoryProduct(t::Tuple) = _CategoryProduct(t) -CategoryProduct(cats::AbstractCategory...) = CategoryProduct((cats...,)) +function sym_categories_insert_unspecified(nt1::NamedTuple, nt2::NamedTuple) + return categories_insert_unspecified(nt1, nt2), categories_insert_unspecified(nt2, nt1) +end -categories_equal(o1::Tuple, o2::Tuple) = (o1 == o2) +function categories_insert_unspecified(nt1::NamedTuple, nt2::NamedTuple) + diff1 = categories_trivial(typeof(setdiff_keys(nt2, nt1))) + return sort_keys(union_keys(nt1, diff1)) +end -function categories_fusion_rule(o1::Tuple, o2::Tuple) - N = length(o1) - length(o2) == N || - throw(DimensionMismatch("Ordered CategoryProduct must have same size in ⊗")) - os = [o1] - replace(o, n, val) = ntuple(m -> (m == n) ? val : o[m], length(o)) - for n in 1:N - os = [replace(o, n, f) for f in ⊗(o1[n], o2[n]) for o in os] +categories_product(l1::NamedTuple, ::Tuple{}) = l1 +categories_product(::Tuple{}, l2::NamedTuple) = l2 +function categories_product(l1::NamedTuple, l2::NamedTuple) + if length(intersect_keys(l1, l2)) > 0 + throw(ArgumentError("Cannot define product of shared keys")) end - return os + return union_keys(l1, l2) end -sector(args...; kws...) = CategoryProduct(args...; kws...) +function categories_trivial(type::Type{<:NamedTuple{Keys}}) where {Keys} + return NamedTuple{Keys}(trivial.(fieldtypes(type))) +end + +function find_common(nt1::NamedTuple, nt2::NamedTuple) + return intersect_keys(nt1, nt2), intersect_keys(nt2, nt1) +end + +find_diff(nt1::NamedTuple, nt2::NamedTuple) = symdiff_keys(nt1, nt2) diff --git a/NDTensors/src/lib/Sectors/src/symmetry_style.jl b/NDTensors/src/lib/Sectors/src/symmetry_style.jl new file mode 100644 index 0000000000..f65a65fc91 --- /dev/null +++ b/NDTensors/src/lib/Sectors/src/symmetry_style.jl @@ -0,0 +1,31 @@ +# This file defines SymmetryStyle, a trait to distinguish abelian groups, non-abelian groups +# and non-group fusion categories. + +using BlockArrays + +using NDTensors.LabelledNumbers +using NDTensors.GradedAxes + +abstract type SymmetryStyle end + +struct AbelianGroup <: SymmetryStyle end +struct NonAbelianGroup <: SymmetryStyle end +struct NonGroupCategory <: SymmetryStyle end +struct EmptyCategory <: SymmetryStyle end # CategoryProduct with zero category inside + +combine_styles(::AbelianGroup, ::AbelianGroup) = AbelianGroup() +combine_styles(::AbelianGroup, ::NonAbelianGroup) = NonAbelianGroup() +combine_styles(::AbelianGroup, ::NonGroupCategory) = NonGroupCategory() +combine_styles(::NonAbelianGroup, ::AbelianGroup) = NonAbelianGroup() +combine_styles(::NonAbelianGroup, ::NonAbelianGroup) = NonAbelianGroup() +combine_styles(::NonAbelianGroup, ::NonGroupCategory) = NonGroupCategory() +combine_styles(::NonGroupCategory, ::SymmetryStyle) = NonGroupCategory() +combine_styles(::NonGroupCategory, ::EmptyCategory) = NonGroupCategory() +combine_styles(::EmptyCategory, s::SymmetryStyle) = s +combine_styles(s::SymmetryStyle, ::EmptyCategory) = s +combine_styles(::EmptyCategory, ::EmptyCategory) = EmptyCategory() + +SymmetryStyle(l::LabelledNumbers.LabelledInteger) = SymmetryStyle(LabelledNumbers.label(l)) + +# crash for empty g. Currently impossible to construct. +SymmetryStyle(g::AbstractUnitRange) = SymmetryStyle(first(g)) diff --git a/NDTensors/src/lib/Sectors/test/test_category.jl b/NDTensors/src/lib/Sectors/test/test_category.jl deleted file mode 100644 index e14582afc8..0000000000 --- a/NDTensors/src/lib/Sectors/test/test_category.jl +++ /dev/null @@ -1,192 +0,0 @@ -@eval module $(gensym()) -using NDTensors.Sectors: - Fib, - Ising, - SU, - SU2, - U1, - Z, - ⊗, - ⊕, - dimension, - dual, - istrivial, - trivial, - fundamental, - adjoint -using Test: @test, @testset -@testset "Test Category Types" begin - @testset "U(1)" begin - q1 = U1(1) - q2 = U1(2) - q3 = U1(3) - - @test dimension(q1) == 1 - @test dimension(q2) == 1 - - @test q1 ⊗ q1 == [q2] - @test q1 ⊗ q2 == [q3] - @test q2 ⊗ q1 == [q3] - - @test trivial(U1) == U1(0) - @test istrivial(U1(0)) - - @test dual(U1(2)) == U1(-2) - @test isless(U1(1), U1(2)) - @test !isless(U1(2), U1(1)) - end - - @testset "Z₂" begin - z0 = Z{2}(0) - z1 = Z{2}(1) - - @test trivial(Z{2}) == Z{2}(0) - @test istrivial(Z{2}(0)) - - @test dimension(z0) == 1 - @test dimension(z1) == 1 - - @test dual(z0) == z0 - @test dual(z1) == z1 - - @test z0 ⊗ z0 == [z0] - @test z0 ⊗ z1 == [z1] - @test z1 ⊗ z1 == [z0] - - @test dual(Z{2}(1)) == Z{2}(1) - @test isless(Z{2}(0), Z{2}(1)) - @test !isless(Z{2}(1), Z{2}(0)) - end - - @testset "SU2" begin - j1 = SU2(0) - j2 = SU2(1//2) - j3 = SU2(1) - j4 = SU2(3//2) - j5 = SU2(2) - - @test trivial(SU2) == SU2(0) - @test istrivial(SU2(0)) - - @test fundamental(SU2) == SU2(1//2) - @test adjoint(SU2) == SU2(1) - - @test dimension(j1) == 1 - @test dimension(j2) == 2 - @test dimension(j3) == 3 - @test dimension(j4) == 4 - - @test dual(j1) == j1 - @test dual(j2) == j2 - @test dual(j3) == j3 - @test dual(j4) == j4 - - @test j1 ⊗ j2 == [j2] - @test j2 ⊗ j2 == j1 ⊕ j3 - @test j2 ⊗ j3 == j2 ⊕ j4 - @test j3 ⊗ j3 == j1 ⊕ j3 ⊕ j5 - end - - @testset "SU(2)" begin - j1 = SU{2}(1) - j2 = SU{2}(2) - j3 = SU{2}(3) - j4 = SU{2}(4) - j5 = SU{2}(5) - - @test trivial(SU{2}) == SU{2}(1) - @test istrivial(SU{2}(1)) - - @test fundamental(SU{2}) == SU{2}(2) - @test adjoint(SU{2}) == SU{2}(3) - - @test dimension(j1) == 1 - @test dimension(j2) == 2 - @test dimension(j3) == 3 - @test dimension(j4) == 4 - - @test dual(j1) == j1 - @test dual(j2) == j2 - @test dual(j3) == j3 - @test dual(j4) == j4 - - @test j1 ⊗ j2 == [j2] - @test j2 ⊗ j2 == j1 ⊕ j3 - @test j2 ⊗ j3 == j2 ⊕ j4 - @test j3 ⊗ j3 == j1 ⊕ j3 ⊕ j5 - end - - @testset "SU(N)" begin - f3 = SU{3}((1, 0, 0)) - f4 = SU{4}((1, 0, 0, 0)) - ad3 = SU{3}((2, 1, 0)) - ad4 = SU{4}((2, 1, 1, 0)) - - @test trivial(SU{3}) == SU{3}((0, 0, 0)) - @test istrivial(SU{3}((0, 0, 0))) - @test trivial(SU{4}) == SU{4}((0, 0, 0, 0)) - @test istrivial(SU{4}((0, 0, 0, 0))) - - @test fundamental(SU{3}) == f3 - @test adjoint(SU{3}) == ad3 - @test fundamental(SU{4}) == f4 - @test adjoint(SU{4}) == ad4 - - @test dual(f3) == SU{3}((1, 1, 0)) - @test dual(f4) == SU{4}((1, 1, 1, 0)) - @test dual(ad3) == ad3 - @test dual(ad4) == ad4 - - @test dimension(f3) == 3 - @test dimension(f4) == 4 - @test dimension(ad3) == 8 - @test dimension(ad4) == 15 - @test dimension(SU{3}((4, 2, 0))) == 27 - @test dimension(SU{3}((3, 3, 0))) == 10 - @test dimension(SU{3}((3, 0, 0))) == 10 - @test dimension(SU{3}((0, 0, 0))) == 1 - end - - @testset "Fibonacci" begin - ı = Fib("1") - τ = Fib("τ") - - @test trivial(Fib) == ı - @test istrivial(ı) - - @test dual(ı) == ı - @test dual(τ) == τ - - @test dimension(ı) == 1 - @test dimension(τ) == ((1 + √5) / 2) - - @test ı ⊗ ı == [ı] - @test ı ⊗ τ == [τ] - @test τ ⊗ ı == [τ] - @test τ ⊗ τ == ı ⊕ τ - end - - @testset "Ising" begin - ı = Ising("1") - σ = Ising("σ") - ψ = Ising("ψ") - - @test trivial(Ising) == ı - @test istrivial(ı) - - @test dual(ı) == ı - @test dual(σ) == σ - @test dual(ψ) == ψ - - @test ı ⊗ ı == [ı] - @test ı ⊗ σ == [σ] - @test σ ⊗ ı == [σ] - @test ı ⊗ ψ == [ψ] - @test ψ ⊗ ı == [ψ] - @test σ ⊗ σ == ı ⊕ ψ - @test σ ⊗ ψ == [σ] - @test ψ ⊗ σ == [σ] - @test ψ ⊗ ψ == [ı] - end -end -end diff --git a/NDTensors/src/lib/Sectors/test/test_category_product.jl b/NDTensors/src/lib/Sectors/test/test_category_product.jl index 58f095c110..8f70741730 100644 --- a/NDTensors/src/lib/Sectors/test/test_category_product.jl +++ b/NDTensors/src/lib/Sectors/test/test_category_product.jl @@ -1,135 +1,582 @@ @eval module $(gensym()) -using NDTensors.Sectors: Fib, Ising, SU, SU2, U1, Z, ⊗, ⊕, ×, sector -using Test: @test, @testset, @test_throws +using NDTensors.Sectors: + ×, + ⊗, + Fib, + Ising, + SU, + SU2, + U1, + Z, + block_dimensions, + categories, + quantum_dimension, + sector, + trivial +using NDTensors.GradedAxes: dual, fusion_product, gradedisequal, gradedrange +using Test: @inferred, @test, @testset, @test_throws + +macro inferred_latest(ex) + if VERSION < v"1.10" + return esc(:($ex)) + end + return esc(:(@inferred $ex)) +end + +@testset "Test Ordered Products" begin + @testset "Ordered Constructor" begin + s = sector(U1(1)) + @test length(categories(s)) == 1 + @test (@inferred quantum_dimension(s)) == 1 + @test (@inferred dual(s)) == sector(U1(-1)) + @test categories(s)[1] == U1(1) + @test (@inferred_latest trivial(s)) == sector(U1(0)) + + s = sector(U1(1), U1(2)) + @test length(categories(s)) == 2 + @test (@inferred quantum_dimension(s)) == 1 + @test (@inferred dual(s)) == sector(U1(-1), U1(-2)) + @test categories(s)[1] == U1(1) + @test categories(s)[2] == U1(2) + @test (@inferred_latest trivial(s)) == sector(U1(0), U1(0)) + + s = U1(1) × SU2(1//2) × U1(3) + @test length(categories(s)) == 3 + @test (@inferred quantum_dimension(s)) == 2 + @test (@inferred dual(s)) == U1(-1) × SU2(1//2) × U1(-3) + @test categories(s)[1] == U1(1) + @test categories(s)[2] == SU2(1//2) + @test categories(s)[3] == U1(3) + @test (@inferred_latest trivial(s)) == sector(U1(0), SU2(0), U1(0)) + @test (@inferred sector(typeof(categories(s)), categories(s))) == s + @test (@inferred sector(typeof(s), categories(s))) == s + + s = U1(3) × SU2(1//2) × Fib("τ") + @test length(categories(s)) == 3 + @test (@inferred quantum_dimension(s)) == 1.0 + √5 + @test dual(s) == U1(-3) × SU2(1//2) × Fib("τ") + @test categories(s)[1] == U1(3) + @test categories(s)[2] == SU2(1//2) + @test categories(s)[3] == Fib("τ") + @test (@inferred_latest trivial(s)) == sector(U1(0), SU2(0), Fib("1")) + end + + @testset "Ordered comparisons" begin + # convention: categories must have same length to evaluate as equal + @test sector(U1(1), SU2(1)) == sector(U1(1), SU2(1)) + @test sector(U1(1), SU2(0)) != sector(U1(1), SU2(1)) + @test sector(U1(0), SU2(1)) != sector(U1(1), SU2(1)) + @test sector(U1(1)) != sector(U1(1), U1(0)) + + # convention: categories must have same length to be compared + @test sector(U1(0)) < sector((U1(1))) + @test sector(U1(0), U1(2)) < sector((U1(1)), U1(0)) + @test_throws ArgumentError sector(U1(0)) < sector(U1(1), U1(2)) + end + + @testset "Quantum dimension and GradedUnitRange" begin + g = gradedrange([(U1(0) × Z{2}(0)) => 1, (U1(1) × Z{2}(0)) => 2]) # abelian + @test (@inferred quantum_dimension(g)) == 3 + + g = gradedrange([ # non-abelian + (SU2(0) × SU2(0)) => 1, + (SU2(1) × SU2(0)) => 1, + (SU2(0) × SU2(1)) => 1, + (SU2(1) × SU2(1)) => 1, + ]) + @test (@inferred quantum_dimension(g)) == 16 + @test (@inferred block_dimensions(g)) == [1, 3, 3, 9] + + # mixed group + g = gradedrange([(U1(2) × SU2(0) × Z{2}(0)) => 1, (U1(2) × SU2(1) × Z{2}(0)) => 1]) + @test (@inferred quantum_dimension(g)) == 4 + @test (@inferred block_dimensions(g)) == [1, 3] + g = gradedrange([(SU2(0) × U1(0) × SU2(1//2)) => 1, (SU2(0) × U1(1) × SU2(1//2)) => 1]) + @test (@inferred quantum_dimension(g)) == 4 + @test (@inferred block_dimensions(g)) == [2, 2] + + # NonGroupCategory + g_fib = gradedrange([(Fib("1") × Fib("1")) => 1]) + g_ising = gradedrange([(Ising("1") × Ising("1")) => 1]) + @test (@inferred quantum_dimension((Fib("1") × Fib("1")))) == 1.0 + @test (@inferred quantum_dimension(g_fib)) == 1.0 + @test (@inferred quantum_dimension(g_ising)) == 1.0 + @test (@inferred quantum_dimension((Ising("1") × Ising("1")))) == 1.0 + @test (@inferred block_dimensions(g_fib)) == [1.0] + @test (@inferred block_dimensions(g_ising)) == [1.0] + + @test (@inferred quantum_dimension(U1(1) × Fib("1"))) == 1.0 + @test (@inferred quantum_dimension(gradedrange([U1(1) × Fib("1") => 1]))) == 1.0 + + # mixed product Abelian / NonAbelian / NonGroup + g = gradedrange([ + (U1(2) × SU2(0) × Ising(1)) => 1, + (U1(2) × SU2(1) × Ising(1)) => 1, + (U1(2) × SU2(0) × Ising("ψ")) => 1, + (U1(2) × SU2(1) × Ising("ψ")) => 1, + ]) + @test (@inferred quantum_dimension(g)) == 8.0 + @test (@inferred block_dimensions(g)) == [1.0, 3.0, 1.0, 3.0] + + ϕ = (1 + √5) / 2 + g = gradedrange([ + (Fib("1") × SU2(0) × U1(2)) => 1, + (Fib("1") × SU2(1) × U1(2)) => 1, + (Fib("τ") × SU2(0) × U1(2)) => 1, + (Fib("τ") × SU2(1) × U1(2)) => 1, + ]) + @test (@inferred quantum_dimension(g)) == 4.0 + 4.0ϕ + @test (@inferred block_dimensions(g)) == [1.0, 3.0, 1.0ϕ, 3.0ϕ] + end + + @testset "Fusion of Abelian products" begin + p1 = sector(U1(1)) + p2 = sector(U1(2)) + @test (@inferred p1 ⊗ p2) == sector(U1(3)) + + p11 = U1(1) × U1(1) + @test (@inferred p11 ⊗ p11) == U1(2) × U1(2) + + p123 = U1(1) × U1(2) × U1(3) + @test (@inferred p123 ⊗ p123) == U1(2) × U1(4) × U1(6) + + s1 = sector(U1(1), Z{2}(1)) + s2 = sector(U1(0), Z{2}(0)) + @test (@inferred s1 ⊗ s2) == U1(1) × Z{2}(1) + end + + @testset "Fusion of NonAbelian products" begin + p0 = sector(SU2(0)) + ph = sector(SU2(1//2)) + @test gradedisequal((@inferred p0 ⊗ ph), gradedrange([sector(SU2(1//2)) => 1])) + + phh = SU2(1//2) × SU2(1//2) + @test gradedisequal( + phh ⊗ phh, + gradedrange([ + (SU2(0) × SU2(0)) => 1, + (SU2(1) × SU2(0)) => 1, + (SU2(0) × SU2(1)) => 1, + (SU2(1) × SU2(1)) => 1, + ]), + ) + @test gradedisequal( + (@inferred phh ⊗ phh), + gradedrange([ + (SU2(0) × SU2(0)) => 1, + (SU2(1) × SU2(0)) => 1, + (SU2(0) × SU2(1)) => 1, + (SU2(1) × SU2(1)) => 1, + ]), + ) + end + + @testset "Fusion of NonGroupCategory products" begin + ı = Fib("1") + τ = Fib("τ") + s = ı × ı + @test gradedisequal((@inferred s ⊗ s), gradedrange([s => 1])) + + s = τ × τ + @test gradedisequal( + (@inferred s ⊗ s), + gradedrange([(ı × ı) => 1, (τ × ı) => 1, (ı × τ) => 1, (τ × τ) => 1]), + ) + + σ = Ising("σ") + ψ = Ising("ψ") + s = τ × σ + g = gradedrange([ + (ı × Ising("1")) => 1, (τ × Ising("1")) => 1, (ı × ψ) => 1, (τ × ψ) => 1 + ]) + @test gradedisequal((@inferred s ⊗ s), g) + end + + @testset "Fusion of mixed Abelian and NonAbelian products" begin + p2h = U1(2) × SU2(1//2) + p1h = U1(1) × SU2(1//2) + @test gradedisequal( + (@inferred p2h ⊗ p1h), gradedrange([(U1(3) × SU2(0)) => 1, (U1(3) × SU2(1)) => 1]) + ) + + p1h1 = U1(1) × SU2(1//2) × Z{2}(1) + @test gradedisequal( + (@inferred p1h1 ⊗ p1h1), + gradedrange([(U1(2) × SU2(0) × Z{2}(0)) => 1, (U1(2) × SU2(1) × Z{2}(0)) => 1]), + ) + end + + @testset "Fusion of fully mixed products" begin + s = U1(1) × SU2(1//2) × Ising("σ") + @test gradedisequal( + (@inferred s ⊗ s), + gradedrange([ + (U1(2) × SU2(0) × Ising("1")) => 1, + (U1(2) × SU2(1) × Ising("1")) => 1, + (U1(2) × SU2(0) × Ising("ψ")) => 1, + (U1(2) × SU2(1) × Ising("ψ")) => 1, + ]), + ) + + ı = Fib("1") + τ = Fib("τ") + s = SU2(1//2) × U1(1) × τ + @test gradedisequal( + (@inferred s ⊗ s), + gradedrange([ + (SU2(0) × U1(2) × ı) => 1, + (SU2(1) × U1(2) × ı) => 1, + (SU2(0) × U1(2) × τ) => 1, + (SU2(1) × U1(2) × τ) => 1, + ]), + ) + + s = U1(1) × ı × τ + @test gradedisequal( + (@inferred s ⊗ s), gradedrange([(U1(2) × ı × ı) => 1, (U1(2) × ı × τ) => 1]) + ) + end + + @testset "Fusion of different length Categories" begin + @test sector(U1(1) × U1(0)) ⊗ sector(U1(1)) == sector(U1(2) × U1(0)) + @test gradedisequal( + (@inferred sector(SU2(0) × SU2(0)) ⊗ sector(SU2(1))), + gradedrange([sector(SU2(1) × SU2(0)) => 1]), + ) + + @test gradedisequal( + (@inferred sector(SU2(1) × U1(1)) ⊗ sector(SU2(0))), + gradedrange([sector(SU2(1) × U1(1)) => 1]), + ) + @test gradedisequal( + (@inferred sector(U1(1) × SU2(1)) ⊗ sector(U1(2))), + gradedrange([sector(U1(3) × SU2(1)) => 1]), + ) + + # check incompatible categories + p12 = Z{2}(1) × U1(2) + z12 = Z{2}(1) × Z{2}(1) + @test_throws MethodError p12 ⊗ z12 + end + + @testset "GradedUnitRange fusion rules" begin + s1 = U1(1) × SU2(1//2) × Ising("σ") + s2 = U1(0) × SU2(1//2) × Ising("1") + g1 = gradedrange([s1 => 2]) + g2 = gradedrange([s2 => 1]) + @test gradedisequal( + (@inferred fusion_product(g1, g2)), + gradedrange([U1(1) × SU2(0) × Ising("σ") => 2, U1(1) × SU2(1) × Ising("σ") => 2]), + ) + end +end + @testset "Test Named Category Products" begin @testset "Construct from × of NamedTuples" begin + s = (A=U1(1),) × (B=Z{2}(0),) + @test length(categories(s)) == 2 + @test categories(s)[:A] == U1(1) + @test categories(s)[:B] == Z{2}(0) + @test (@inferred quantum_dimension(s)) == 1 + @test (@inferred dual(s)) == (A=U1(-1),) × (B=Z{2}(0),) + @test (@inferred_latest trivial(s)) == (A=U1(0),) × (B=Z{2}(0),) + s = (A=U1(1),) × (B=SU2(2),) - @test length(s) == 2 - @test s[:A] == U1(1) - @test s[:B] == SU2(2) + @test length(categories(s)) == 2 + @test categories(s)[:A] == U1(1) + @test categories(s)[:B] == SU2(2) + @test (@inferred quantum_dimension(s)) == 5 + @test (@inferred dual(s)) == (A=U1(-1),) × (B=SU2(2),) + @test (@inferred_latest trivial(s)) == (A=U1(0),) × (B=SU2(0),) + @test (@inferred sector(typeof(categories(s)), Tuple(categories(s)))) == s + @test (@inferred sector(typeof(s), Tuple(categories(s)))) == s s = s × (C=Ising("ψ"),) - @test length(s) == 3 - @test s[:C] == Ising("ψ") + @test length(categories(s)) == 3 + @test categories(s)[:C] == Ising("ψ") + @test (@inferred_latest quantum_dimension(s)) == 5.0 + @test (@inferred dual(s)) == (A=U1(-1),) × (B=SU2(2),) × (C=Ising("ψ"),) + + s1 = (A=U1(1),) × (B=Z{2}(0),) + s2 = (A=U1(1),) × (C=Z{2}(0),) + @test_throws ArgumentError s1 × s2 end @testset "Construct from Pairs" begin s = sector("A" => U1(2)) - @test length(s) == 1 - @test s[:A] == U1(2) + @test length(categories(s)) == 1 + @test categories(s)[:A] == U1(2) @test s == sector(; A=U1(2)) + @test (@inferred quantum_dimension(s)) == 1 + @test (@inferred dual(s)) == sector("A" => U1(-2)) + @test (@inferred_latest trivial(s)) == sector(; A=U1(0)) s = sector("B" => Ising("ψ"), :C => Z{2}(1)) - @test length(s) == 2 - @test s[:B] == Ising("ψ") - @test s[:C] == Z{2}(1) + @test length(categories(s)) == 2 + @test categories(s)[:B] == Ising("ψ") + @test categories(s)[:C] == Z{2}(1) + @test (@inferred quantum_dimension(s)) == 1.0 + end + + @testset "Comparisons with unspecified labels" begin + # convention: categories evaluate as equal if unmatched labels are trivial + # this is different from ordered tuple convention + q2 = sector(; N=U1(2)) + q20 = (N=U1(2),) × (J=SU2(0),) + @test q20 == q2 + @test !(q20 < q2) + @test !(q2 < q20) + + q21 = (N=U1(2),) × (J=SU2(1),) + @test q21 != q2 + @test q20 < q21 + @test q2 < q21 + + a = (A=U1(0),) × (B=U1(2),) + b = (B=U1(2),) × (C=U1(0),) + @test a == b + c = (B=U1(2),) × (C=U1(1),) + @test a != c + end + + @testset "Quantum dimension and GradedUnitRange" begin + g = gradedrange([sector(; A=U1(0), B=Z{2}(0)) => 1, sector(; A=U1(1), B=Z{2}(0)) => 2]) # abelian + @test (@inferred quantum_dimension(g)) == 3 + + g = gradedrange([ # non-abelian + sector(; A=SU2(0), B=SU2(0)) => 1, + sector(; A=SU2(1), B=SU2(0)) => 1, + sector(; A=SU2(0), B=SU2(1)) => 1, + sector(; A=SU2(1), B=SU2(1)) => 1, + ]) + @test (@inferred quantum_dimension(g)) == 16 + + # mixed group + g = gradedrange([ + sector(; A=U1(2), B=SU2(0), C=Z{2}(0)) => 1, + sector(; A=U1(2), B=SU2(1), C=Z{2}(0)) => 1, + ]) + @test (@inferred quantum_dimension(g)) == 4 + g = gradedrange([ + sector(; A=SU2(0), B=Z{2}(0), C=SU2(1//2)) => 1, + sector(; A=SU2(0), B=Z{2}(1), C=SU2(1//2)) => 1, + ]) + @test (@inferred quantum_dimension(g)) == 4 + + # non group categories + g_fib = gradedrange([sector(; A=Fib("1"), B=Fib("1")) => 1]) + g_ising = gradedrange([sector(; A=Ising("1"), B=Ising("1")) => 1]) + @test (@inferred quantum_dimension(g_fib)) == 1.0 + @test (@inferred quantum_dimension(g_ising)) == 1.0 + + # mixed product Abelian / NonAbelian / NonGroup + g = gradedrange([ + sector(; A=U1(2), B=SU2(0), C=Ising(1)) => 1, + sector(; A=U1(2), B=SU2(1), C=Ising(1)) => 1, + sector(; A=U1(2), B=SU2(0), C=Ising("ψ")) => 1, + sector(; A=U1(2), B=SU2(1), C=Ising("ψ")) => 1, + ]) + @test (@inferred_latest quantum_dimension(g)) == 8.0 + + g = gradedrange([ + sector(; A=Fib("1"), B=SU2(0), C=U1(2)) => 1, + sector(; A=Fib("1"), B=SU2(1), C=U1(2)) => 1, + sector(; A=Fib("τ"), B=SU2(0), C=U1(2)) => 1, + sector(; A=Fib("τ"), B=SU2(1), C=U1(2)) => 1, + ]) + @test (@inferred quantum_dimension(g)) == 4.0 + 4.0quantum_dimension(Fib("τ")) end - @testset "Multiple U(1)'s" begin + @testset "Fusion of Abelian products" begin q00 = sector() q10 = sector(; A=U1(1)) q01 = sector(; B=U1(1)) q11 = sector(; A=U1(1), B=U1(1)) - @test q00 ⊗ q00 == [q00] - @test q01 ⊗ q00 == [q01] - @test q00 ⊗ q01 == [q01] - @test q10 ⊗ q01 == [q11] + @test (@inferred q10 ⊗ q10) == sector(; A=U1(2)) + @test (@inferred q01 ⊗ q00) == q01 + @test (@inferred q00 ⊗ q01) == q01 + @test (@inferred q10 ⊗ q01) == q11 + @test (@inferred q11 ⊗ q11) == sector(; A=U1(2), B=U1(2)) + + s11 = sector(; A=U1(1), B=Z{2}(1)) + s10 = sector(; A=U1(1)) + s01 = sector(; B=Z{2}(1)) + @test (@inferred s01 ⊗ q00) == s01 + @test (@inferred q00 ⊗ s01) == s01 + @test (@inferred s10 ⊗ s01) == s11 + @test (@inferred s11 ⊗ s11) == sector(; A=U1(2), B=Z{2}(0)) end - @testset "U(1) ⊗ SU(2) conventional" begin - q0 = sector() + @testset "Fusion of NonAbelian products" begin + p0 = sector() + pha = sector(; A=SU2(1//2)) + phb = sector(; B=SU2(1//2)) + phab = sector(; A=SU2(1//2), B=SU2(1//2)) + + @test gradedisequal( + (@inferred pha ⊗ pha), gradedrange([sector(; A=SU2(0)) => 1, sector(; A=SU2(1)) => 1]) + ) + @test gradedisequal((@inferred pha ⊗ p0), gradedrange([pha => 1])) + @test gradedisequal((@inferred p0 ⊗ phb), gradedrange([phb => 1])) + @test gradedisequal((@inferred pha ⊗ phb), gradedrange([phab => 1])) + + @test gradedisequal( + (@inferred phab ⊗ phab), + gradedrange([ + sector(; A=SU2(0), B=SU2(0)) => 1, + sector(; A=SU2(1), B=SU2(0)) => 1, + sector(; A=SU2(0), B=SU2(1)) => 1, + sector(; A=SU2(1), B=SU2(1)) => 1, + ]), + ) + end + + @testset "Fusion of NonGroupCategory products" begin + ı = Fib("1") + τ = Fib("τ") + s = sector(; A=ı, B=ı) + @test gradedisequal((@inferred s ⊗ s), gradedrange([s => 1])) + + s = sector(; A=τ, B=τ) + @test gradedisequal( + (@inferred s ⊗ s), + gradedrange([ + sector(; A=ı, B=ı) => 1, + sector(; A=τ, B=ı) => 1, + sector(; A=ı, B=τ) => 1, + sector(; A=τ, B=τ) => 1, + ]), + ) + + σ = Ising("σ") + ψ = Ising("ψ") + s = sector(; A=τ, B=σ) + g = gradedrange([ + sector(; A=ı, B=Ising("1")) => 1, + sector(; A=τ, B=Ising("1")) => 1, + sector(; A=ı, B=ψ) => 1, + sector(; A=τ, B=ψ) => 1, + ]) + @test gradedisequal((@inferred s ⊗ s), g) + end + + @testset "Fusion of mixed Abelian and NonAbelian products" begin q0h = sector(; J=SU2(1//2)) q10 = (N=U1(1),) × (J=SU2(0),) # Put names in reverse order sometimes: q1h = (J=SU2(1//2),) × (N=U1(1),) q11 = (N=U1(1),) × (J=SU2(1),) - q20 = sector(; N=U1(2)) + q20 = (N=U1(2),) × (J=SU2(0),) # julia 1.6 does not accept gradedrange without J q2h = (N=U1(2),) × (J=SU2(1//2),) q21 = (N=U1(2),) × (J=SU2(1),) q22 = (N=U1(2),) × (J=SU2(2),) - @test q1h ⊗ q1h == q20 ⊕ q21 - @test q10 ⊗ q1h == [q2h] - @test q0h ⊗ q1h == q10 ⊕ q11 - @test q11 ⊗ q11 == q20 ⊕ q21 ⊕ q22 + @test gradedisequal((@inferred q1h ⊗ q1h), gradedrange([q20 => 1, q21 => 1])) + @test gradedisequal((@inferred q10 ⊗ q1h), gradedrange([q2h => 1])) + @test gradedisequal((@inferred q0h ⊗ q1h), gradedrange([q10 => 1, q11 => 1])) + @test gradedisequal((@inferred q11 ⊗ q11), gradedrange([q20 => 1, q21 => 1, q22 => 1])) end - @testset "U(1) ⊗ SU(2)" begin - q0 = sector() - q0h = sector(; J=SU{2}(2)) - q10 = (N=U1(1),) × (J=SU{2}(1),) - # Put names in reverse order sometimes: - q1h = (J=SU{2}(2),) × (N=U1(1),) - q11 = (N=U1(1),) × (J=SU{2}(3),) - q20 = sector(; N=U1(2)) - q2h = (N=U1(2),) × (J=SU{2}(2),) - q21 = (N=U1(2),) × (J=SU{2}(3),) - q22 = (N=U1(2),) × (J=SU{2}(5),) + @testset "Fusion of fully mixed products" begin + s = sector(; A=U1(1), B=SU2(1//2), C=Ising("σ")) + @test gradedisequal( + (@inferred s ⊗ s), + gradedrange([ + sector(; A=U1(2), B=SU2(0), C=Ising("1")) => 1, + sector(; A=U1(2), B=SU2(1), C=Ising("1")) => 1, + sector(; A=U1(2), B=SU2(0), C=Ising("ψ")) => 1, + sector(; A=U1(2), B=SU2(1), C=Ising("ψ")) => 1, + ]), + ) - @test q1h ⊗ q1h == q20 ⊕ q21 - @test q10 ⊗ q1h == [q2h] - @test q0h ⊗ q1h == q10 ⊕ q11 - @test q11 ⊗ q11 == q20 ⊕ q21 ⊕ q22 - end - - @testset "Comparisons with unspecified labels" begin - q2 = sector(; N=U1(2)) - q20 = (N=U1(2),) × (J=SU{2}(1),) - @test q20 == q2 + ı = Fib("1") + τ = Fib("τ") + s = sector(; A=SU2(1//2), B=U1(1), C=τ) + @test gradedisequal( + (@inferred s ⊗ s), + gradedrange([ + sector(; A=SU2(0), B=U1(2), C=ı) => 1, + sector(; A=SU2(1), B=U1(2), C=ı) => 1, + sector(; A=SU2(0), B=U1(2), C=τ) => 1, + sector(; A=SU2(1), B=U1(2), C=τ) => 1, + ]), + ) - q21 = (N=U1(2),) × (J=SU{2}(3),) - @test q21 != q2 + s = sector(; A=τ, B=U1(1), C=ı) + @test gradedisequal( + (@inferred s ⊗ s), + gradedrange([sector(; B=U1(2), A=ı, C=ı) => 1, sector(; B=U1(2), A=τ, C=ı) => 1]), + ) + end + @testset "GradedUnitRange fusion rules" begin + s1 = sector(; A=U1(1), B=SU2(1//2), C=Ising("σ")) + s2 = sector(; A=U1(0), B=SU2(1//2), C=Ising("1")) + g1 = gradedrange([s1 => 2]) + g2 = gradedrange([s2 => 1]) + s3 = sector(; A=U1(1), B=SU2(0), C=Ising("σ")) + s4 = sector(; A=U1(1), B=SU2(1), C=Ising("σ")) + @test gradedisequal( + (@inferred_latest fusion_product(g1, g2)), gradedrange([s3 => 2, s4 => 2]) + ) - a = (A=U1(0),) × (B=U1(2),) - b = (B=U1(2),) × (C=U1(0),) - @test a == b - c = (B=U1(2),) × (C=U1(1),) - @test a != c + sA = sector(; A=U1(1)) + sB = sector(; B=SU2(1//2)) + sAB = sector(; A=U1(1), B=SU2(1//2)) + gA = gradedrange([sA => 2]) + gB = gradedrange([sB => 1]) + @test gradedisequal((@inferred fusion_product(gA, gB)), gradedrange([sAB => 2])) end end -@testset "Test Ordered Products" begin - @testset "Ordered Constructor" begin - s = sector(U1(1), U1(2)) - @test length(s) == 2 - @test s[1] == U1(1) - @test s[2] == U1(2) +@testset "Empty category" begin + s = sector() + @test s == s + @test (@inferred dual(s)) == s + @test (@inferred s × s) == s + @test (@inferred s ⊗ s) == s + @test (@inferred quantum_dimension(s)) == 1 + @test (@inferred_latest trivial(s)) == s + @test typeof(s) == typeof(sector(())) + @test typeof(s) == typeof(sector((;))) # empty NamedTuple is cast to Tuple{} - s = U1(1) × SU2(1//2) × U1(3) - @test length(s) == 3 - @test s[1] == U1(1) - @test s[2] == SU2(1//2) - @test s[3] == U1(3) - end + @test (@inferred s × U1(1)) == sector(U1(1)) + @test (@inferred s × sector(U1(1))) == sector(U1(1)) + @test (@inferred s × sector(; A=U1(1))) == sector(; A=U1(1)) + @test (@inferred U1(1) × s) == sector(U1(1)) + @test (@inferred sector(U1(1)) × s) == sector(U1(1)) + @test (@inferred sector(; A=U1(1)) × s) == sector(; A=U1(1)) - @testset "Fusion of U1 products" begin - p11 = U1(1) × U1(1) - @test p11 ⊗ p11 == [U1(2) × U1(2)] + # Empty acts as trivial + @test (@inferred U1(1) ⊗ s) == U1(1) + @test (@inferred SU2(0) ⊗ s) == gradedrange([SU2(0) => 1]) + @test (@inferred Fib("τ") ⊗ s) == gradedrange([Fib("τ") => 1]) + @test (@inferred s ⊗ U1(1)) == U1(1) + @test (@inferred s ⊗ SU2(0)) == gradedrange([SU2(0) => 1]) + @test (@inferred s ⊗ Fib("τ")) == gradedrange([Fib("τ") => 1]) - p123 = U1(1) × U1(2) × U1(3) - @test p123 ⊗ p123 == [U1(2) × U1(4) × U1(6)] - end - - @testset "Enforce same number of spaces" begin - p12 = U1(1) × U1(2) - p123 = U1(1) × U1(2) × U1(3) - @test_throws DimensionMismatch p12 ⊗ p123 - end + @test (@inferred sector(U1(1)) ⊗ s) == sector(U1(1)) + @test (@inferred sector(SU2(0)) ⊗ s) == gradedrange([sector(SU2(0)) => 1]) + @test (@inferred sector(Fib("τ"), SU2(1), U1(2)) ⊗ s) == + gradedrange([sector(Fib("τ"), SU2(1), U1(2)) => 1]) - @testset "Fusion of SU2 products" begin - phh = SU2(1//2) × SU2(1//2) - @test phh ⊗ phh == - (SU2(0) × SU2(0)) ⊕ (SU2(1) × SU2(0)) ⊕ (SU2(0) × SU2(1)) ⊕ (SU2(1) × SU2(1)) - end + @test (@inferred sector(; A=U1(1)) ⊗ s) == sector(; A=U1(1)) + @test (@inferred sector(; A=SU2(0)) ⊗ s) == gradedrange([sector(; A=SU2(0)) => 1]) + @test (@inferred sector(; A=Fib("τ"), B=SU2(1), C=U1(2)) ⊗ s) == + gradedrange([sector(; A=Fib("τ"), B=SU2(1), C=U1(2)) => 1]) - @testset "Fusion of mixed U1 and SU2 products" begin - p2h = U1(2) × SU2(1//2) - p1h = U1(1) × SU2(1//2) - @test p2h ⊗ p1h == (U1(3) × SU2(0)) ⊕ (U1(3) × SU2(1)) + # Empty behaves as empty NamedTuple + @test s != U1(0) + @test s != sector(U1(0)) + @test s != sector(; A=U1(1)) + @test s == sector(; A=U1(0)) + @test sector(; A=U1(0)) == s - p1h1 = U1(1) × SU2(1//2) × Z{2}(1) - @test p1h1 ⊗ p1h1 == (U1(2) × SU2(0) × Z{2}(0)) ⊕ (U1(2) × SU2(1) × Z{2}(0)) - end + @test !(s < s) + @test_throws ArgumentError s < sector(U1(0)) + @test s < sector(; A=U1(1)) + @test s > sector(; A=U1(-1)) + @test !(s < sector(; A=U1(0))) + @test !(s > sector(; A=U1(0))) end end diff --git a/NDTensors/src/lib/Sectors/test/test_fusion_rules.jl b/NDTensors/src/lib/Sectors/test/test_fusion_rules.jl new file mode 100644 index 0000000000..c689bc7c3e --- /dev/null +++ b/NDTensors/src/lib/Sectors/test/test_fusion_rules.jl @@ -0,0 +1,242 @@ +@eval module $(gensym()) +using NDTensors.GradedAxes: + dual, fusion_product, gradedisequal, gradedrange, flip, tensor_product +using NDTensors.Sectors: + ⊗, Fib, Ising, SU, SU2, U1, Z, block_dimensions, quantum_dimension, trivial +using Test: @inferred, @test, @testset, @test_throws + +@testset "Simple object fusion rules" begin + @testset "Z{2} fusion rules" begin + z0 = Z{2}(0) + z1 = Z{2}(1) + + @test z0 ⊗ z0 == z0 + @test z0 ⊗ z1 == z1 + @test z1 ⊗ z1 == z0 + @test (@inferred z0 ⊗ z0) == z0 # no better way, see Julia PR 23426 + + # using GradedAxes interface + @test gradedisequal(fusion_product(z0, z0), gradedrange([z0 => 1])) + @test gradedisequal(fusion_product(z0, z1), gradedrange([z1 => 1])) + + # test different input number + @test gradedisequal(fusion_product(z0), gradedrange([z0 => 1])) + @test gradedisequal(fusion_product(z0, z0, z0), gradedrange([z0 => 1])) + @test gradedisequal(fusion_product(z0, z0, z0, z0), gradedrange([z0 => 1])) + @test (@inferred block_dimensions(gradedrange([z1 => 1]))) == [1] + end + @testset "U(1) fusion rules" begin + q1 = U1(1) + q2 = U1(2) + q3 = U1(3) + + @test q1 ⊗ q1 == U1(2) + @test q1 ⊗ q2 == U1(3) + @test q2 ⊗ q1 == U1(3) + @test (@inferred q1 ⊗ q2) == q3 # no better way, see Julia PR 23426 + end + @testset "SU2 fusion rules" begin + j1 = SU2(0) + j2 = SU2(1//2) + j3 = SU2(1) + j4 = SU2(3//2) + j5 = SU2(2) + + @test gradedisequal(j1 ⊗ j2, gradedrange([j2 => 1])) + @test gradedisequal(j2 ⊗ j2, gradedrange([j1 => 1, j3 => 1])) + @test gradedisequal(j2 ⊗ j3, gradedrange([j2 => 1, j4 => 1])) + @test gradedisequal(j3 ⊗ j3, gradedrange([j1 => 1, j3 => 1, j5 => 1])) + @test gradedisequal((@inferred j1 ⊗ j2), gradedrange([j2 => 1])) + @test (@inferred quantum_dimension(j1 ⊗ j2)) == 2 + @test (@inferred block_dimensions(j1 ⊗ j2)) == [2] + + @test gradedisequal(fusion_product(j2), gradedrange([j2 => 1])) + @test gradedisequal(fusion_product(j2, j1), gradedrange([j2 => 1])) + @test gradedisequal(fusion_product(j2, j1, j1), gradedrange([j2 => 1])) + end + + @testset "Fibonacci fusion rules" begin + ı = Fib("1") + τ = Fib("τ") + + @test gradedisequal(ı ⊗ ı, gradedrange([ı => 1])) + @test gradedisequal(ı ⊗ τ, gradedrange([τ => 1])) + @test gradedisequal(τ ⊗ ı, gradedrange([τ => 1])) + @test gradedisequal((@inferred τ ⊗ τ), gradedrange([ı => 1, τ => 1])) + @test (@inferred quantum_dimension(gradedrange([ı => 1, ı => 1]))) == 2.0 + end + + @testset "Ising fusion rules" begin + ı = Ising("1") + σ = Ising("σ") + ψ = Ising("ψ") + + @test gradedisequal(ı ⊗ ı, gradedrange([ı => 1])) + @test gradedisequal(ı ⊗ σ, gradedrange([σ => 1])) + @test gradedisequal(σ ⊗ ı, gradedrange([σ => 1])) + @test gradedisequal(ı ⊗ ψ, gradedrange([ψ => 1])) + @test gradedisequal(ψ ⊗ ı, gradedrange([ψ => 1])) + @test gradedisequal(σ ⊗ σ, gradedrange([ı => 1, ψ => 1])) + @test gradedisequal(σ ⊗ ψ, gradedrange([σ => 1])) + @test gradedisequal(ψ ⊗ σ, gradedrange([σ => 1])) + @test gradedisequal(ψ ⊗ ψ, gradedrange([ı => 1])) + @test gradedisequal((@inferred ψ ⊗ ψ), gradedrange([ı => 1])) + @test (@inferred quantum_dimension(σ ⊗ σ)) == 2.0 + end +end +@testset "Reducible object fusion rules" begin + @testset "Trivial GradedUnitRange" begin + g1 = gradedrange([U1(0) => 1]) + g2 = gradedrange([SU2(0) => 1]) + @test gradedisequal(trivial(g1), g1) + @test gradedisequal(trivial(dual(g1)), g1) # trivial returns nondual + @test gradedisequal(trivial(typeof(g2)), g2) + end + @testset "GradedUnitRange abelian tensor/fusion product" begin + g1 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 2]) + g2 = gradedrange([U1(-2) => 2, U1(0) => 1, U1(1) => 2]) + + @test gradedisequal(flip(dual(g1)), gradedrange([U1(1) => 1, U1(0) => 1, U1(-1) => 2])) + @test (@inferred block_dimensions(g1)) == [1, 1, 2] + + gt = gradedrange([ + U1(-3) => 2, + U1(-2) => 2, + U1(-1) => 4, + U1(-1) => 1, + U1(0) => 1, + U1(1) => 2, + U1(0) => 2, + U1(1) => 2, + U1(2) => 4, + ]) + gf = gradedrange([ + U1(-3) => 2, U1(-2) => 2, U1(-1) => 5, U1(0) => 3, U1(1) => 4, U1(2) => 4 + ]) + @test gradedisequal((@inferred tensor_product(g1, g2)), gt) + @test gradedisequal((@inferred fusion_product(g1, g2)), gf) + + gtd1 = gradedrange([ + U1(-1) => 2, + U1(-2) => 2, + U1(-3) => 4, + U1(1) => 1, + U1(0) => 1, + U1(-1) => 2, + U1(2) => 2, + U1(1) => 2, + U1(0) => 4, + ]) + gfd1 = gradedrange([ + U1(-3) => 4, U1(-2) => 2, U1(-1) => 4, U1(0) => 5, U1(1) => 3, U1(2) => 2 + ]) + @test gradedisequal((@inferred tensor_product(dual(g1), g2)), gtd1) + @test gradedisequal((@inferred fusion_product(dual(g1), g2)), gfd1) + + gtd2 = gradedrange([ + U1(1) => 2, + U1(2) => 2, + U1(3) => 4, + U1(-1) => 1, + U1(0) => 1, + U1(1) => 2, + U1(-2) => 2, + U1(-1) => 2, + U1(0) => 4, + ]) + gfd2 = gradedrange([ + U1(-2) => 2, U1(-1) => 3, U1(0) => 5, U1(1) => 4, U1(2) => 2, U1(3) => 4 + ]) + @test gradedisequal((@inferred tensor_product(g1, dual(g2))), gtd2) + @test gradedisequal((@inferred fusion_product(g1, dual(g2))), gfd2) + + gtd = gradedrange([ + U1(3) => 2, + U1(2) => 2, + U1(1) => 4, + U1(1) => 1, + U1(0) => 1, + U1(-1) => 2, + U1(0) => 2, + U1(-1) => 2, + U1(-2) => 4, + ]) + gfd = gradedrange([ + U1(-2) => 4, U1(-1) => 4, U1(0) => 3, U1(1) => 5, U1(2) => 2, U1(3) => 2 + ]) + @test gradedisequal((@inferred tensor_product(dual(g1), dual(g2))), gtd) + @test gradedisequal((@inferred fusion_product(dual(g1), dual(g2))), gfd) + + # test different (non-product) categories cannot be fused + @test_throws MethodError fusion_product(gradedrange([Z{2}(0) => 1]), g1) + @test_throws MethodError tensor_product(gradedrange([Z{2}(0) => 1]), g2) + end + + @testset "GradedUnitRange non-abelian fusion rules" begin + g3 = gradedrange([SU2(0) => 1, SU2(1//2) => 2, SU2(1) => 1]) + g4 = gradedrange([SU2(1//2) => 1, SU2(1) => 2]) + g34 = gradedrange([ + SU2(1//2) => 1, + SU2(0) => 2, + SU2(1) => 2, + SU2(1//2) => 1, + SU2(3//2) => 1, + SU2(1) => 2, + SU2(1//2) => 4, + SU2(3//2) => 4, + SU2(0) => 2, + SU2(1) => 2, + SU2(2) => 2, + ]) + + @test gradedisequal(tensor_product(g3, g4), g34) + + @test gradedisequal(dual(flip(g3)), g3) # trivial for SU(2) + @test gradedisequal( + (@inferred fusion_product(g3, g4)), + gradedrange([SU2(0) => 4, SU2(1//2) => 6, SU2(1) => 6, SU2(3//2) => 5, SU2(2) => 2]), + ) + @test (@inferred block_dimensions(g3)) == [1, 4, 3] + + # test dual on non self-conjugate non-abelian representations + s1 = SU{3}((0, 0)) + f3 = SU{3}((1, 0)) + c3 = SU{3}((1, 1)) + ad8 = SU{3}((2, 1)) + + g5 = gradedrange([s1 => 1, f3 => 1]) + g6 = gradedrange([s1 => 1, c3 => 1]) + @test gradedisequal(dual(flip(g5)), g6) + @test gradedisequal( + fusion_product(g5, g6), gradedrange([s1 => 2, f3 => 1, c3 => 1, ad8 => 1]) + ) + @test gradedisequal( + fusion_product(dual(g5), g6), + gradedrange([s1 => 1, f3 => 1, c3 => 2, SU{3}((2, 2)) => 1]), + ) + @test gradedisequal( + fusion_product(g5, dual(g6)), + gradedrange([s1 => 1, f3 => 2, c3 => 1, SU{3}((2, 0)) => 1]), + ) + @test gradedisequal( + fusion_product(dual(g5), dual(g6)), gradedrange([s1 => 2, f3 => 1, c3 => 1, ad8 => 1]) + ) + end + + @testset "Mixed GradedUnitRange - Category fusion rules" begin + g1 = gradedrange([U1(1) => 1, U1(2) => 2]) + g2 = gradedrange([U1(2) => 1, U1(3) => 2]) + @test gradedisequal((@inferred fusion_product(g1, U1(1))), g2) + @test gradedisequal((@inferred fusion_product(U1(1), g1)), g2) + + g3 = gradedrange([SU2(0) => 1, SU2(1//2) => 2]) + g4 = gradedrange([SU2(0) => 2, SU2(1//2) => 1, SU2(1) => 2]) + @test gradedisequal((@inferred fusion_product(g3, SU2(1//2))), g4) + @test gradedisequal((@inferred fusion_product(SU2(1//2), g3)), g4) + + # test different categories cannot be fused + @test_throws MethodError fusion_product(g1, SU2(1)) + @test_throws MethodError fusion_product(U1(1), g3) + end +end +end diff --git a/NDTensors/src/lib/Sectors/test/test_simple_categories.jl b/NDTensors/src/lib/Sectors/test/test_simple_categories.jl new file mode 100644 index 0000000000..7553367b87 --- /dev/null +++ b/NDTensors/src/lib/Sectors/test/test_simple_categories.jl @@ -0,0 +1,141 @@ +@eval module $(gensym()) +using NDTensors.GradedAxes: dual +using NDTensors.Sectors: + Fib, Ising, SU, SU2, U1, Z, adjoint, quantum_dimension, fundamental, istrivial, trivial +using Test: @inferred, @test, @testset, @test_throws +@testset "Test Category Types" begin + @testset "U(1)" begin + q1 = U1(1) + q2 = U1(2) + q3 = U1(3) + + @test quantum_dimension(q1) == 1 + @test quantum_dimension(q2) == 1 + @test (@inferred quantum_dimension(q1)) == 1 + + @test trivial(q1) == U1(0) + @test trivial(U1) == U1(0) + @test istrivial(U1(0)) + + @test dual(U1(2)) == U1(-2) + @test isless(U1(1), U1(2)) + @test !isless(U1(2), U1(1)) + end + + @testset "Z₂" begin + z0 = Z{2}(0) + z1 = Z{2}(1) + + @test trivial(Z{2}) == Z{2}(0) + @test istrivial(Z{2}(0)) + + @test quantum_dimension(z0) == 1 + @test quantum_dimension(z1) == 1 + @test (@inferred quantum_dimension(z0)) == 1 + + @test dual(z0) == z0 + @test dual(z1) == z1 + + @test dual(Z{2}(1)) == Z{2}(1) + @test isless(Z{2}(0), Z{2}(1)) + @test !isless(Z{2}(1), Z{2}(0)) + end + + @testset "SU2" begin + j1 = SU2(0) + j2 = SU2(1//2) # Rational will be cast to HalfInteger + j3 = SU2(1) + j4 = SU2(3//2) + + # alternative constructors + @test j2 == SU{2}((1,)) # tuple SU(N)-like constructor + @test j2 == SU{2,1}((1,)) # tuple constructor with explicit {N,N-1} + @test j2 == SU((1,)) # infer N from tuple length + @test j2 == SU{2}((Int8(1),)) # any Integer type accepted + @test j2 == SU{2}((UInt32(1),)) # any Integer type accepted + @test j2 == SU2(1 / 2) # Float will be cast to HalfInteger + @test_throws MethodError SU2((1,)) # avoid confusion between tuple and half-integer interfaces + @test_throws MethodError SU{2,1}(1) # avoid confusion + + @test trivial(SU{2}) == SU2(0) + @test istrivial(SU2(0)) + + @test fundamental(SU{2}) == SU2(1//2) + @test adjoint(SU{2}) == SU2(1) + + @test quantum_dimension(j1) == 1 + @test quantum_dimension(j2) == 2 + @test quantum_dimension(j3) == 3 + @test quantum_dimension(j4) == 4 + @test (@inferred quantum_dimension(j1)) == 1 + + @test dual(j1) == j1 + @test dual(j2) == j2 + @test dual(j3) == j3 + @test dual(j4) == j4 + end + + @testset "SU(N)" begin + f3 = SU{3}((1, 0)) + f4 = SU{4}((1, 0, 0)) + ad3 = SU{3}((2, 1)) + ad4 = SU{4}((2, 1, 1)) + + @test trivial(SU{3}) == SU{3}((0, 0)) + @test istrivial(SU{3}((0, 0))) + @test trivial(SU{4}) == SU{4}((0, 0, 0)) + @test istrivial(SU{4}((0, 0, 0))) + + @test fundamental(SU{3}) == f3 + @test adjoint(SU{3}) == ad3 + @test fundamental(SU{4}) == f4 + @test adjoint(SU{4}) == ad4 + + @test dual(f3) == SU{3}((1, 1)) + @test dual(f4) == SU{4}((1, 1, 1)) + @test dual(ad3) == ad3 + @test dual(ad4) == ad4 + + @test quantum_dimension(f3) == 3 + @test quantum_dimension(f4) == 4 + @test quantum_dimension(ad3) == 8 + @test quantum_dimension(ad4) == 15 + @test quantum_dimension(SU{3}((4, 2))) == 27 + @test quantum_dimension(SU{3}((3, 3))) == 10 + @test quantum_dimension(SU{3}((3, 0))) == 10 + @test quantum_dimension(SU{3}((0, 0))) == 1 + @test (@inferred quantum_dimension(f3)) == 3 + end + + @testset "Fibonacci" begin + ı = Fib("1") + τ = Fib("τ") + + @test trivial(Fib) == ı + @test istrivial(ı) + + @test dual(ı) == ı + @test dual(τ) == τ + + @test (@inferred quantum_dimension(ı)) == 1.0 + @test (@inferred quantum_dimension(τ)) == ((1 + √5) / 2) + end + + @testset "Ising" begin + ı = Ising("1") + σ = Ising("σ") + ψ = Ising("ψ") + + @test trivial(Ising) == ı + @test istrivial(ı) + + @test dual(ı) == ı + @test dual(σ) == σ + @test dual(ψ) == ψ + + @test (@inferred quantum_dimension(ı)) == 1.0 + @test (@inferred quantum_dimension(σ)) == √2 + @test (@inferred quantum_dimension(ψ)) == 1.0 + end +end +end