diff --git a/Project.toml b/Project.toml index a542085..341005b 100644 --- a/Project.toml +++ b/Project.toml @@ -4,25 +4,33 @@ authors = ["Boris De Vos, Laurens Lootens and Lukas Devos"] version = "0.1.0" [deps] -Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" -JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f" [compat] Aqua = "0.8.9" -Artifacts = "1.10, 1" -JSON3 = "1.14.1" +DelimitedFiles = "1.9" +LinearAlgebra = "1" +Pkg = "1" +Random = "1" SafeTestsets = "0.1" -TensorKitSectors = "0.1.2" +TensorKit = "0.16" +TensorKitSectors = "0.3.7" +TensorOperations = "5" Test = "1.10" TestExtras = "0.3" julia = "1.10" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a" [targets] -test = ["Aqua", "SafeTestsets", "Test", "TestExtras"] +test = ["Aqua", "LinearAlgebra", "Random", "SafeTestsets", "TensorOperations", "Test", "TestExtras"] \ No newline at end of file diff --git a/src/MultiTensorKit.jl b/src/MultiTensorKit.jl index 34d66bc..635ac69 100644 --- a/src/MultiTensorKit.jl +++ b/src/MultiTensorKit.jl @@ -2,10 +2,14 @@ module MultiTensorKit export BimoduleSector, A4Object -using JSON3 -using Artifacts +using DelimitedFiles +using Pkg +using Pkg.Artifacts using TensorKitSectors +using TensorKit +import TensorKit: hasblock, dim + include("bimodulesector.jl") end diff --git a/src/bimodulesector.jl b/src/bimodulesector.jl index fd52d31..ea6cbc1 100644 --- a/src/bimodulesector.jl +++ b/src/bimodulesector.jl @@ -1,80 +1,118 @@ +const ImplementedBimoduleSectors = [:A4] + struct BimoduleSector{Name} <: Sector i::Int j::Int label::Int + function BimoduleSector{Name}(i::Int, j::Int, label::Int) where {Name} + Name ∈ ImplementedBimoduleSectors || + throw(ArgumentError("BimoduleSector $Name not implemented")) + i <= size(BimoduleSector{Name}) && j <= size(BimoduleSector{Name}) || + throw(DomainError("object outside the matrix $Name")) + return label <= _numlabels(BimoduleSector{Name}, i, j) ? new{Name}(i, j, label) : + throw(DomainError("label outside category $Name($i, $j)")) + end end -BimoduleSector{Name}(data::NTuple{3,Int}) where {Name} = BimoduleSector{Name}(data...) + +BimoduleSector{Name}(data::NTuple{3, Int}) where {Name} = BimoduleSector{Name}(data...) +BimoduleSectorName(::Type{BimoduleSector{Name}}) where {Name} = Name const A4Object = BimoduleSector{:A4} +Base.convert(::Type{<:BimoduleSector{Name}}, labels::NTuple{3, Int}) where {Name} = BimoduleSector{Name}(labels...) + +function Base.show(io::IO, a::BimoduleSector{Name}) where {Name} + if get(io, :typeinfo, nothing) === typeof(a) + print(io, (a.i, a.j, a.label)) + else + print(io, typeof(a), (a.i, a.j, a.label)) + end + return nothing +end + # Utility implementations # ----------------------- -function Base.isless(a::I, b::I) where {I<:BimoduleSector} +function Base.isless(a::I, b::I) where {I <: BimoduleSector} return isless((a.i, a.j, a.label), (b.i, b.j, b.label)) end Base.hash(a::BimoduleSector, h::UInt) = hash(a.i, hash(a.j, hash(a.label, h))) -function Base.convert(::Type{BimoduleSector{Name}}, d::NTuple{3,Int}) where {Name} +function Base.convert(::Type{BimoduleSector{Name}}, d::NTuple{3, Int}) where {Name} return BimoduleSector{Name}(d...) end -Base.IteratorSize(::Type{SectorValues{<:BimoduleSector}}) = Base.SizeUnknown() +Base.size(::Type{A4Object}) = 7 + +Base.IteratorSize(::Type{<:SectorValues{<:BimoduleSector}}) = Base.SizeUnknown() -# TODO: generalize? -# TODO: numlabels? -function Base.iterate(iter::SectorValues{A4Object}, (I, label)=(1, 1)) - I > 12 * 12 && return nothing - i, j = CartesianIndices((12, 12))[I].I - maxlabel = numlabels(A4Object, i, j) +function Base.iterate(iter::SectorValues{<:BimoduleSector}, (I, label) = (1, 1)) + A = eltype(iter) + s = size(A) + I > s * s && return nothing + i, j = CartesianIndices((s, s))[I].I + maxlabel = _numlabels(A, i, j) return if label > maxlabel iterate(iter, (I + 1, 1)) else - A4Object(i, j, label), (I, label + 1) + A(i, j, label), (I, label + 1) end end +function Base.length(::SectorValues{I}) where {I <: BimoduleSector} + s = size(I) + return sum(_numlabels(I, i, j) for i in 1:s, j in 1:s) +end + TensorKitSectors.FusionStyle(::Type{A4Object}) = GenericFusion() -TensorKitSectors.BraidingStyle(::Type{A4Object}) = NoBraiding() +TensorKitSectors.BraidingStyle(::Type{<:BimoduleSector}) = NoBraiding() +TensorKitSectors.sectorscalartype(::Type{A4Object}) = ComplexF64 -function TensorKitSectors.:βŠ—(a::A4Object, b::A4Object) - @assert a.j == b.i - Ncache = _get_Ncache(A4Object)[a.i, a.j, b.j] - return A4Object[A4Object(a.i, b.j, c_l) - for (a_l, b_l, c_l) in keys(Ncache) - if (a_l == a.label && b_l == b.label)] +function TensorKitSectors.:βŠ—(a::I, b::I) where {I <: BimoduleSector} + # @assert a.j == b.i + a.j == b.i || return I[] + Ncache = _get_Ncache(I)[a.i, a.j, b.j] + return I[ + I(a.i, b.j, c_l) for (a_l, b_l, c_l) in keys(Ncache) + if (a_l == a.label && b_l == b.label) + ] end -function _numlabels(::Type{A4Object}, i, j) - return Ncache = _get_Ncache(A4Object) +function _numlabels(::Type{T}, i, j) where {T <: BimoduleSector} + return length(_get_dual_cache(T)[2][i, j]) end +# User-friendly functions +# ------------------- +#TODO: add functions to identify categories + # Data from files # --------------- const artifact_path = joinpath(artifact"fusiondata", "MultiTensorKit.jl-data-v0.1.5") -function extract_Nsymbol(::Type{A4Object}) - filename = joinpath(artifact_path, "A4", "Nsymbol.json") - isfile(filename) || throw(LoadError(filename, 0, "Nsymbol file not found for $Name")) - json_string = read(filename, String) - Narray = copy(JSON3.read(json_string)) - return map(reshape(Narray, 12, 12, 12)) do x - y = Dict{NTuple{3,Int},Int}() - for (k, v) in x - a, b, c = parse.(Int, split(string(k)[2:(end - 1)], ", ")) - y[(a, b, c)] = v - end - return y +function extract_Nsymbol(::Type{I}) where {I <: BimoduleSector} + name = string(BimoduleSectorName(I)) + filename = joinpath(artifact_path, name, "Nsymbol.txt") + isfile(filename) || throw(LoadError(filename, 0, "Nsymbol file not found for $name")) + Narray = readdlm(filename) # matrix with 7 columns + + data_dict = Dict{NTuple{3, Int}, Dict{NTuple{3, Int}, Int}}() + for row in eachrow(Narray) + i, j, k, a, b, c, N = Int.(@view(row[1:size(I)])) + colordict = get!(data_dict, (i, j, k), Dict{NTuple{3, Int}, Int}()) + push!(colordict, (a, b, c) => N) end + + return data_dict end -const Ncache = IdDict{Type{<:BimoduleSector},Array{Dict{NTuple{3,Int},Int},3}}() +const Ncache = IdDict{Type{<:BimoduleSector}, Dict{NTuple{3, Int}, Dict{NTuple{3, Int}, Int}}}() -function _get_Ncache(::Type{T}) where {T<:BimoduleSector} +function _get_Ncache(::Type{T}) where {T <: BimoduleSector} global Ncache return get!(Ncache, T) do return extract_Nsymbol(T) end end -function TensorKitSectors.Nsymbol(a::I, b::I, c::I) where {I<:A4Object} +function TensorKitSectors.Nsymbol(a::I, b::I, c::I) where {I <: BimoduleSector} # TODO: should this error or return 0? (a.j == b.i && a.i == c.i && b.j == c.j) || throw(ArgumentError("invalid fusion channel")) @@ -82,107 +120,196 @@ function TensorKitSectors.Nsymbol(a::I, b::I, c::I) where {I<:A4Object} return get(_get_Ncache(I)[i, j, k], (a.label, b.label, c.label), 0) end -# TODO: can we define dual for modules? -const Dualcache = IdDict{Type{<:BimoduleSector},Vector{Tuple{Int,Vector{Int}}}}() +const Dualcache = IdDict{Type{<:BimoduleSector}, Tuple{Vector{Int64}, Matrix{Vector{Int64}}}}() -function _get_dual_cache(::Type{T}) where {T<:BimoduleSector} +function _get_dual_cache(::Type{T}) where {T <: BimoduleSector} global Dualcache return get!(Dualcache, T) do return extract_dual(T) end end -function extract_dual(::Type{A4Object}) - N = _get_Ncache(A4Object) - ncats = size(N, 1) - return map(1:ncats) do i +function extract_dual(::Type{I}) where {I <: BimoduleSector} + N = _get_Ncache(I) + ncats = size(I) + Is = zeros(Int, ncats) + + map(1:ncats) do i Niii = N[i, i, i] - nobj = maximum(first, keys(Niii)) - - # find identity object: - # I x I -> I, a x I -> a, I x a -> a - I = findfirst(1:nobj) do u - get(Niii, (u, u, u), 0) == 1 || return false - for j in 1:nobj - get(Niii, (j, u, j), 0) == 1 || return false - get(Niii, (u, j, j), 0) == 1 || return false + nobji = maximum(first, keys(N[i, i, i])) + # want to return a leftunit and rightunit for each entry in multifusion cat + # leftunit/rightunit needs to at least be the unit object within a fusion cat + Is[i] = findfirst(1:nobji) do a + get(Niii, (a, a, a), 0) == 1 || return false # I x I -> I + for othera in 1:nobji + get(Niii, (othera, a, othera), 0) == 1 || return false # a x I -> a + get(Niii, (a, othera, othera), 0) == 1 || return false # I x a -> a + end + + # check leftunit + map(1:ncats) do j + nobjj = maximum(first, keys(N[j, j, j])) + for b in 1:nobjj + get(N[i, j, j], (a, b, b), 0) == 1 || return false # I = leftunit(b) + end + end + + # check rightunit + map(1:ncats) do k + nobjk = maximum(first, keys(N[k, k, k])) + for c in 1:nobjk + get(N[k, i, k], (c, a, c), 0) == 1 || return false # I = rightunit(c) + end end return true end + end - # find duals - # a x abar -> I - duals = map(1:nobj) do j - return findfirst(1:nobj) do k - return get(Niii, (j, k, I), 0) == 1 + allduals = Matrix{Vector{Int}}(undef, ncats, ncats) # ncats square matrix of vectors + for i in 1:ncats + nobji = maximum(first, keys(N[i, i, i])) + for j in 1:ncats + allduals[i, j] = Int[] + + nobjj = maximum(first, keys(N[j, j, j])) + # the nested vectors contain the duals of the objects in π’ž_ij, which are in C_ji + Niji = N[i, j, i] # π’ž_ij x π’ž_ji -> C_ii + Njij = N[j, i, j] # π’ž_ji x π’ž_ij -> C_jj + for i_ob in 1:nobji, j_ob in 1:nobjj + get(Niji, (i_ob, j_ob, Is[i]), 0) == 1 || continue # leftunit(c_ij) ∈ c_ij x c_ji + get(Njij, (j_ob, i_ob, Is[j]), 0) == 1 || continue # rightunit(c_ij) ∈ c_ji x c_ij + push!(allduals[i, j], j_ob) end end - - return I, duals end + return Is, allduals end -function Base.one(a::BimoduleSector) - a.i == a.j || error("don't know how to define one for modules") - return A4Object(a.i, a.i, _get_dual_cache(typeof(a))[a.i][1]) -end - -function Base.conj(a::BimoduleSector) - a.i == a.j || error("don't know how to define dual for modules") - return A4Object(a.i, a.i, _get_dual_cache(typeof(a))[a.i][2][a.label]) -end - -function extract_Fsymbol(::Type{A4Object}) - return mapreduce((x, y) -> cat(x, y; dims=1), 1:12) do i - filename = joinpath(artifact_path, "A4", "Fsymbol_$i.json") - @assert isfile(filename) "cannot find $filename" - json_string = read(filename, String) - Farray_part = copy(JSON3.read(json_string)) - return map(enumerate(Farray_part[Symbol(i)])) do (I, x) - j, k, l = Tuple(CartesianIndices((12, 12, 12))[I]) - y = Dict{NTuple{6,Int},Array{ComplexF64,4}}() - for (key, v) in x - a, b, c, d, e, f = parse.(Int, split(string(key)[2:(end - 1)], ", ")) - a_ob, b_ob, c_ob, d_ob, e_ob, f_ob = A4Object.(((i, j, a), (j, k, b), - (k, l, c), (i, l, d), - (i, k, e), (j, l, f))) - result = Array{ComplexF64,4}(undef, - (Nsymbol(a_ob, b_ob, e_ob), - Nsymbol(e_ob, c_ob, d_ob), - Nsymbol(b_ob, c_ob, f_ob), - Nsymbol(a_ob, f_ob, d_ob))) - map!(result, reshape(v, size(result))) do cmplxdict - return complex(cmplxdict[:re], cmplxdict[:im]) - end +function TensorKitSectors.unit(a::BimoduleSector) + a.i == a.j || throw(DomainError("unit of module category ($(a.i), $(a.j)) of $(typeof(a)) is ill-defined")) + return typeof(a)(a.i, a.i, _get_dual_cache(typeof(a))[1][a.i]) +end + +function TensorKitSectors.allunits(::Type{I}) where {I <: BimoduleSector} + s = size(I) + return I[I(i, i, _get_dual_cache(I)[1][i]) for i in 1:s] +end + +function TensorKitSectors.unit(::Type{<:BimoduleSector}) + throw(ArgumentError("unit of Type BimoduleSector doesn't exist")) +end + +function TensorKitSectors.leftunit(a::BimoduleSector) + return typeof(a)(a.i, a.i, _get_dual_cache(typeof(a))[1][a.i]) +end - y[(a, b, c, d, e, f)] = result +function TensorKitSectors.rightunit(a::BimoduleSector) + return typeof(a)(a.j, a.j, _get_dual_cache(typeof(a))[1][a.j]) +end + +function TensorKitSectors.dual(a::BimoduleSector) + return typeof(a)(a.j, a.i, _get_dual_cache(typeof(a))[2][a.i, a.j][a.label]) +end + +function extract_Fsymbol(::Type{I}) where {I <: BimoduleSector} + result = Dict{NTuple{4, Int}, Dict{NTuple{6, Int}, Array{ComplexF64, 4}}}() + name = string(BimoduleSectorName(I)) + filename = joinpath(artifact_path, name, "Fsymbol.txt") + @assert isfile(filename) "cannot find $filename" + + Farray = readdlm(filename) + for ((i, j, k, l), colordict) in convert_Fs(Farray) + result[(i, j, k, l)] = Dict{NTuple{6, Int}, Array{ComplexF64, 4}}() + for ((a, b, c, d, e, f), Fvals) in colordict + a_ob, b_ob, c_ob, d_ob, e_ob, f_ob = I.( + ((i, j, a), (j, k, b), (k, l, c), (i, l, d), (i, k, e), (j, l, f)) + ) + result[(i, j, k, l)][(a, b, c, d, e, f)] = zeros( + ComplexF64, Nsymbol(a_ob, b_ob, e_ob), Nsymbol(e_ob, c_ob, d_ob), + Nsymbol(b_ob, c_ob, f_ob), Nsymbol(a_ob, f_ob, d_ob) + ) + for (K, v) in Fvals + result[(i, j, k, l)][(a, b, c, d, e, f)][K] = v end end end + return result +end + +function convert_Fs(Farray_part::Matrix{Float64}) # Farray_part is a matrix with 16 columns + data_dict = Dict{NTuple{4, Int}, Dict{NTuple{6, Int}, Vector{Pair{CartesianIndex{4}, ComplexF64}}}}() + # want to make a Dict with keys (i,j,k,l) and vals + # a Dict with keys (a,b,c,d,e,f) and vals + # a pair of (mu, nu, rho, sigma) and the F value + for row in eachrow(Farray_part) + i, j, k, l, a, b, c, d, e, f, mu, nu, rho, sigma = Int.(@view(row[1:14])) + v = complex(row[15], row[16]) + colordict = get!( + data_dict, (i, j, k, l), Dict{NTuple{6, Int}, Vector{Pair{CartesianIndex{4}, ComplexF64}}}() + ) + Fdict = get!( + colordict, (a, b, c, d, e, f), Vector{Pair{CartesianIndex{4}, ComplexF64}}() + ) + push!(Fdict, CartesianIndex(mu, nu, rho, sigma) => v) + end + return data_dict end -const Fcache = IdDict{Type{<:BimoduleSector}, - Array{Dict{NTuple{6,Int},Array{ComplexF64,4}},4}}() +const Fcache = IdDict{Type{<:BimoduleSector}, Dict{NTuple{4, Int64}, Dict{NTuple{6, Int64}, Array{ComplexF64, 4}}}}() -function _get_Fcache(::Type{T}) where {T<:BimoduleSector} +function _get_Fcache(::Type{T}) where {T <: BimoduleSector} global Fcache return get!(Fcache, T) do return extract_Fsymbol(T) end end -function TensorKitSectors.Fsymbol(a::I, b::I, c::I, d::I, e::I, - f::I) where {I<:A4Object} - # TODO: should this error or return 0? - (a.j == b.i && b.j == c.i && a.i == d.i && c.j == d.j && - a.i == e.i && b.j == e.j && f.i == a.j && f.j == c.j) || - throw(ArgumentError("invalid fusion channel")) +function TensorKitSectors.Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I <: BimoduleSector} + # required to keep track of multiplicities where F-move is partially unallowed + # also deals with invalid fusion channels + Nabe = Nsymbol(a, b, e) + Necd = Nsymbol(e, c, d) + Nbcf = Nsymbol(b, c, f) + Nafd = Nsymbol(a, f, d) + + zero_array = zeros(sectorscalartype(I), Nabe, Necd, Nbcf, Nafd) + Nabe > 0 && Necd > 0 && Nbcf > 0 && Nafd > 0 || + return zero_array i, j, k, l = a.i, a.j, b.j, c.j - return get(_get_Fcache(I)[i, j, k, l], - (a.label, b.label, c.label, d.label, e.label, f.label)) do - return zeros(sectorscalartype(A4Object), - (Nsymbol(a, b, e), Nsymbol(e, c, d), Nsymbol(b, c, f), - Nsymbol(a, f, d))) + colordict = _get_Fcache(I)[i, j, k, l] + return get!(colordict, (a.label, b.label, c.label, d.label, e.label, f.label), zero_array) +end + +# interface with TensorKit where necessary +#----------------------------------------- + +# TODO: can remove this once the otimes assert is removed +function TensorKit.fuse(V₁::GradedSpace{I}, Vβ‚‚::GradedSpace{I}) where {I <: BimoduleSector} + dims = TensorKit.SectorDict{I, Int}() + for a in sectors(V₁), b in sectors(Vβ‚‚) + a.j == b.i || continue # skip if not compatible + for c in a βŠ— b + dims[c] = get(dims, c, 0) + Nsymbol(a, b, c) * dim(V₁, a) * dim(Vβ‚‚, b) + end end + return typeof(V₁)(dims) end + +#TODO: these might not be necessary anymore after TensorKit#291 +# check after BlockTensorKit#38 + +# function TensorKit.unitspace(S::SumSpace{<:GradedSpace{<:BimoduleSector}}) +# @assert !isempty(S) "Cannot determine type of empty space" +# return SumSpace(oneunit(first(S.spaces))) # assuming diagonal SumSpace (like in MPSKit) +# end + +# function rightunitspace(S::SumSpace{<:GradedSpace{<:BimoduleSector}}) +# @assert !isempty(S) "Cannot determine type of empty space" +# return SumSpace(rightunitspace(first(S.spaces))) +# end + +# function leftunitspace(S::SumSpace{<:GradedSpace{<:BimoduleSector}}) +# @assert !isempty(S) "Cannot determine type of empty space" +# return SumSpace(leftunitspace(first(S.spaces))) +# end diff --git a/test/TK_compat.jl b/test/TK_compat.jl new file mode 100644 index 0000000..f3e0c98 --- /dev/null +++ b/test/TK_compat.jl @@ -0,0 +1,1452 @@ +## start of TensorKit tests ### + +using TensorKit +using LinearAlgebra: LinearAlgebra +const TK = TensorKit + +println("---------------------------------") +println("| Multifusion space tests |") +println("---------------------------------") + +@timedtestset "Multifusion spaces " verbose = true begin + @timedtestset "GradedSpace: $(TK.type_repr(Vect[I]))" begin + gen = (values(I)[k] => (k + 1) for k in 1:length(values(I))) + + V = GradedSpace(gen) + @test eval(Meta.parse(type_repr(typeof(V)))) == typeof(V) + @test eval_show(V) == V + @test eval_show(V') == V' + @test V' == GradedSpace(gen; dual = true) + @test V == @constinferred GradedSpace(gen...) + @test V' == @constinferred GradedSpace(gen...; dual = true) + @test V == @constinferred GradedSpace(tuple(gen...)) + @test V' == @constinferred GradedSpace(tuple(gen...); dual = true) + @test V == @constinferred GradedSpace(Dict(gen)) + @test V' == @constinferred GradedSpace(Dict(gen); dual = true) + @test V == @inferred Vect[I](gen) + @test V' == @constinferred Vect[I](gen; dual = true) + @test V == @constinferred Vect[I](gen...) + @test V' == @constinferred Vect[I](gen...; dual = true) + @test V == @constinferred Vect[I](Dict(gen)) + @test V' == @constinferred Vect[I](Dict(gen); dual = true) + @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V') + @test V == GradedSpace(reverse(collect(gen))...) + @test eval_show(V) == V + @test eval_show(typeof(V)) == typeof(V) + + @test dim(@constinferred(zerospace(V))) == 0 + + W = @constinferred GradedSpace(unit => 1 for unit in allunits(I)) + dict = Dict(unit => 1 for unit in allunits(I)) + @test W == GradedSpace(dict) + @test W == GradedSpace(push!(dict, randsector(I) => 0)) + @test @constinferred(zerospace(V)) == GradedSpace(unit => 0 for unit in allunits(I)) + randunit = rand(collect(allunits(I))) + @test_throws ArgumentError("Sector $(randunit) appears multiple times") GradedSpace(randunit => 1, randunit => 3) + + @test isunitspace(W) + @test @constinferred(unitspace(V)) == W == unitspace(typeof(V)) + @test_throws ArgumentError leftunitspace(V) + @test_throws ArgumentError rightunitspace(V) + @test eval_show(W) == W + + @test isa(V, VectorSpace) + @test isa(V, ElementarySpace) + @test isa(InnerProductStyle(V), HasInnerProduct) + @test isa(InnerProductStyle(V), EuclideanInnerProduct) + @test isa(V, GradedSpace) + @test isa(V, GradedSpace{I}) + @test @constinferred(dual(V)) == @constinferred(conj(V)) == @constinferred(adjoint(V)) != V + @test @constinferred(field(V)) == β„‚ + @test @constinferred(sectortype(V)) == I + slist = @constinferred sectors(V) + @test @constinferred(hassector(V, first(slist))) + @test @constinferred(dim(V)) == sum(dim(s) * dim(V, s) for s in slist) + @test @constinferred(reduceddim(V)) == sum(dim(V, s) for s in slist) + @constinferred dim(V, first(slist)) + + @test @constinferred(βŠ•(V, zerospace(V))) == V + @test @constinferred(βŠ•(V, V)) == Vect[I](c => 2dim(V, c) for c in sectors(V)) + @test @constinferred(βŠ•(V, V, V, V)) == Vect[I](c => 4dim(V, c) for c in sectors(V)) + @test @constinferred(βŠ•(V, unitspace(V))) == Vect[I](c => isunit(c) + dim(V, c) for c in sectors(V)) + @test @constinferred(fuse(V, unitspace(V))) == V + + @testset "$Istr ($i, $j) spaces" for i in 1:r, j in 1:r #TODO: look at these tests better + # space with a single sector + Wleft = @constinferred Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)) + Wright = @constinferred Vect[I]((j, j, label) => 1 for label in 1:MTK._numlabels(I, j, j)) + WM = @constinferred Vect[I]((i, j, label) => 1 for label in 1:MTK._numlabels(I, i, j)) + WMop = @constinferred Vect[I]((j, i, label) => 1 for label in 1:MTK._numlabels(I, j, i)) + + @test leftunitspace(Wleft) == rightunitspace(Wleft) + @test leftunitspace(Wright) == rightunitspace(Wright) + @test @constinferred(leftunitspace(βŠ•(Wleft, WM))) == leftunitspace(Wleft) + @test @constinferred(leftunitspace(βŠ•(Wright, WMop))) == leftunitspace(Wright) + @test @constinferred(rightunitspace(βŠ•(Wright, WM))) == rightunitspace(Wright) + @test @constinferred(rightunitspace(βŠ•(Wleft, WMop))) == rightunitspace(Wleft) + + if i != j # some tests specialised for modules + + # sensible direct sums and fuses + ul, ur = unit(I(i, i, 1)), unit(I(j, j, 1)) + @test @constinferred(βŠ•(Wleft, WM)) == + Vect[I](c => 1 for c in sectors(V) if leftunit(c) == ul == rightunit(c) || (c.i == i && c.j == j)) + @test @constinferred(βŠ•(Wright, WMop)) == + Vect[I](c => 1 for c in sectors(V) if leftunit(c) == ur == rightunit(c) || (c.i == j && c.j == i)) + @test @constinferred(βŠ•(Wright, WM)) == + Vect[I](c => 1 for c in sectors(V) if rightunit(c) == ur == leftunit(c) || (c.i == i && c.j == j)) + @test @constinferred(βŠ•(Wleft, WMop)) == + Vect[I](c => 1 for c in sectors(V) if rightunit(c) == ul == leftunit(c) || (c.i == j && c.j == i)) + # round needed below because of numerical F-symbols not being integer when they should be + # although this test might be stupid, because I'm assuming integer qdims bc everything's a group or irrep on the diagonal + @test @constinferred(fuse(Wleft, WM)) == Vect[I](c => round(Int, dim(Wleft)) for c in sectors(WM)) # this might be wrong + @test @constinferred(fuse(Wright, WMop)) == Vect[I](c => round(Int, dim(Wright)) for c in sectors(WMop)) # same + + # less sensible fuse + @test @constinferred(fuse(Wleft, WMop)) == fuse(Wright, WM) == + Vect[I](c => 0 for c in sectors(V)) + + for W in [WM, WMop, Wright] + @test infimum(Wleft, W) == Vect[I](c => 0 for c in sectors(V)) + end + else + @test @constinferred(βŠ•(Wleft, Wright)) == + Vect[I](c => 2 for c in sectors(V) if c.i == c.j == i) + @test @constinferred(fuse(Wleft, WMop)) == fuse(Wright, WM) + end + + for W in [Wleft, Wright] + @test @constinferred(βŠ•(W, rightunitspace(W))) == + Vect[I](c => isunit(c) + dim(W, c) for c in sectors(W)) + @test @constinferred(fuse(W, rightunitspace(W))) == W + end + end + + d = Dict{I, Int}() + for a in sectors(V), b in sectors(V) + a.j == b.i || continue # skip if not compatible + for c in a βŠ— b + d[c] = get(d, c, 0) + dim(V, a) * dim(V, b) * Nsymbol(a, b, c) + end + end + @test @constinferred(fuse(V, V)) == GradedSpace(d) + @test @constinferred(flip(V)) == + Vect[I](conj(c) => dim(V, c) for c in sectors(V))' + @test flip(V) β‰… V + @test flip(V) β‰Ύ V + @test flip(V) β‰Ώ V + @test @constinferred(βŠ•(V, V)) == @constinferred supremum(V, βŠ•(V, V)) + @test V == @constinferred infimum(V, βŠ•(V, V)) + @test V β‰Ί βŠ•(V, V) + @test !(V ≻ βŠ•(V, V)) + + randlen = rand(1:length(values(I))) + s = rand(collect(values(I))[randlen:end]) # such that dim(V, s) > randlen + @test infimum(V, GradedSpace(s => randlen)) == GradedSpace(s => randlen) + @test_throws SpaceMismatch (βŠ•(V, V')) + end + + @timedtestset "HomSpace with $(TK.type_repr(Vect[I])) involving ($i, $j)" for i in 1:r, j in 1:r + V1, V2, V3, V4, V5 = ( + Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), + Vect[I]((i, j, label) => 1 for label in 1:MTK._numlabels(I, i, j)), + Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), # same as V1 + Vect[I]((i, j, 1) => 3), + Vect[I]((j, j, label) => 1 for label in 1:MTK._numlabels(I, j, j)), + ) + W = HomSpace(V1 βŠ— V2, V3 βŠ— V4 βŠ— V5) + + @test W == (V3 βŠ— V4 βŠ— V5 β†’ V1 βŠ— V2) + @test W == (V1 βŠ— V2 ← V3 βŠ— V4 βŠ— V5) + @test W' == (V1 βŠ— V2 β†’ V3 βŠ— V4 βŠ— V5) + @test eval(Meta.parse(sprint(show, W))) == W + @test eval(Meta.parse(sprint(show, typeof(W)))) == typeof(W) + @test spacetype(W) == typeof(V1) + @test sectortype(W) == sectortype(V1) + @test W[1] == V1 + @test W[2] == V2 + @test W[3] == V3' + @test W[4] == V4' + @test W[5] == V5' + + @test @constinferred(hash(W)) == hash(deepcopy(W)) != hash(W') + @test W == deepcopy(W) + @test W == @constinferred permute(W, ((1, 2), (3, 4, 5))) + @test permute(W, ((2, 4, 5), (3, 1))) == (V2 βŠ— V4' βŠ— V5' ← V3 βŠ— V1') + @test (V1 βŠ— V2 ← V1 βŠ— V2) == @constinferred TK.compose(W, W') + + @test (V1 βŠ— V2 ← V3 βŠ— V4 βŠ— V5 βŠ— rightunitspace(V5)) == + @constinferred(insertleftunit(W)) == + @constinferred(insertrightunit(W)) + @test @constinferred(removeunit(insertleftunit(W), $(numind(W) + 1))) == W + @test_throws BoundsError insertrightunit(W, 6) + @test_throws BoundsError insertleftunit(W, 0) + + @test (V1 βŠ— V2 βŠ— rightunitspace(V2) ← V3 βŠ— V4 βŠ— V5) == + @constinferred(insertrightunit(W, 2)) + @test (V1 βŠ— V2 ← leftunitspace(V3) βŠ— V3 βŠ— V4 βŠ— V5) == + @constinferred(insertleftunit(W, 3)) + @test @constinferred(removeunit(insertleftunit(W, 3), 3)) == W + @test_throws ArgumentError @constinferred(insertrightunit(one(V1) ← V1, 0)) # should I specify it's the other error? + @test_throws ArgumentError insertleftunit(one(V1) ← V1, 0) + end +end + +println("---------------------------------------") +println("| Multifusion fusion tree tests |") +println("---------------------------------------") + +@timedtestset "Fusion trees for $(TK.type_repr(I)) involving ($i, $j)" verbose = true for i in 1:r, j in 1:r + N = 6 + out = random_fusion(I, i, j, Val(N)) + isdual = ntuple(n -> rand(Bool), N) + in = rand(collect(βŠ—(out...))) # will be in π’žβ±Όβ±Ό with this choice of out + + numtrees = length(fusiontrees(out, in, isdual)) # will be 1 for i != j + @test numtrees == count(n -> true, fusiontrees(out, in, isdual)) + + it = @constinferred fusiontrees(out, in, isdual) + @constinferred Nothing iterate(it) + f, s = iterate(it) + @constinferred Nothing iterate(it, s) + @test f == @constinferred first(it) + @testset "Fusion tree $Istr: printing" begin + @test eval(Meta.parse(sprint(show, f))) == f + end + + C0, D0 = unit(I(i, i, 1)), unit(I(j, j, 1)) + @testset "Fusion tree $Istr: constructor properties" for u in (C0, D0) + @constinferred FusionTree((), u, (), (), ()) + @constinferred FusionTree((u,), u, (false,), (), ()) + @constinferred FusionTree((u, u), u, (false, false), (), (1,)) + @constinferred FusionTree((u, u, u), u, (false, false, false), (u,), (1, 1)) + @constinferred FusionTree( + (u, u, u, u), u, (false, false, false, false), (u, u), (1, 1, 1) + ) + @test_throws MethodError FusionTree((u, u, u), u, (false, false), (u,), (1, 1)) + @test_throws MethodError FusionTree( + (u, u, u), u, (false, false, false), (u, u), (1, 1) + ) + @test_throws MethodError FusionTree( + (u, u, u), u, (false, false, false), (u,), (1, 1, 1) + ) + @test_throws MethodError FusionTree((u, u, u), u, (false, false, false), (), (1,)) + + f = FusionTree((u, u, u), u, (false, false, false), (u,), (1, 1)) + @test sectortype(f) == I + @test length(f) == 3 + @test FusionStyle(f) == FusionStyle(I) + @test BraidingStyle(f) == BraidingStyle(I) + + if FusionStyle(I) isa UniqueFusion + @constinferred FusionTree((), u, ()) + @constinferred FusionTree((u,), u, (false,)) + @constinferred FusionTree((u, u), u, (false, false)) + @constinferred FusionTree((u, u, u), u) + if UnitStyle(I) isa SimpleUnit + @constinferred FusionTree((u, u, u, u)) + else + @test_throws ArgumentError FusionTree((u, u, u, u)) + end + @test_throws MethodError FusionTree((u, u), u, (false, false, false)) + else + @test_throws ArgumentError FusionTree((), u, ()) + @test_throws ArgumentError FusionTree((u,), u, (false,)) + @test_throws ArgumentError FusionTree((u, u), u, (false, false)) + @test_throws ArgumentError FusionTree((u, u, u), u) + if I <: ProductSector && UnitStyle(I) isa GenericUnit + @test_throws DomainError FusionTree((u, u, u, u)) + else + @test_throws ArgumentError FusionTree((u, u, u, u)) + end + end + end + + @testset "Fusion tree $Istr: insertat" begin + N = 4 + out2 = random_fusion(I, i, j, Val(N)) + in2 = rand(collect(βŠ—(out2...))) + isdual2 = ntuple(n -> rand(Bool), N) + f2 = rand(collect(fusiontrees(out2, in2, isdual2))) + for k in 1:N + out1 = random_fusion(I, i, j, Val(N)) # guaranteed good fusion + out1 = Base.setindex(out1, in2, i) # can lead to poor fusion + while isempty(βŠ—(out1...)) # TODO: better way to do this? + out1 = random_fusion(I, i, j, Val(N)) + out1 = Base.setindex(out1, in2, i) + end + in1 = rand(collect(βŠ—(out1...))) + isdual1 = ntuple(n -> rand(Bool), N) + isdual1 = Base.setindex(isdual1, false, k) + f1 = rand(collect(fusiontrees(out1, in1, isdual1))) + + trees = @constinferred TK.insertat(f1, k, f2) + @test norm(values(trees)) β‰ˆ 1 + + f1a, f1b = @constinferred TK.split(f1, $k) + @test length(TK.insertat(f1b, 1, f1a)) == 1 + @test first(TK.insertat(f1b, 1, f1a)) == (f1 => 1) + + # no braid tests for non-hardcoded example + end + end + # no planar trace tests + + @testset "Fusion tree $Istr: elementary artin braid" begin + N = length(out) + isdual = ntuple(n -> rand(Bool), N) + # no general artin braid test + + # not sure how useful this test is, it does the trivial braiding (choice of out) + f = rand(collect(it)) # in this case the 1 tree + d1 = TK.artin_braid(f, 2) # takes unit C0 with current out + d2 = empty(d1) + for (f1, coeff1) in d1 + for (f2, coeff2) in TK.artin_braid(f1, 3) + d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 + end + end + d1 = d2 + d2 = empty(d1) + for (f1, coeff1) in d1 + for (f2, coeff2) in TK.artin_braid(f1, 3; inv = true) + d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 + end + end + d1 = d2 + d2 = empty(d1) + for (f1, coeff1) in d1 + for (f2, coeff2) in TK.artin_braid(f1, 2; inv = true) + d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 + end + end + d1 = d2 + for (f1, coeff1) in d1 + if f1 == f + @test coeff1 β‰ˆ 1 + else + @test isapprox(coeff1, 0; atol = 1.0e-12, rtol = 1.0e-12) + end + end + end + + # no braiding and permuting test + @testset "Fusion tree $Istr: merging" begin + N = 3 + out1 = random_fusion(I, i, j, N) + out2 = random_fusion(I, i, j, N) + in1 = rand(collect(βŠ—(out1...))) + in2 = rand(collect(βŠ—(out2...))) + tp = βŠ—(in1, in2) # messy solution but it works + while isempty(tp) + out1 = random_fusion(I, i, j, Val(N)) + out2 = random_fusion(I, i, j, Val(N)) + in1 = rand(collect(βŠ—(out1...))) + in2 = rand(collect(βŠ—(out2...))) + tp = βŠ—(in1, in2) + end + + f1 = rand(collect(fusiontrees(out1, in1))) + f2 = rand(collect(fusiontrees(out2, in2))) + + + @test dim(in1) * dim(in2) β‰ˆ sum( + abs2(coeff) * dim(c) for c in in1 βŠ— in2 + for ΞΌ in 1:Nsymbol(in1, in2, c) + for (f, coeff) in TK.merge(f1, f2, c, ΞΌ) + ) + # no merge and braid interplay tests + end + + # double fusion tree tests + N = 4 + out = random_fusion(I, i, j, Val(N)) + out2 = random_fusion(I, i, j, Val(N)) + tp = βŠ—(out...) + tp2 = βŠ—(out2...) + while isempty(intersect(tp, tp2)) # guarantee fusion to same coloring + out2 = random_fusion(I, i, j, Val(N)) + tp2 = βŠ—(out2...) + end + @test_throws ArgumentError fusiontrees((out..., map(dual, out)...)) + incoming = rand(collect(intersect(tp, tp2))) + f1 = rand(collect(fusiontrees(out, incoming, ntuple(n -> rand(Bool), N)))) + f2 = rand(collect(fusiontrees(out2, incoming, ntuple(n -> rand(Bool), N)))) # no permuting + + @testset "Double fusion tree $Istr: repartitioning" begin + for n in 0:(2 * N) + d = @constinferred TK.repartition(f1, f2, $n) + @test dim(incoming) β‰ˆ + sum(abs2(coef) * dim(f1.coupled) for ((f1, f2), coef) in d) + d2 = Dict{typeof((f1, f2)), valtype(d)}() + for ((f1β€², f2β€²), coeff) in d + for ((f1β€²β€², f2β€²β€²), coeff2) in TK.repartition(f1β€², f2β€², N) + d2[(f1β€²β€², f2β€²β€²)] = get(d2, (f1β€²β€², f2β€²β€²), zero(coeff)) + coeff2 * coeff + end + end + for ((f1β€², f2β€²), coeff2) in d2 + if f1 == f1β€² && f2 == f2β€² + @test coeff2 β‰ˆ 1 + else + @test isapprox(coeff2, 0; atol = 1.0e-12, rtol = 1.0e-12) + end + end + end + end + + # no double fusion tree permutation tests + + # very slow for (1, 6), (3, 4), (3, 5), (3, 6), (5, 6), (6, 1), (6, 5), (7, 1), (7, 4), (7, 6) + @testset "Double fusion tree $Istr: transposition" begin + for n in 0:(2N) + i0 = rand(1:(2N)) + p = mod1.(i0 .+ (1:(2N)), 2N) + ip = mod1.(-i0 .+ (1:(2N)), 2N) + pβ€² = tuple(getindex.(Ref(vcat(1:N, (2N):-1:(N + 1))), p)...) + p1, p2 = pβ€²[1:n], pβ€²[(2N):-1:(n + 1)] + ipβ€² = tuple(getindex.(Ref(vcat(1:n, (2N):-1:(n + 1))), ip)...) + ip1, ip2 = ipβ€²[1:N], ipβ€²[(2N):-1:(N + 1)] + + d = @constinferred transpose(f1, f2, p1, p2) + @test dim(incoming) β‰ˆ + sum(abs2(coef) * dim(f1.coupled) for ((f1, f2), coef) in d) + d2 = Dict{typeof((f1, f2)), valtype(d)}() + for ((f1β€², f2β€²), coeff) in d + dβ€² = transpose(f1β€², f2β€², ip1, ip2) + for ((f1β€²β€², f2β€²β€²), coeff2) in dβ€² + d2[(f1β€²β€², f2β€²β€²)] = get(d2, (f1β€²β€², f2β€²β€²), zero(coeff)) + coeff2 * coeff + end + end + for ((f1β€², f2β€²), coeff2) in d2 + if f1 == f1β€² && f2 == f2β€² + @test coeff2 β‰ˆ 1 + else + @test abs(coeff2) < 1.0e-12 + end + end + end + end + + @testset "Double fusion tree $Istr: planar trace" begin + d1 = transpose(f1, f1, (N + 1, 1:N..., ((2N):-1:(N + 3))...), (N + 2,)) + f1front, = TK.split(f1, N - 1) + T = sectorscalartype(I) + d2 = Dict{typeof((f1front, f1front)), T}() + for ((f1β€², f2β€²), coeffβ€²) in d1 + for ((f1β€²β€², f2β€²β€²), coeffβ€²β€²) in + TK.planar_trace( + f1β€², f2β€², (2:N...,), (1, ((2N):-1:(N + 3))...), (N + 1,), + (N + 2,) + ) + coeff = coeffβ€² * coeffβ€²β€² + d2[(f1β€²β€², f2β€²β€²)] = get(d2, (f1β€²β€², f2β€²β€²), zero(coeff)) + coeff + end + end + for ((f1_, f2_), coeff) in d2 + if (f1_, f2_) == (f1front, f1front) + @test coeff β‰ˆ dim(f1.coupled) / dim(f1front.coupled) + else + @test abs(coeff) < 1.0e-12 + end + end + end +end + +println("-------------------------------------------") +println("| Multifusion diagonal tensor tests |") +println("-------------------------------------------") + +V = Vect[I](values(I)[k] => 1 for k in 1:length(values(I))) + +@timedtestset "DiagonalTensor" begin + @timedtestset "Basic properties and algebra" begin + for T in (Float32, Float64, ComplexF32, ComplexF64, BigFloat) + # constructors + t = @constinferred DiagonalTensorMap{T}(undef, V) + t = @constinferred DiagonalTensorMap(rand(T, reduceddim(V)), V) + t2 = @constinferred DiagonalTensorMap{T}(undef, space(t)) + @test space(t2) == space(t) + @test_throws ArgumentError DiagonalTensorMap{T}(undef, V^2 ← V) + t2 = @constinferred DiagonalTensorMap{T}(undef, domain(t)) + @test space(t2) == space(t) + @test_throws ArgumentError DiagonalTensorMap{T}(undef, V^2) + # properties + @test @constinferred(hash(t)) == hash(deepcopy(t)) + @test scalartype(t) == T + @test codomain(t) == ProductSpace(V) + @test domain(t) == ProductSpace(V) + @test space(t) == (V ← V) + @test space(t') == (V ← V) + @test dim(t) == dim(space(t)) + # blocks + bs = @constinferred blocks(t) + (c, b1), state = @constinferred Nothing iterate(bs) + @test c == first(blocksectors(V ← V)) + next = @constinferred Nothing iterate(bs, state) + b2 = @constinferred block(t, first(blocksectors(t))) + @test b1 == b2 + @test eltype(bs) === Pair{typeof(c), typeof(b1)} + @test typeof(b1) === TensorKit.blocktype(t) + # basic linear algebra + @test isa(@constinferred(norm(t)), real(T)) + @test norm(t)^2 β‰ˆ dot(t, t) + Ξ± = rand(T) + @test norm(Ξ± * t) β‰ˆ abs(Ξ±) * norm(t) + @test norm(t + t, 2) β‰ˆ 2 * norm(t, 2) + @test norm(t + t, 1) β‰ˆ 2 * norm(t, 1) + @test norm(t + t, Inf) β‰ˆ 2 * norm(t, Inf) + p = 3 * rand(Float64) + @test norm(t + t, p) β‰ˆ 2 * norm(t, p) + @test norm(t) β‰ˆ norm(t') + + @test t == @constinferred(TensorMap(t)) + @test norm(t + TensorMap(t)) β‰ˆ 2 * norm(t) + + @test norm(zerovector!(t)) == 0 + @test norm(one!(t)) β‰ˆ sqrt(dim(V)) + @test one!(t) == id(V) + if T != BigFloat # seems broken for now + @test norm(one!(t) - id(V)) == 0 + end + + t1 = DiagonalTensorMap(rand(T, reduceddim(V)), V) + t2 = DiagonalTensorMap(rand(T, reduceddim(V)), V) + t3 = DiagonalTensorMap(rand(T, reduceddim(V)), V) + Ξ± = rand(T) + Ξ² = rand(T) + @test @constinferred(dot(t1, t2)) β‰ˆ conj(dot(t2, t1)) + @test dot(t2, t1) β‰ˆ conj(dot(t2', t1')) + @test dot(t3, Ξ± * t1 + Ξ² * t2) β‰ˆ Ξ± * dot(t3, t1) + Ξ² * dot(t3, t2) + end + end + + @timedtestset "Basic linear algebra: test via conversion" begin + for T in (Float32, ComplexF64) + t1 = DiagonalTensorMap(rand(T, reduceddim(V)), V) + t2 = DiagonalTensorMap(rand(T, reduceddim(V)), V) + @test norm(t1, 2) β‰ˆ norm(convert(TensorMap, t1), 2) + @test dot(t2, t1) β‰ˆ dot(convert(TensorMap, t2), convert(TensorMap, t1)) + Ξ± = rand(T) + @test convert(TensorMap, Ξ± * t1) β‰ˆ Ξ± * convert(TensorMap, t1) + @test convert(TensorMap, t1') β‰ˆ convert(TensorMap, t1)' + @test convert(TensorMap, t1 + t2) β‰ˆ convert(TensorMap, t1) + convert(TensorMap, t2) + end + end + @timedtestset "Real and imaginary parts" begin + for T in (Float64, ComplexF64, ComplexF32) + t = DiagonalTensorMap(rand(T, reduceddim(V)), V) + + tr = @constinferred real(t) + @test scalartype(tr) <: Real + @test real(convert(TensorMap, t)) == convert(TensorMap, tr) + + ti = @constinferred imag(t) + @test scalartype(ti) <: Real + @test imag(convert(TensorMap, t)) == convert(TensorMap, ti) + + tc = @inferred complex(t) + @test scalartype(tc) <: Complex + @test complex(convert(TensorMap, t)) == convert(TensorMap, tc) + + tc2 = @inferred complex(tr, ti) + @test tc2 β‰ˆ tc + end + end + @timedtestset "Tensor conversion" begin + t = @constinferred DiagonalTensorMap(undef, V) + rand!(t.data) + # element type conversion + tc = complex(t) + @test convert(typeof(tc), t) == tc + @test typeof(convert(typeof(tc), t)) == typeof(tc) + # to and from generic TensorMap + td = DiagonalTensorMap(TensorMap(t)) + @test t == td + @test typeof(td) == typeof(t) + end + @timedtestset "Trace, Multiplication and inverse" begin + t1 = DiagonalTensorMap(rand(Float64, reduceddim(V)), V) + t2 = DiagonalTensorMap(rand(ComplexF64, reduceddim(V)), V) + @test tr(TensorMap(t1)) == @constinferred tr(t1) + @test tr(TensorMap(t2)) == @constinferred tr(t2) + @test TensorMap(@constinferred t1 * t2) β‰ˆ TensorMap(t1) * TensorMap(t2) + @test TensorMap(@constinferred t1 \ t2) β‰ˆ TensorMap(t1) \ TensorMap(t2) + @test TensorMap(@constinferred t1 / t2) β‰ˆ TensorMap(t1) / TensorMap(t2) + @test TensorMap(@constinferred inv(t1)) β‰ˆ inv(TensorMap(t1)) + @test TensorMap(@constinferred pinv(t1)) β‰ˆ pinv(TensorMap(t1)) + @test all( + Base.Fix2(isa, DiagonalTensorMap), (t1 * t2, t1 \ t2, t1 / t2, inv(t1), pinv(t1)) + ) + # no V * V' * V ← V or V^2 ← V tests due to Nsymbol erroring where fusion is forbidden + end + @timedtestset "Tensor contraction " for i in 1:r + W = Vect[I]((i, i, label) => 2 for label in 1:MTK._numlabels(I, i, i)) + + d = DiagonalTensorMap(rand(ComplexF64, reduceddim(W)), W) + t = TensorMap(d) + A = randn(ComplexF64, W βŠ— W' βŠ— W, W) + B = randn(ComplexF64, W βŠ— W' βŠ— W, W βŠ— W') # empty for modules so untested + + @planar E1[-1 -2 -3; -4 -5] := B[-1 -2 -3; 1 -5] * d[1; -4] + @planar E2[-1 -2 -3; -4 -5] := B[-1 -2 -3; 1 -5] * t[1; -4] + @test E1 β‰ˆ E2 + @planar E1[-1 -2 -3; -4 -5] = B[-1 -2 -3; -4 1] * d'[-5; 1] + @planar E2[-1 -2 -3; -4 -5] = B[-1 -2 -3; -4 1] * t'[-5; 1] + @test E1 β‰ˆ E2 + @planar E1[-1 -2 -3; -4 -5] = B[1 -2 -3; -4 -5] * d[-1; 1] + @planar E2[-1 -2 -3; -4 -5] = B[1 -2 -3; -4 -5] * t[-1; 1] + @test E1 β‰ˆ E2 + @planar E1[-1 -2 -3; -4 -5] = B[-1 1 -3; -4 -5] * d[1; -2] + @planar E2[-1 -2 -3; -4 -5] = B[-1 1 -3; -4 -5] * t[1; -2] + @test E1 β‰ˆ E2 + @planar E1[-1 -2 -3; -4 -5] = B[-1 -2 1; -4 -5] * d'[-3; 1] + @planar E2[-1 -2 -3; -4 -5] = B[-1 -2 1; -4 -5] * t'[-3; 1] + @test E1 β‰ˆ E2 + end + @timedtestset "Tensor functions" begin + for T in (Float64, ComplexF64) + d = DiagonalTensorMap(rand(T, reduceddim(V)), V) + # rand is important for positive numbers in the real case, for log and sqrt + t = TensorMap(d) + @test @constinferred exp(d) β‰ˆ exp(t) + @test @constinferred log(d) β‰ˆ log(t) + @test @constinferred sqrt(d) β‰ˆ sqrt(t) + @test @constinferred sin(d) β‰ˆ sin(t) + @test @constinferred cos(d) β‰ˆ cos(t) + @test @constinferred tan(d) β‰ˆ tan(t) + @test @constinferred cot(d) β‰ˆ cot(t) + @test @constinferred sinh(d) β‰ˆ sinh(t) + @test @constinferred cosh(d) β‰ˆ cosh(t) + @test @constinferred tanh(d) β‰ˆ tanh(t) + @test @constinferred coth(d) β‰ˆ coth(t) + @test @constinferred asin(d) β‰ˆ asin(t) + @test @constinferred acos(d) β‰ˆ acos(t) + @test @constinferred atan(d) β‰ˆ atan(t) + @test @constinferred acot(d) β‰ˆ acot(t) + @test @constinferred asinh(d) β‰ˆ asinh(t) + @test @constinferred acosh(one(d) + d) β‰ˆ acosh(one(t) + t) + @test @constinferred atanh(d) β‰ˆ atanh(t) + @test @constinferred acoth(one(t) + d) β‰ˆ acoth(one(d) + t) + end + end +end + +# no conversion tests because no fusion tensor +# no permute tests: NoBraiding() + +println("---------------------------------------") +println("Tensors with symmetry: $Istr") +println("---------------------------------------") + +@timedtestset "Tensors with symmetry involving $Istr ($i, $j)" verbose = true for i in 1:r, j in 1:r + isdiag = i == j + + VC = ( + Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), + Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), # avoids OOMs? + Vect[I](unit(I(i, i, 1)) => 2, rand_object(I, i, i) => 1), + Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), + Vect[I](unit(I(i, i, 1)) => 2, rand_object(I, i, i) => 3), + ) + + VM = Vect[I]((i, j, label) => 1 for label in 1:MTK._numlabels(I, i, j)) # all module objects + + VM1 = ( + Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), # written so V1 βŠ— V2 ← V3 βŠ— V4 βŠ— V5 works + Vect[I](rand_object(I, i, j) => 2), # generally less blocksectors + Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), + VM, # important that V4 is module-graded + Vect[I](unit(I(j, j, 1)) => 2, rand_object(I, j, j) => 1), + ) + + VM2 = ( + Vect[I](rand_object(I, i, j) => 2), # second set where module is V1 here + Vect[I](unit(I(j, j, 1)) => 1, rand_object(I, j, j) => 1), + Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), + VM, + Vect[I](unit(I(j, j, 1)) => 2, rand_object(I, j, j) => 1), + ) + + Vcol = isdiag ? (VC,) : (VM1, VM2) # avoid duplicate runs + + for V in Vcol # TODO: add enumerate to keep track of potential erroring space + V1, V2, V3, V4, V5 = V + @timedtestset "Basic tensor properties" begin + W = isdiag ? V1 βŠ— V2 βŠ— V3 βŠ— V4 βŠ— V5 : V3 βŠ— V4 βŠ— V5 # fusion matters + for T in (Int, Float32, Float64, ComplexF32, ComplexF64, BigFloat) + t = @constinferred zeros(T, W) # empty for i != j b/c blocks are module-graded + @test @constinferred(hash(t)) == hash(deepcopy(t)) + @test scalartype(t) == T + @test norm(t) == 0 + @test codomain(t) == W + @test space(t) == (W ← one(W)) + @test domain(t) == one(W) + @test typeof(t) == TensorMap{T, spacetype(t), length(W), 0, Vector{T}} + # blocks + bs = @constinferred blocks(t) + if !isempty(bs) + (c, b1), state = @constinferred Nothing iterate(bs) # errors if fusion gives empty data + # @test c == first(blocksectors(W)) # unit doesn't have label 1 + next = @constinferred Nothing iterate(bs, state) + b2 = @constinferred block(t, first(blocksectors(t))) + @test b1 == b2 + @test eltype(bs) === Pair{typeof(c), typeof(b1)} + @test typeof(b1) === TK.blocktype(t) + @test typeof(c) === sectortype(t) + end + end + end + + @timedtestset "Tensor Dict conversion" begin + W = V1 βŠ— V2 ← V3 βŠ— V4 βŠ— V5 # rewritten to be compatible with module fusion + for T in (Int, Float32, ComplexF64) + t = @constinferred rand(T, W) + d = convert(Dict, t) + @test t == convert(TensorMap, d) + end + end + + @timedtestset "Basic linear algebra" begin + W = V1 βŠ— V2 ← V3 βŠ— V4 βŠ— V5 + for T in (Float32, ComplexF64) + t = @constinferred rand(T, W) # fusion matters here + @test scalartype(t) == T + @test space(t) == W + @test space(t') == W' + @test dim(t) == dim(space(t)) + @test codomain(t) == codomain(W) + @test domain(t) == domain(W) + # blocks for adjoint + bs = @constinferred blocks(t') + (c, b1), state = @constinferred Nothing iterate(bs) + @test c == first(blocksectors(W')) + next = @constinferred Nothing iterate(bs, state) + b2 = @constinferred block(t', first(blocksectors(t'))) + @test b1 == b2 + @test eltype(bs) === Pair{typeof(c), typeof(b1)} + @test typeof(b1) === TensorKit.blocktype(t') + @test typeof(c) === sectortype(t) + # linear algebra + @test isa(@constinferred(norm(t)), real(T)) + @test norm(t)^2 β‰ˆ dot(t, t) + Ξ± = rand(T) + @test norm(Ξ± * t) β‰ˆ abs(Ξ±) * norm(t) + @test norm(t + t, 2) β‰ˆ 2 * norm(t, 2) + @test norm(t + t, 1) β‰ˆ 2 * norm(t, 1) + @test norm(t + t, Inf) β‰ˆ 2 * norm(t, Inf) + p = 3 * rand(Float64) + @test norm(t + t, p) β‰ˆ 2 * norm(t, p) + @test norm(t) β‰ˆ norm(t') + + t2 = @constinferred rand!(similar(t)) + Ξ² = rand(T) + @test @constinferred(dot(Ξ² * t2, Ξ± * t)) β‰ˆ conj(Ξ²) * Ξ± * conj(dot(t, t2)) + @test dot(t2, t) β‰ˆ conj(dot(t, t2)) + @test dot(t2, t) β‰ˆ conj(dot(t2', t')) + @test dot(t2, t) β‰ˆ dot(t', t2') + + if !isempty(blocksectors(V2 βŠ— V1)) + i1 = @constinferred(isomorphism(T, V1 βŠ— V2, V2 βŠ— V1)) # can't reverse fusion here when modules are involved + i2 = @constinferred(isomorphism(Vector{T}, V2 βŠ— V1, V1 βŠ— V2)) + @test i1 * i2 == @constinferred(id(T, V1 βŠ— V2)) + @test i2 * i1 == @constinferred(id(Vector{T}, V2 βŠ— V1)) + end + + w = @constinferred isometry(T, V1 βŠ— (rightunitspace(V1) βŠ• rightunitspace(V1)), V1) + @test dim(w) == 2 * dim(V1 ← V1) + @test w' * w == id(Vector{T}, V1) + @test w * w' == (w * w')^2 + end + end + + @timedtestset "Trivial space insertion and removal" begin + W = V1 βŠ— V2 ← V3 βŠ— V4 βŠ— V5 + for T in (Float32, ComplexF64) + t = @constinferred rand(T, W) # fusion matters here + t2 = @constinferred insertleftunit(t) + @test t2 == @constinferred insertrightunit(t) + @test space(t2) == insertleftunit(space(t)) + @test @constinferred(removeunit(t2, $(numind(t2)))) == t + t3 = @constinferred insertleftunit(t; copy = true) + @test t3 == @constinferred insertrightunit(t; copy = true) + @test @constinferred(removeunit(t3, $(numind(t3)))) == t + + @test numind(t2) == numind(t) + 1 + @test scalartype(t2) === T + @test t.data === t2.data + + @test t.data !== t3.data + for (c, b) in blocks(t) + @test b == block(t3, c) + end + + t4 = @constinferred insertrightunit(t, 3; dual = true) + @test numin(t4) == numin(t) + 1 && numout(t4) == numout(t) + for (c, b) in blocks(t) + @test b == block(t4, c) + end + @test @constinferred(removeunit(t4, 4)) == t + + t5 = @constinferred insertleftunit(t, 4; dual = true) + @test numin(t5) == numin(t) + 1 && numout(t5) == numout(t) + for (c, b) in blocks(t) + @test b == block(t5, c) + end + @test @constinferred(removeunit(t5, 4)) == t + end + end + + @timedtestset "Tensor conversion" begin + W = V1 βŠ— V2 + t = @constinferred randn(W ← W) # fusion matters here + @test typeof(convert(TensorMap, t')) == typeof(t) + tc = complex(t) + @test convert(typeof(tc), t) == tc + @test typeof(convert(typeof(tc), t)) == typeof(tc) + @test typeof(convert(typeof(tc), t')) == typeof(tc) + @test Base.promote_typeof(t, tc) == typeof(tc) + @test Base.promote_typeof(tc, t) == typeof(tc + t) + end + + @timedtestset "Full trace: test self-consistency" begin + t = rand(ComplexF64, V1 βŠ— V2 ← V1 βŠ— V2) # avoid permutes + ss = @constinferred tr(t) + @test conj(ss) β‰ˆ tr(t') + @planar s2 = t[a b; a b] + @planar t3[a; b] := t[a c; b c] + @planar s3 = t3[a; a] + + @test ss β‰ˆ s2 + @test ss β‰ˆ s3 + end + + @timedtestset "Partial trace: test self-consistency" begin + t = rand(ComplexF64, V3 βŠ— V4 βŠ— V5 ← V3 βŠ— V4 βŠ— V5) # compatible with module fusion + @planar t2[a; b] := t[c a d; c b d] + @planar t4[a b; c d] := t[e a b; e c d] + @planar t5[a; b] := t4[a c; b c] + @test t2 β‰ˆ t5 + end + + @timedtestset "Trace and contraction" begin #TODO: find some version of this that works for off-diagonal case + t1 = rand(ComplexF64, V3 βŠ— V4 βŠ— V5) + t2 = rand(ComplexF64, V3 βŠ— V4 βŠ— V5) + t3 = t1 βŠ— t2' + # if all(a.i != a.j for a in blocksectors(t3)) + # replace!(x -> rand(ComplexF64), t3.data) # otherwise full of zeros in off-diagonal case + # end + if all(a.i == a.j for a in blocksectors(t3)) + @planar ta[b; a] := conj(t2[x, a, y]) * t1[x, b, y] # works for diagonal case + @planar tb[a; b] := t3[x a y; x b y] + @test ta β‰ˆ tb + end + end + + @timedtestset "Multiplication of isometries: test properties" begin + W2 = V4 βŠ— V5 + W1 = W2 βŠ— (rightunitspace(V5) βŠ• rightunitspace(V5)) + for T in (Float64, ComplexF64) + t1 = randisometry(T, W1, W2) + t2 = randisometry(T, W2 ← W2) + @test isisometric(t1) + @test isunitary(t2) + P = t1 * t1' + @test P * P β‰ˆ P + end + end + + @timedtestset "Multiplication and inverse: test compatibility" begin + W1 = V1 βŠ— V2 + W2 = V3 βŠ— V4 βŠ— V5 + for T in (Float64, ComplexF64) + t1 = rand(T, W1, W1) + t2 = rand(T, W2 ← W2) + t = rand(T, W1, W2) + @test t1 * (t1 \ t) β‰ˆ t + @test (t / t2) * t2 β‰ˆ t + @test t1 \ one(t1) β‰ˆ inv(t1) + @test one(t1) / t1 β‰ˆ pinv(t1) + @test_throws SpaceMismatch inv(t) + @test_throws SpaceMismatch t2 \ t + @test_throws SpaceMismatch t / t1 + tp = pinv(t) * t + @test tp β‰ˆ tp * tp + end + end + + @timedtestset "diag/diagm" begin + W = V1 βŠ— V2 ← V3 βŠ— V4 βŠ— V5 + t = randn(ComplexF64, W) + d = LinearAlgebra.diag(t) + D = LinearAlgebra.diagm(codomain(t), domain(t), d) + @test LinearAlgebra.isdiag(D) + @test LinearAlgebra.diag(D) == d + end + + @timedtestset "Sylvester equation" begin + for T in (Float32, ComplexF64) + tA = rand(T, V1 βŠ— V2, V1 βŠ— V2) # rewritten for modules + tB = rand(T, V4 βŠ— V5, V4 βŠ— V5) + tA = 3 // 2 * leftorth(tA; alg = TK.Polar())[1] + tB = 1 // 5 * leftorth(tB; alg = TK.Polar())[1] + tC = rand(T, V1 βŠ— V2, V4 βŠ— V5) + t = @constinferred sylvester(tA, tB, tC) + @test codomain(t) == V1 βŠ— V2 + @test domain(t) == V4 βŠ— V5 + @test norm(tA * t + t * tB + tC) < + (norm(tA) + norm(tB) + norm(tC)) * eps(real(T))^(2 / 3) + end + end + + @timedtestset "Tensor product: test via norm preservation" begin # OOMs over here with full spaces + for T in (Float32, ComplexF64) + if !isempty(blocksectors(V2 βŠ— V1)) + t1 = rand(T, V2 βŠ— V3 βŠ— V1, V1 βŠ— V2) + t2 = rand(T, V2 βŠ— V1 βŠ— V3, V1 βŠ— V1) + else + t1 = rand(T, V3 βŠ— V4 βŠ— V5, V1 βŠ— V2) + t2 = rand(T, V5' βŠ— V4' βŠ— V3', V2' βŠ— V1') + end + t = @constinferred (t1 βŠ— t2) + @test norm(t) β‰ˆ norm(t1) * norm(t2) + end + end + + # TODO: should this test exist? + @timedtestset "Tensor product: test via tensor contraction" begin + # W = V3 βŠ— V4 βŠ— V5 ← V1 βŠ— V2 + W = V4 ← V1 βŠ— V2 # less costly + for T in (Float32, ComplexF64) + if !isdiag + t1 = rand(T, W) + t2 = rand(T, V4' ← V2' βŠ— V1') + # t2 = rand(T, V5' βŠ— V4' βŠ— V3', V2' βŠ— V1') # same as previous test + # @planar tβ€²[1 2 3 6 7 8; 4 5 9 10] := t1[1 2 3; 4 5] * t2[6 7 8; 9 10] + @planar tβ€²[1 4; 2 3 5 6] := t1[1; 2 3] * t2[4; 5 6] + else + t1 = rand(T, V2 βŠ— V3, V1) + t2 = rand(T, V2, V1 βŠ— V3) + @planar tβ€²[1 2 4; 3 5 6] := t1[1 2; 3] * t2[4; 5 6] + end + t = @constinferred (t1 βŠ— t2) + @test t β‰ˆ tβ€² + end + end + end + + @timedtestset "Tensor absorption" begin + # absorbing small into large + if !isempty(blocksectors(V2 βŠ— V3)) + t1 = zeros(V1 βŠ• V1, V2 βŠ— V3) + t2 = rand(V1, V2 βŠ— V3) + else + t1 = zeros(V1 βŠ• V2, V3 βŠ— V4 βŠ— V5) + t2 = rand(V1, V3 βŠ— V4 βŠ— V5) + end + t3 = @constinferred absorb(t1, t2) + @test norm(t3) β‰ˆ norm(t2) + @test norm(t1) == 0 + t4 = @constinferred absorb!(t1, t2) + @test t1 === t4 + @test t3 β‰ˆ t4 + + # absorbing large into small + if !isempty(blocksectors(V2 βŠ— V3)) + t1 = rand(V1 βŠ• V1, V2 βŠ— V3) + t2 = zeros(V1, V2 βŠ— V3) + else + t1 = rand(V1 βŠ• V2, V3 βŠ— V4 βŠ— V5) + t2 = zeros(V1, V3 βŠ— V4 βŠ— V5) + end + t3 = @constinferred absorb(t2, t1) + @test norm(t3) < norm(t1) + @test norm(t2) == 0 + t4 = @constinferred absorb!(t2, t1) + @test t2 === t4 + @test t3 β‰ˆ t4 + end +end + +println("---------------------------------------") +println("Factorizations with symmetry: $Istr") +println("---------------------------------------") + +@timedtestset "Factorizations with symmetry involving $Istr ($i, $j)" verbose = true for i in 1:r, j in 1:r + isdiag = i == j + + VC = ( + Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), + Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), # avoids OOMs? + Vect[I](unit(I(i, i, 1)) => 2, rand_object(I, i, i) => 1), + Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), + Vect[I](unit(I(i, i, 1)) => 2, rand_object(I, i, i) => 3), + ) + + VM = Vect[I]((i, j, label) => 1 for label in 1:MTK._numlabels(I, i, j)) # all module objects + + VM1 = ( + Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), # written so V1 βŠ— V2 ← V3 βŠ— V4 βŠ— V5 works + Vect[I](rand_object(I, i, j) => 2), # generally less blocksectors + Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), + VM, # important that V4 is module-graded + Vect[I](unit(I(j, j, 1)) => 2, rand_object(I, j, j) => 1), + ) + + VM2 = ( + Vect[I](rand_object(I, i, j) => 2), # second set where module is V1 here + Vect[I](unit(I(j, j, 1)) => 1, rand_object(I, j, j) => 1), + Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), + VM, + Vect[I](unit(I(j, j, 1)) => 2, rand_object(I, j, j) => 1), + ) + + Vs = isdiag ? (VC,) : (VM1, VM2) # avoid duplicate runs + + # some fail for (2, 2), (3, 3), (6, 6) + # rightorth RQ(pos) and Polar (fail) for 2nd space + # leftorth with QL(pos) and Polar for 1st space + # leftnull QR for 1st space + # cond and rank leftnull for 1st space + + # factorization tests require equal objects in blocksectors of domain and codomain, so just put them all + # FIXME: not sure if still needed + # VC_all = fill(Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), 5) + + # VM1_all = (Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), + # VM, + # Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), + # VM, + # Vect[I](unit(I(j, j, 1)) => 1, rand_object(I, j, j) => 2) + # ) + + # VM2_all = (VM, + # Vect[I](unit(I(j, j, 1)) => 1, rand_object(I, j, j) => 1), + # Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), + # VM, + # Vect[I](unit(I(j, j, 1)) => 2, rand_object(I, j, j) => 2) + # ) + + # fact_Vs = (i != j) ? (VM1_all, VM2_all) : (VC_all,) + + @timedtestset "Factorization" for V in Vs + V1, V2, V3, V4, V5 = V + W = V1 βŠ— V2 + @assert !isempty(blocksectors(W)) + @assert !isempty(intersect(blocksectors(V4), blocksectors(W))) + + @testset "QR decomposition" begin + for T in eltypes, + t in ( + rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)', + DiagonalTensorMap(rand(T, reduceddim(V1)), V1), + ) + + Q, R = @constinferred qr_full(t) + @test Q * R β‰ˆ t + @test isunitary(Q) + + Q, R = @constinferred qr_compact(t) + @test Q * R β‰ˆ t + @test isisometric(Q) + + Q, R = @constinferred left_orth(t) + @test Q * R β‰ˆ t + @test isisometric(Q) + + N = @constinferred qr_null(t) + @test isisometric(N) + @test norm(N' * t) β‰ˆ 0 atol = 100 * eps(norm(t)) + + N = @constinferred left_null(t) + @test isisometric(N) + @test norm(N' * t) β‰ˆ 0 atol = 100 * eps(norm(t)) + end + + # empty tensor + for T in eltypes + t = rand(T, V1 βŠ— V2, zerospace(V1)) + + Q, R = @constinferred qr_full(t) + @test Q * R β‰ˆ t + @test isunitary(Q) + @test dim(R) == dim(t) == 0 + + Q, R = @constinferred qr_compact(t) + @test Q * R β‰ˆ t + @test isisometric(Q) + @test dim(Q) == dim(R) == dim(t) + + Q, R = @constinferred left_orth(t) + @test Q * R β‰ˆ t + @test isisometric(Q) + @test dim(Q) == dim(R) == dim(t) + + N = @constinferred qr_null(t) + @test isunitary(N) + @test norm(N' * t) β‰ˆ 0 atol = 100 * eps(norm(t)) + end + end + + @testset "LQ decomposition" begin + for T in eltypes, + t in ( + rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)', + DiagonalTensorMap(rand(T, reduceddim(V1)), V1), + ) + + L, Q = @constinferred lq_full(t) + @test L * Q β‰ˆ t + @test isunitary(Q) + + L, Q = @constinferred lq_compact(t) + @test L * Q β‰ˆ t + @test isisometric(Q; side = :right) + + L, Q = @constinferred right_orth(t) + @test L * Q β‰ˆ t + @test isisometric(Q; side = :right) + + Nα΄΄ = @constinferred lq_null(t) + @test isisometric(Nα΄΄; side = :right) + @test norm(t * Nα΄΄') β‰ˆ 0 atol = 100 * eps(norm(t)) + end + + for T in eltypes + # empty tensor + t = rand(T, zerospace(V1), V1 βŠ— V2) + + L, Q = @constinferred lq_full(t) + @test L * Q β‰ˆ t + @test isunitary(Q) + @test dim(L) == dim(t) == 0 + + L, Q = @constinferred lq_compact(t) + @test L * Q β‰ˆ t + @test isisometric(Q; side = :right) + @test dim(Q) == dim(L) == dim(t) + + L, Q = @constinferred right_orth(t) + @test L * Q β‰ˆ t + @test isisometric(Q; side = :right) + @test dim(Q) == dim(L) == dim(t) + + Nα΄΄ = @constinferred lq_null(t) + @test isunitary(Nα΄΄) + @test norm(t * Nα΄΄') β‰ˆ 0 atol = 100 * eps(norm(t)) + end + end + + @testset "Polar decomposition" begin + for T in eltypes, + t in ( + rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)', + DiagonalTensorMap(rand(T, reduceddim(V1)), V1), + ) + + @assert domain(t) β‰Ύ codomain(t) + w, p = @constinferred left_polar(t) + @test w * p β‰ˆ t + @test isisometric(w) + @test isposdef(p) + + w, p = @constinferred left_orth(t; alg = :polar) + @test w * p β‰ˆ t + @test isisometric(w) + end + + for T in eltypes, + t in (rand(T, W, W), rand(T, W, W)', rand(T, V4, W), rand(T, W, V4)') + + @assert codomain(t) β‰Ύ domain(t) + p, wα΄΄ = @constinferred right_polar(t) + @test p * wα΄΄ β‰ˆ t + @test isisometric(wα΄΄; side = :right) + @test isposdef(p) + + p, wα΄΄ = @constinferred right_orth(t; alg = :polar) + @test p * wα΄΄ β‰ˆ t + @test isisometric(wα΄΄; side = :right) + end + end + + @testset "SVD" begin + for T in eltypes, + t in ( + rand(T, W, W), rand(T, W, W)', + rand(T, W, V4), rand(T, V4, W), + rand(T, W, V4)', rand(T, V4, W)', + DiagonalTensorMap(rand(T, reduceddim(V1)), V1), + ) + + u, s, vα΄΄ = @constinferred svd_full(t) + @test u * s * vα΄΄ β‰ˆ t + @test isunitary(u) + @test isunitary(vα΄΄) + + u, s, vα΄΄ = @constinferred svd_compact(t) + @test u * s * vα΄΄ β‰ˆ t + @test isisometric(u) + @test isposdef(s) + @test isisometric(vα΄΄; side = :right) + + sβ€² = @constinferred svd_vals(t) + @test sβ€² β‰ˆ diagview(s) + @test sβ€² isa TensorKit.SectorVector + + v, c = @constinferred left_orth(t; alg = :svd) + @test v * c β‰ˆ t + @test isisometric(v) + + c, vα΄΄ = @constinferred right_orth(t; alg = :svd) + @test c * vα΄΄ β‰ˆ t + @test isisometric(vα΄΄; side = :right) + + N = @constinferred left_null(t; alg = :svd) + @test isisometric(N) + @test norm(N' * t) β‰ˆ 0 atol = 100 * eps(norm(t)) + + N = @constinferred left_null(t; trunc = (; atol = 100 * eps(norm(t)))) + @test isisometric(N) + @test norm(N' * t) β‰ˆ 0 atol = 100 * eps(norm(t)) + + Nα΄΄ = @constinferred right_null(t; alg = :svd) + @test isisometric(Nα΄΄; side = :right) + @test norm(t * Nα΄΄') β‰ˆ 0 atol = 100 * eps(norm(t)) + + Nα΄΄ = @constinferred right_null(t; trunc = (; atol = 100 * eps(norm(t)))) + @test isisometric(Nα΄΄; side = :right) + @test norm(t * Nα΄΄') β‰ˆ 0 atol = 100 * eps(norm(t)) + end + + # empty tensor + for T in eltypes, t in (rand(T, W, zerospace(V1)), rand(T, zerospace(V1), W)) + U, S, Vα΄΄ = @constinferred svd_full(t) + @test U * S * Vα΄΄ β‰ˆ t + @test isunitary(U) + @test isunitary(Vα΄΄) + + U, S, Vα΄΄ = @constinferred svd_compact(t) + @test U * S * Vα΄΄ β‰ˆ t + @test dim(U) == dim(S) == dim(Vα΄΄) == dim(t) == 0 + end + end + + @testset "truncated SVD" begin + for T in eltypes, + t in ( + randn(T, W, W), randn(T, W, W)', + randn(T, W, V4), randn(T, V4, W), + randn(T, W, V4)', randn(T, V4, W)', + DiagonalTensorMap(randn(T, reduceddim(V1)), V1), + ) + + @constinferred normalize!(t) + + U, S, Vα΄΄, Ο΅ = @constinferred svd_trunc(t; trunc = notrunc()) + @test U * S * Vα΄΄ β‰ˆ t + @test Ο΅ β‰ˆ 0 + @test isisometric(U) + @test isisometric(Vα΄΄; side = :right) + + # dimension of S is a float for IsingBimodule + nvals = round(Int, dim(domain(S)) / 2) + trunc = truncrank(nvals) + U1, S1, Vα΄΄1, Ο΅1 = @constinferred svd_trunc(t; trunc) + @test t * Vα΄΄1' β‰ˆ U1 * S1 + @test isisometric(U1) + @test isisometric(Vα΄΄1; side = :right) + @test norm(t - U1 * S1 * Vα΄΄1) β‰ˆ Ο΅1 atol = eps(real(T))^(4 / 5) + @test dim(domain(S1)) <= nvals + + Ξ» = minimum(diagview(S1)) + trunc = trunctol(; atol = Ξ» - 10eps(Ξ»)) + U2, S2, Vα΄΄2, Ο΅2 = @constinferred svd_trunc(t; trunc) + @test t * Vα΄΄2' β‰ˆ U2 * S2 + @test isisometric(U2) + @test isisometric(Vα΄΄2; side = :right) + @test norm(t - U2 * S2 * Vα΄΄2) β‰ˆ Ο΅2 atol = eps(real(T))^(4 / 5) + @test minimum(diagview(S1)) >= Ξ» + @test U2 β‰ˆ U1 + @test S2 β‰ˆ S1 + @test Vα΄΄2 β‰ˆ Vα΄΄1 + @test Ο΅1 β‰ˆ Ο΅2 + + trunc = truncspace(space(S2, 1)) + U3, S3, Vα΄΄3, Ο΅3 = @constinferred svd_trunc(t; trunc) + @test t * Vα΄΄3' β‰ˆ U3 * S3 + @test isisometric(U3) + @test isisometric(Vα΄΄3; side = :right) + @test norm(t - U3 * S3 * Vα΄΄3) β‰ˆ Ο΅3 atol = eps(real(T))^(4 / 5) + @test space(S3, 1) β‰Ύ space(S2, 1) + + for trunc in (truncerror(; atol = Ο΅2), truncerror(; rtol = Ο΅2 / norm(t))) + U4, S4, Vα΄΄4, Ο΅4 = @constinferred svd_trunc(t; trunc) + @test t * Vα΄΄4' β‰ˆ U4 * S4 + @test isisometric(U4) + @test isisometric(Vα΄΄4; side = :right) + @test norm(t - U4 * S4 * Vα΄΄4) β‰ˆ Ο΅4 atol = eps(real(T))^(4 / 5) + @test Ο΅4 ≀ Ο΅2 + end + + trunc = truncrank(nvals) & trunctol(; atol = Ξ» - 10eps(Ξ»)) + U5, S5, Vα΄΄5, Ο΅5 = @constinferred svd_trunc(t; trunc) + @test t * Vα΄΄5' β‰ˆ U5 * S5 + @test isisometric(U5) + @test isisometric(Vα΄΄5; side = :right) + @test norm(t - U5 * S5 * Vα΄΄5) β‰ˆ Ο΅5 atol = eps(real(T))^(4 / 5) + @test minimum(diagview(S5)) >= Ξ» + @test dim(domain(S5)) ≀ nvals + end + end + + @testset "Eigenvalue decomposition" begin + for T in eltypes, + t in ( + rand(T, V1, V1), rand(T, W, W), rand(T, W, W)', + DiagonalTensorMap(rand(T, reduceddim(V1)), V1), + ) + + d, v = @constinferred eig_full(t) + @test t * v β‰ˆ v * d + + dβ€² = @constinferred eig_vals(t) + @test dβ€² β‰ˆ diagview(d) + @test dβ€² isa TensorKit.SectorVector + + vdv = project_hermitian!(v' * v) + @test @constinferred isposdef(vdv) + t isa DiagonalTensorMap || @test !isposdef(t) # unlikely for non-hermitian map + + nvals = round(Int, dim(domain(t)) / 2) + d, v = @constinferred eig_trunc(t; trunc = truncrank(nvals)) + @test t * v β‰ˆ v * d + @test dim(domain(d)) ≀ nvals + + t2 = @constinferred project_hermitian(t) + D, V = eigen(t2) + @test isisometric(V) + DΜƒ, VΜƒ = @constinferred eigh_full(t2) + @test D β‰ˆ DΜƒ + @test V β‰ˆ VΜƒ + Ξ» = minimum(real, diagview(D)) + @test cond(VΜƒ) β‰ˆ one(real(T)) + @test isposdef(t2) == isposdef(Ξ») + @test isposdef(t2 - Ξ» * one(t2) + 0.1 * one(t2)) + @test !isposdef(t2 - Ξ» * one(t2) - 0.1 * one(t2)) + + d, v = @constinferred eigh_full(t2) + @test t2 * v β‰ˆ v * d + @test isunitary(v) + + dβ€² = @constinferred eigh_vals(t2) + @test dβ€² β‰ˆ diagview(d) + @test dβ€² isa TensorKit.SectorVector + + Ξ» = minimum(real, diagview(d)) + @test cond(v) β‰ˆ one(real(T)) + @test isposdef(t2) == isposdef(Ξ») + @test isposdef(t2 - Ξ» * one(t) + 0.1 * one(t2)) + @test !isposdef(t2 - Ξ» * one(t) - 0.1 * one(t2)) + + d, v = @constinferred eigh_trunc(t2; trunc = truncrank(nvals)) + @test t2 * v β‰ˆ v * d + @test dim(domain(d)) ≀ nvals + end + end + + @testset "Condition number and rank" begin + for T in eltypes, + t in ( + rand(T, W, W), rand(T, W, W)', + rand(T, W, V4), rand(T, V4, W), + rand(T, W, V4)', rand(T, V4, W)', + DiagonalTensorMap(rand(T, reduceddim(V1)), V1), + ) + + d1, d2 = dim(codomain(t)), dim(domain(t)) + r = rank(t) + @test r == min(d1, d2) + @test typeof(r) == typeof(d1) + M = left_null(t) + @test @constinferred(rank(M)) + r β‰ˆ d1 + Mα΄΄ = right_null(t) + @test rank(Mα΄΄) + r β‰ˆ d2 + end + for T in eltypes + u = unitary(T, V1 βŠ— V2, V1 βŠ— V2) + @test @constinferred(cond(u)) β‰ˆ one(real(T)) + @test @constinferred(rank(u)) == dim(V1 βŠ— V2) + + t = rand(T, zerospace(V1), W) + @test rank(t) == 0 + t2 = rand(T, zerospace(V1) * zerospace(V2), zerospace(V1) * zerospace(V2)) + @test rank(t2) == 0 + @test cond(t2) == 0.0 + end + for T in eltypes, t in (rand(T, W, W), rand(T, W, W)') + project_hermitian!(t) + vals = @constinferred LinearAlgebra.eigvals(t) + Ξ»max = maximum(s -> maximum(abs, s), values(vals)) + Ξ»min = minimum(s -> minimum(abs, s), values(vals)) + @test cond(t) β‰ˆ Ξ»max / Ξ»min + end + end + + @testset "Hermitian projections" begin + for T in eltypes, + t in ( + rand(T, V1, V1), rand(T, W, W), rand(T, W, W)', + DiagonalTensorMap(rand(T, reduceddim(V1)), V1), + ) + normalize!(t) + noisefactor = eps(real(T))^(3 / 4) + + th = (t + t') / 2 + ta = (t - t') / 2 + tc = copy(t) + + thβ€² = @constinferred project_hermitian(t) + @test ishermitian(thβ€²) + @test thβ€² β‰ˆ th + @test t == tc + th_approx = th + noisefactor * ta + @test !ishermitian(th_approx) || (T <: Real && t isa DiagonalTensorMap) + @test ishermitian(th_approx; atol = 10 * noisefactor) + + taβ€² = project_antihermitian(t) + @test isantihermitian(taβ€²) + @test taβ€² β‰ˆ ta + @test t == tc + ta_approx = ta + noisefactor * th + @test !isantihermitian(ta_approx) + @test isantihermitian(ta_approx; atol = 10 * noisefactor) || (T <: Real && t isa DiagonalTensorMap) + end + end + + @testset "Isometric projections" begin + for T in eltypes, + t in ( + randn(T, W, W), randn(T, W, W)', + randn(T, W, V4), randn(T, V4, W)', + ) + t2 = project_isometric(t) + @test isisometric(t2) + t3 = project_isometric(t2) + @test t3 β‰ˆ t2 # stability of the projection + @test t2 * (t2' * t) β‰ˆ t + + tc = similar(t) + t3 = @constinferred project_isometric!(copy!(tc, t), t2) + @test t3 === t2 + @test isisometric(t2) + + # test that t2 is closer to A then any other isometry + for k in 1:10 + Ξ΄t = randn!(similar(t)) + t3 = project_isometric(t + Ξ΄t / 100) + @test norm(t - t3) > norm(t - t2) + end + end + end + end +end diff --git a/test/test_A4.jl b/test/test_A4.jl index 696ca5b..59fa329 100644 --- a/test/test_A4.jl +++ b/test/test_A4.jl @@ -1,54 +1,163 @@ -using MultiTensorKit using TensorKitSectors + +testsuite_path = joinpath( + dirname(dirname(pathof(TensorKitSectors))), # TensorKitSectors root + "test", "testsuite.jl" +) +include(testsuite_path) +using .SectorTestSuite + +using MultiTensorKit using Test, TestExtras +const MTK = MultiTensorKit + I = A4Object +Istr = TensorKitSectors.type_repr(I) +r = size(I) -@testset "Basic type properties" begin - Istr = TensorKitSectors.type_repr(I) - @test eval(Meta.parse(sprint(show, I))) == I - @test eval(Meta.parse(TensorKitSectors.type_repr(I))) == I +println("----------------------") +println("| Sector tests |") +println("----------------------") +@testset "Sector test suite" verbose = true begin + @time SectorTestSuite.test_sector(I) end -@testset "Fusion Category $i" for i in 1:12 - objects = A4Object.(i, i, MultiTensorKit._get_dual_cache(I)[i][2]) - - @testset "Basic properties" begin - s = rand(objects, 3) - @test eval(Meta.parse(sprint(show, s[1]))) == s[1] - @test @constinferred(hash(s[1])) == hash(deepcopy(s[1])) - @test isone(@constinferred(one(s[1]))) - @constinferred dual(s[1]) - @constinferred dim(s[1]) - @constinferred frobeniusschur(s[1]) - @constinferred Bsymbol(s...) - @constinferred Fsymbol(s..., s...) +@testset "$Istr ($i, $j) basic properties" for i in 1:r, j in 1:r + Cii_obs = I.(i, i, MTK._get_dual_cache(I)[2][i, i]) + Cij_obs = I.(i, j, MTK._get_dual_cache(I)[2][i, j]) + Cji_obs = I.(j, i, MTK._get_dual_cache(I)[2][j, i]) + Cjj_obs = I.(j, j, MTK._get_dual_cache(I)[2][j, j]) + c, m, mop, d = rand(Cii_obs), rand(Cij_obs), rand(Cji_obs), rand(Cjj_obs) + + if i == j + @testset "Basic fusion properties" begin + @test isunit(@testinferred(unit(c))) + u = I.(i, i, MTK._get_dual_cache(I)[1][i]) + @test u == @testinferred(leftunit(u)) == @testinferred(rightunit(u)) == + @testinferred(unit(u)) + @test isunit(@testinferred(unit(c))) + @test dual(c) == I.(i, i, MTK._get_dual_cache(I)[2][i, i][c.label]) + end + else + @testset "Basic module properties" begin + @test isunit(m) == false + @test isunit(mop) == false + @test (isunit(@testinferred(leftunit(m))) && isunit(@testinferred(rightunit(m)))) + @test unit(c) == leftunit(m) == rightunit(mop) + @test unit(d) == rightunit(m) == leftunit(mop) + @test_throws DomainError unit(m) + @test_throws DomainError unit(mop) + @test dual(m) == I.(j, i, MTK._get_dual_cache(I)[2][i, j][m.label]) + end + + @testset "$Istr Fusion rules" begin + argerr = ArgumentError("invalid fusion channel") + # forbidden fusions + for obs in [(c, d), (d, c), (m, m), (mop, mop), (d, m), (m, c), (mop, d), (c, mop)] + @test isempty(βŠ—(obs...)) + @test_throws argerr Nsymbol(obs..., rand([c, m, mop, d])) + end + + # allowed fusions + for obs in [(c, c), (d, d), (m, mop), (mop, m), (c, m), (mop, c), (m, d), (d, mop)] + @test !isempty(βŠ—(obs...)) + end + + @test Nsymbol(c, unit(c), c) == Nsymbol(d, unit(d), d) == 1 + + @test_throws argerr Nsymbol(m, mop, d) + @test_throws argerr Nsymbol(mop, m, c) + @test_throws argerr Fsymbol(m, mop, m, mop, c, d) + end end +end + +println("-----------------------------") +println("| F-symbol data tests |") +println("-----------------------------") - @testset "Unitarity of F-move" begin - for a in objects, b in objects, c in objects - for d in βŠ—(a, b, c) - es = collect(intersect(βŠ—(a, b), map(dual, βŠ—(c, dual(d))))) - fs = collect(intersect(βŠ—(b, c), map(dual, βŠ—(dual(d), a)))) - Fblocks = Vector{Any}() - for e in es - for f in fs - Fs = Fsymbol(a, b, c, d, e, f) - push!(Fblocks, - reshape(Fs, - (size(Fs, 1) * size(Fs, 2), - size(Fs, 3) * size(Fs, 4)))) - end +# explicitly test everything related to F-symbols +# other option is to edit smallset to sample more +for i in 1:r, j in 1:r + @testset "Unitarity of $Istr F-move ($i, $j)" begin + if i == j + @testset "Unitarity of fusion F-move ($i, $j)" begin + fusion_objects = I.(i, i, MTK._get_dual_cache(I)[2][i, i]) + for a in fusion_objects, b in fusion_objects, c in fusion_objects + @test SectorTestSuite.F_unitarity_test(a, b, c) end - F = hvcat(length(fs), Fblocks...) - @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12) end end + + i != j || continue # do this part only when off-diagonal + mod_objects = I.(i, j, MTK._get_dual_cache(I)[2][i, j]) + left_fusion_objects = I.(i, i, MTK._get_dual_cache(I)[2][i, i]) + right_fusion_objects = I.(j, j, MTK._get_dual_cache(I)[2][j, j]) + + # C x C x M -> M or D x D x Mop -> Mop + @testset "Unitarity of left module F-move ($i, $j)" begin + for a in left_fusion_objects, b in left_fusion_objects, A in mod_objects + @test SectorTestSuite.F_unitarity_test(a, b, A) + end + end + + # M x D x D -> M or Mop x C x C -> Mop + @testset "Unitarity of right module F-move ($i, $j)" begin + for A in mod_objects, b in right_fusion_objects, c in right_fusion_objects + @test SectorTestSuite.F_unitarity_test(A, b, c) + end + end + + # C x M x D -> M or D x Mop x C -> Mop + @testset "Unitarity of bimodule F-move ($i, $j)" begin + for a in left_fusion_objects, A in mod_objects, Ξ± in right_fusion_objects + @test SectorTestSuite.F_unitarity_test(a, A, Ξ±) + end + end + + @testset "Unitarity of mixed module F-move ($i, $j) and opposite ($j, $i)" begin + modop_objects = I.(j, i, MTK._get_dual_cache(I)[2][j, i]) + + # C x M x Mop -> C or D x Mop x M -> D + # M x Mop x C -> C or Mop x M x D -> D + # Mop x C x M -> D or M x D x Mop -> C + for a in left_fusion_objects, A in mod_objects, Aop in modop_objects + @test SectorTestSuite.F_unitarity_test(a, A, Aop) + @test SectorTestSuite.F_unitarity_test(A, Aop, a) + @test SectorTestSuite.F_unitarity_test(Aop, a, A) + end + + # M x Mop x M -> M or Mop x M x Mop -> Mop + for A in mod_objects, Aop in modop_objects + @test SectorTestSuite.F_unitarity_test(A, Aop, A) + end + end + end +end + +@testset "Triangle equation" begin + objects = collect(values(I)) + for a in objects, b in objects + a.j == b.i || continue # skip if not compatible + @test triangle_equation(a, b; atol = 1.0e-12, rtol = 1.0e-12) end +end - @testset "Pentagon equation" begin - for a in objects, b in objects, c in objects, d in objects - @test pentagon_equation(a, b, c, d; atol=1e-12, rtol=1e-12) +@testset "$Istr Pentagon equation" begin + objects = collect(values(I)) + for a in objects + for b in objects + a.j == b.i || continue # skip if not compatible + for c in objects + b.j == c.i || continue # skip if not compatible + for d in objects + c.j == d.i || continue # skip if not compatible + @test pentagon_equation(a, b, c, d; atol = 1.0e-12, rtol = 1.0e-12) + end + end end end end + +# include("TK_compat.jl") \ No newline at end of file diff --git a/test/test_aqua.jl b/test/test_aqua.jl index 14e5b7c..f57c8ba 100644 --- a/test/test_aqua.jl +++ b/test/test_aqua.jl @@ -3,5 +3,5 @@ using Aqua: Aqua using Test: @testset @testset "Code quality (Aqua.jl)" begin - Aqua.test_all(MultiTensorKit) + Aqua.test_all(MultiTensorKit) end