From b30f2055eb4661e87c97024fad03beb62f0f8d58 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 11 Aug 2026 08:26:39 +0200 Subject: [PATCH 01/36] specialise dim to storage type --- src/TensorKit.jl | 2 +- src/spaces/gradedspace.jl | 17 ++++++++++++++--- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 87a3a2380..929dd804e 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -44,7 +44,7 @@ export infimum, supremum, isisomorphic, ismonomorphic, isepimorphic export sectortype, sectors, hassector export unit, rightunit, leftunit, allunits, isunit, otimes, deligneproduct, timereversed export Nsymbol, Fsymbol, Rsymbol, Bsymbol, frobenius_schur_phase, frobenius_schur_indicator, twist, fusiontensor -export sectorscalartype, fusionscalartype, braidingscalartype +export sectorscalartype, fusionscalartype, braidingscalartype, dimscalartype # Export methods for fusion trees export fusiontrees, braid, permute, transpose diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 53f48dafd..d469648cc 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -89,9 +89,20 @@ GradedSpace(g::AbstractDict; dual::Bool = false) = GradedSpace(g...; dual = dual field(::Type{<:GradedSpace}) = ℂ InnerProductStyle(::Type{<:GradedSpace}) = EuclideanInnerProduct() -function dim(V::GradedSpace) - init = 0 * dim(first(allunits(sectortype(V)))) - return sum(c -> dim(c) * dim(V, c), sectors(V); init = init) +function dim(V::GradedSpace{I, <:AbstractDict}) where {I <: Sector} + init = zero(dimscalartype(I)) + return sum(((c, d),) -> dim(c) * d, V.dims; init) +end +function dim(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} + init = zero(dimscalartype(I)) + D = init + vals = values(I) + @inbounds for n in 1:N + d = V.dims[n] + iszero(d) && continue + D += dim(vals[n]) * d # dim(c) = dim(dual(c)) + end + return D end function dim(V::GradedSpace{I, <:AbstractDict}, c::I) where {I <: Sector} return get(V.dims, isdual(V) ? dual(c) : c, 0) From 41f43fdf997580b4128ee22cf9291570967f826e Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 11 Aug 2026 12:13:17 +0200 Subject: [PATCH 02/36] oplus and ominus --- src/spaces/gradedspace.jl | 72 +++++++++++++++++++++++++++++++-------- 1 file changed, 57 insertions(+), 15 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index d469648cc..6064eeaef 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -137,25 +137,67 @@ function unitspace(S::Type{<:GradedSpace{I}}) where {I <: Sector} end zerospace(S::Type{<:GradedSpace}) = S() -# TODO: the following methods can probably be implemented more efficiently for -# `FiniteGradedSpace`, but we don't expect them to be used often in hot loops, so -# these generic definitions (which are still quite efficient) are good for now. -function ⊕(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector} +function ⊕(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} + dual1 = isdual(V₁) + dual1 == isdual(V₂) || throw(SpaceMismatch("Direct sum of a vector space and a dual space does not exist")) + k1, k2 = V₁.dims.keys, V₂.dims.keys # already sorted + v1, v2 = V₁.dims.values, V₂.dims.values + n1, n2 = length(k1), length(k2) + ks, vs = Vector{I}(), Vector{Int}() + sizehint!(ks, n1 + n2) + sizehint!(vs, n1 + n2) + i, j = 1, 1 + @inbounds while i <= n1 && j <= n2 + if k1[i] == k2[j] + push!(ks, k1[i]); push!(vs, v1[i] + v2[j]); i += 1; j += 1 + elseif k1[i] < k2[j] + push!(ks, k1[i]); push!(vs, v1[i]); i += 1 + else + push!(ks, k2[j]); push!(vs, v2[j]); j += 1 + end + end + @inbounds while i <= n1 + push!(ks, k1[i]); push!(vs, v1[i]); i += 1 + end + @inbounds while j <= n2 + push!(ks, k2[j]); push!(vs, v2[j]); j += 1 + end + return typeof(V₁)(SectorDict{I, Int}(ks, vs), dual1) +end +function ⊕(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} dual1 = isdual(V₁) dual1 == isdual(V₂) || throw(SpaceMismatch("Direct sum of a vector space and a dual space does not exist")) - dims = SectorDict{I, Int}() - for c in union(sectors(V₁), sectors(V₂)) - cout = ifelse(dual1, dual(c), c) - dims[cout] = dim(V₁, c) + dim(V₂, c) + newdims = map(+, V₁.dims, V₂.dims) + return typeof(V₁)(newdims, dual1) +end +function ⊖(V::GradedSpace{I, <: Tuple}, W::GradedSpace{I, <: Tuple}) where {I <: Sector} + dualV = isdual(V) + V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)")) + newdims = map(-, V.dims, W.dims) + return typeof(V)(newdims, dualV) +end +function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) where {I <: Sector} + dualV = isdual(V) + V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)")) + kv, kw = V.dims.keys, W.dims.keys # already sorted + vv, vw = V.dims.values, W.dims.values + ks, vs = Vector{I}(), Vector{Int}() + nv, nw = length(kv), length(kw) + sizehint!(ks, nv) + sizehint!(vs, nv) + j = 1 + @inbounds for i in eachindex(kv) # keys(W) ⊆ keys(V) + d = vv[i] + if j <= nw && kw[j] == kv[i] + d -= vw[j] + j += 1 + end + if !iszero(d) + push!(ks, kv[i]); push!(vs, d) + end end - return typeof(V₁)(dims; dual = dual1) -end -function ⊖(V::GradedSpace{I}, W::GradedSpace{I}) where {I <: Sector} - dual = isdual(V) - V ≿ W && dual == isdual(W) || - throw(SpaceMismatch("$(W) is not a subspace of $(V)")) - return typeof(V)(c => dim(V, c) - dim(W, c) for c in sectors(V); dual) + return typeof(V)(SectorDict{I, Int}(ks, vs), dualV) end function fuse(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector} From 619804db43336769f4031ba76413944befa24d87 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 11 Aug 2026 12:21:28 +0200 Subject: [PATCH 03/36] infimum and supremum --- src/spaces/gradedspace.jl | 76 ++++++++++++++++++++++++++++++--------- 1 file changed, 60 insertions(+), 16 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 6064eeaef..46c94eaf7 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -171,7 +171,7 @@ function ⊕(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where newdims = map(+, V₁.dims, V₂.dims) return typeof(V₁)(newdims, dual1) end -function ⊖(V::GradedSpace{I, <: Tuple}, W::GradedSpace{I, <: Tuple}) where {I <: Sector} +function ⊖(V::GradedSpace{I, <:Tuple}, W::GradedSpace{I, <:Tuple}) where {I <: Sector} dualV = isdual(V) V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)")) newdims = map(-, V.dims, W.dims) @@ -210,23 +210,67 @@ function fuse(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector} return typeof(V₁)(dims) end -function infimum(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector} +function infimum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} Visdual = isdual(V₁) - Visdual == isdual(V₂) || - throw(SpaceMismatch("Infimum of space and dual space does not exist")) - return typeof(V₁)( - (Visdual ? dual(c) : c) => min(dim(V₁, c), dim(V₂, c)) - for c in intersect(sectors(V₁), sectors(V₂)); dual = Visdual - ) -end -function supremum(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector} + Visdual == isdual(V₂) || throw(SpaceMismatch("Infimum of space and dual space does not exist")) + newdims = map(min, V₁.dims, V₂.dims) + return typeof(V₁)(newdims, Visdual) +end +function infimum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} + Visdual = isdual(V₁) + Visdual == isdual(V₂) || throw(SpaceMismatch("Infimum of space and dual space does not exist")) + k1, k2 = V₁.dims.keys, V₂.dims.keys + v1, v2 = V₁.dims.values, V₂.dims.values + n1, n2 = length(k1), length(k2) + ks, vs = Vector{I}(), Vector{Int}() + i, j = 1, 1 + @inbounds while i <= n1 && j <= n2 + if k1[i] == k2[j] + m = min(v1[i], v2[j]) + if !iszero(m) + push!(ks, k1[i]); push!(vs, m) + end + i += 1; j += 1 + elseif k1[i] < k2[j] + i += 1 + else + j += 1 + end + end + return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual) +end +function supremum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} Visdual = isdual(V₁) - Visdual == isdual(V₂) || - throw(SpaceMismatch("Supremum of space and dual space does not exist")) - return typeof(V₁)( - (Visdual ? dual(c) : c) => max(dim(V₁, c), dim(V₂, c)) - for c in union(sectors(V₁), sectors(V₂)); dual = Visdual - ) + Visdual == isdual(V₂) || throw(SpaceMismatch("Supremum of space and dual space does not exist")) + newdims = map(max, V₁.dims, V₂.dims) + return typeof(V₁)(newdims, Visdual) +end +function supremum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} + Visdual = isdual(V₁) + Visdual == isdual(V₂) || throw(SpaceMismatch("Supremum of space and dual space does not exist")) + k1, k2 = V₁.dims.keys, V₂.dims.keys + v1, v2 = V₁.dims.values, V₂.dims.values + n1, n2 = length(k1), length(k2) + ks, vs = Vector{I}(), Vector{Int}() + sizehint!(ks, n1 + n2) + sizehint!(vs, n1 + n2) + i, j = 1, 1 + @inbounds while i <= n1 && j <= n2 + if k1[i] == k2[j] + push!(ks, k1[i]); push!(vs, max(v1[i], v2[j])); i += 1; j += 1 + elseif k1[i] < k2[j] + push!(ks, k1[i]); push!(vs, v1[i]); i += 1 + else + push!(ks, k2[j]); push!(vs, v2[j]); j += 1 + end + end + @inbounds while i <= n1 + push!(ks, k1[i]); push!(vs, v1[i]); i += 1 + end + @inbounds while j <= n2 + push!(ks, k2[j]); push!(vs, v2[j]); j += 1 + end + return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual) end hassector(V::GradedSpace{I}, s::I) where {I <: Sector} = dim(V, s) != 0 From 0967bd071455860e189687dfb347f267736d9e1e Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 11 Aug 2026 12:52:02 +0200 Subject: [PATCH 04/36] fuse --- src/spaces/gradedspace.jl | 47 ++++++++++++++++++++++++++++++++++----- 1 file changed, 41 insertions(+), 6 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 46c94eaf7..d32358405 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -200,14 +200,49 @@ function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) w return typeof(V)(SectorDict{I, Int}(ks, vs), dualV) end -function fuse(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector} - dims = SectorDict{I, Int}() - for a in sectors(V₁), b in sectors(V₂) - for c in a ⊗ b - dims[c] = get(dims, c, 0) + Nsymbol(a, b, c) * dim(V₁, a) * dim(V₂, b) +function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} + dual1, dual2 = isdual(V₁), isdual(V₂) + acc = Dict{I, Int}() # SectorDict `get` within the double for loop accumulates O(N^2) ` findindex` calls -> sort afterwards + k1, k2 = V₁.dims.keys, V₂.dims.keys + v1, v2 = V₁.dims.values, V₂.dims.values + @inbounds for n1 in eachindex(k1) + a0 = k1[n1]; d1 = v1[n1] + a = dual1 ? dual(a0) : a0 + for n2 in eachindex(k2) + b0 = k2[n2]; d2 = v2[n2] + b = dual2 ? dual(b0) : b0 + dab = d1 * d2 + for c in a ⊗ b + acc[c] = get(acc, c, 0) + Nsymbol(a, b, c) * dab + end + end + end + ks = sort!(collect(keys(acc))) + vs = [acc[k] for k in ks] + return typeof(V₁)(SectorDict{I, Int}(ks, vs), false) +end +function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} + vals = values(I) + dual1, dual2 = isdual(V₁), isdual(V₂) + newdims = zeros(Int, N) #TODO: is there a way to avoid dense storage even for sparse results? + @inbounds for n1 in 1:N + d1 = V₁.dims[n1] + iszero(d1) && continue + a0 = vals[n1] # avoid call to sectors(V₁) + a = dual1 ? dual(a0) : a0 + for n2 in 1:N + d2 = V₂.dims[n2] + iszero(d2) && continue + b0 = vals[n2] # idem for V₂ + b = dual2 ? dual(b0) : b0 + dab = d1 * d2 + for c in a ⊗ b + nc = findindex(vals, c) + newdims[nc] += Nsymbol(a, b, c) * dab + end end end - return typeof(V₁)(dims) + return typeof(V₁)(ntuple(i -> newdims[i], Val(N)), false) end function infimum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} From 60186d55213b442c3f46c1d0c3753e2c660e1050 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 11 Aug 2026 13:14:49 +0200 Subject: [PATCH 05/36] truncate_space --- src/factorizations/truncation.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 2facfd2e9..67dac8e56 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -37,6 +37,31 @@ _blocklength(ax::Base.OneTo, ind::AbstractVector{Bool}) = count(ind) function truncate_space(V::ElementarySpace, inds) return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in pairs(inds)) end +function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N} + vals = values(I) + dualV = isdual(V) + newdims = zeros(Int, N) + for (c, ind) in pairs(inds) + n_read = findindex(vals, dualV ? dual(c) : c) # dual-adjusted index for reading V.dims + n_write = findindex(vals, c) # output is never dual, so c is fine as-is + newdims[n_write] = _blocklength(V.dims[n_read], ind) # dim(c) = dim(dual(c)) + end + return typeof(V)(ntuple(i -> newdims[i], Val(N)), false) +end +function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector} + dualV = isdual(V) + ks, vs = Vector{I}(), Vector{Int}() # accumulate and sort once at the end + for (c, ind) in pairs(inds) + d = get(V.dims, dualV ? dual(c) : c, 0) + len = _blocklength(d, ind) + if !iszero(len) + push!(ks, c) + push!(vs, len) + end + end + perm = sortperm(ks) + return typeof(V)(SectorDict{I, Int}(ks[perm], vs[perm]), false) +end function truncate_domain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds) for (c, b) in blocks(tdst) From 8f90be2ddaf82a2d25de10c88abcc859d5d904d2 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 11 Aug 2026 13:20:36 +0200 Subject: [PATCH 06/36] restore binary search for sectordicts --- src/auxiliary/dicts.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl index 626599c01..d495103eb 100644 --- a/src/auxiliary/dicts.jl +++ b/src/auxiliary/dicts.jl @@ -89,14 +89,14 @@ end Base.empty(::SortedVectorDict, ::Type{K}, ::Type{V}) where {K, V} = SortedVectorDict{K, V}() Base.empty!(d::SortedVectorDict) = (empty!(d.keys); empty!(d.values); return d) -# _searchsortedfirst(v::Vector, k) = searchsortedfirst(v, k) -function _searchsortedfirst(v::Vector, k) - i = 1 - @inbounds while i <= length(v) && isless(v[i], k) - i += 1 - end - return i -end +_searchsortedfirst(v::Vector, k) = searchsortedfirst(v, k) +# function _searchsortedfirst(v::Vector, k) +# i = 1 +# @inbounds while i <= length(v) && isless(v[i], k) +# i += 1 +# end +# return i +# end function Base.delete!(d::SortedVectorDict{K}, k) where {K} key = convert(K, k) From ae10f5b565cc93c26caa33eab00d4c4ca8635e2f Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 11 Aug 2026 13:54:36 +0200 Subject: [PATCH 07/36] refactor sorted merge procedure --- src/auxiliary/dicts.jl | 49 ++++++++++++++++++++++ src/spaces/gradedspace.jl | 87 +++------------------------------------ 2 files changed, 55 insertions(+), 81 deletions(-) diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl index d495103eb..17297177e 100644 --- a/src/auxiliary/dicts.jl +++ b/src/auxiliary/dicts.jl @@ -186,6 +186,55 @@ function Base.:(==)(d1::SortedVectorDict, d2::SortedVectorDict) return true end +# merge over two SORTED vector pairs representing keys and values +# - combine(v1,v2): value for a key present in both operands +# - only1(v1) / only2(v2): value for a key present in only one operand; +# pass `nothing` to drop such keys entirely (e.g. for an intersection) +# zero results are dropped (either from `combine` or `only1`/`only2`), matching how GradedSpace never stores an explicit zero dimension +# k1 and k2 originate from GradedSpace.dims.keys, which are guaranteed to be sorted +function _sortedmerge(k1::Vector{I}, v1::Vector{Int}, k2::Vector{I}, v2::Vector{Int}, combine, only1, only2) where {I} + n1, n2 = length(k1), length(k2) + ks, vs = Vector{I}(), Vector{Int}() + sizehint!(ks, n1 + n2) + sizehint!(vs, n1 + n2) + i, j = 1, 1 + @inbounds while i <= n1 && j <= n2 + if k1[i] == k2[j] + d = combine(v1[i], v2[j]) + if !iszero(d) + push!(ks, k1[i]) + push!(vs, d) + end + i += 1 + j += 1 + elseif k1[i] < k2[j] + _mergeonly!(ks, vs, k1[i], v1[i], only1) + i += 1 + else + _mergeonly!(ks, vs, k2[j], v2[j], only2) + j += 1 + end + end + @inbounds while i <= n1 + _mergeonly!(ks, vs, k1[i], v1[i], only1) + i += 1 + end + @inbounds while j <= n2 + _mergeonly!(ks, vs, k2[j], v2[j], only2) + j += 1 + end + return ks, vs +end +@inline _mergeonly!(ks, vs, k, v, ::Nothing) = nothing +@inline function _mergeonly!(ks, vs, k, v, f) + d = f(v) + if !iszero(d) + push!(ks, k) + push!(vs, d) + end + return nothing +end + """ Hashed(value, hashfunction = Base.hash, isequal = Base.isequal) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index d32358405..5ca39f912 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -140,28 +140,7 @@ zerospace(S::Type{<:GradedSpace}) = S() function ⊕(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} dual1 = isdual(V₁) dual1 == isdual(V₂) || throw(SpaceMismatch("Direct sum of a vector space and a dual space does not exist")) - k1, k2 = V₁.dims.keys, V₂.dims.keys # already sorted - v1, v2 = V₁.dims.values, V₂.dims.values - n1, n2 = length(k1), length(k2) - ks, vs = Vector{I}(), Vector{Int}() - sizehint!(ks, n1 + n2) - sizehint!(vs, n1 + n2) - i, j = 1, 1 - @inbounds while i <= n1 && j <= n2 - if k1[i] == k2[j] - push!(ks, k1[i]); push!(vs, v1[i] + v2[j]); i += 1; j += 1 - elseif k1[i] < k2[j] - push!(ks, k1[i]); push!(vs, v1[i]); i += 1 - else - push!(ks, k2[j]); push!(vs, v2[j]); j += 1 - end - end - @inbounds while i <= n1 - push!(ks, k1[i]); push!(vs, v1[i]); i += 1 - end - @inbounds while j <= n2 - push!(ks, k2[j]); push!(vs, v2[j]); j += 1 - end + ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, +, identity, identity) return typeof(V₁)(SectorDict{I, Int}(ks, vs), dual1) end function ⊕(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} @@ -180,23 +159,7 @@ end function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) where {I <: Sector} dualV = isdual(V) V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)")) - kv, kw = V.dims.keys, W.dims.keys # already sorted - vv, vw = V.dims.values, W.dims.values - ks, vs = Vector{I}(), Vector{Int}() - nv, nw = length(kv), length(kw) - sizehint!(ks, nv) - sizehint!(vs, nv) - j = 1 - @inbounds for i in eachindex(kv) # keys(W) ⊆ keys(V) - d = vv[i] - if j <= nw && kw[j] == kv[i] - d -= vw[j] - j += 1 - end - if !iszero(d) - push!(ks, kv[i]); push!(vs, d) - end - end + ks, vs = _sortedmerge(V.dims.keys, V.dims.values, W.dims.keys, W.dims.values, -, identity, nothing) return typeof(V)(SectorDict{I, Int}(ks, vs), dualV) end @@ -217,14 +180,14 @@ function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorD end end end - ks = sort!(collect(keys(acc))) + ks = sort!(collect(keys(acc))) #TODO: sortperm? vs = [acc[k] for k in ks] return typeof(V₁)(SectorDict{I, Int}(ks, vs), false) end function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} vals = values(I) dual1, dual2 = isdual(V₁), isdual(V₂) - newdims = zeros(Int, N) #TODO: is there a way to avoid dense storage even for sparse results? + newdims = zeros(Int, N) @inbounds for n1 in 1:N d1 = V₁.dims[n1] iszero(d1) && continue @@ -254,24 +217,7 @@ end function infimum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} Visdual = isdual(V₁) Visdual == isdual(V₂) || throw(SpaceMismatch("Infimum of space and dual space does not exist")) - k1, k2 = V₁.dims.keys, V₂.dims.keys - v1, v2 = V₁.dims.values, V₂.dims.values - n1, n2 = length(k1), length(k2) - ks, vs = Vector{I}(), Vector{Int}() - i, j = 1, 1 - @inbounds while i <= n1 && j <= n2 - if k1[i] == k2[j] - m = min(v1[i], v2[j]) - if !iszero(m) - push!(ks, k1[i]); push!(vs, m) - end - i += 1; j += 1 - elseif k1[i] < k2[j] - i += 1 - else - j += 1 - end - end + ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, min, nothing, nothing) return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual) end function supremum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} @@ -283,28 +229,7 @@ end function supremum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} Visdual = isdual(V₁) Visdual == isdual(V₂) || throw(SpaceMismatch("Supremum of space and dual space does not exist")) - k1, k2 = V₁.dims.keys, V₂.dims.keys - v1, v2 = V₁.dims.values, V₂.dims.values - n1, n2 = length(k1), length(k2) - ks, vs = Vector{I}(), Vector{Int}() - sizehint!(ks, n1 + n2) - sizehint!(vs, n1 + n2) - i, j = 1, 1 - @inbounds while i <= n1 && j <= n2 - if k1[i] == k2[j] - push!(ks, k1[i]); push!(vs, max(v1[i], v2[j])); i += 1; j += 1 - elseif k1[i] < k2[j] - push!(ks, k1[i]); push!(vs, v1[i]); i += 1 - else - push!(ks, k2[j]); push!(vs, v2[j]); j += 1 - end - end - @inbounds while i <= n1 - push!(ks, k1[i]); push!(vs, v1[i]); i += 1 - end - @inbounds while j <= n2 - push!(ks, k2[j]); push!(vs, v2[j]); j += 1 - end + ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, max, identity, identity) return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual) end From 47a5f3437054d7533c454276ca95f604899e6b0a Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 11 Aug 2026 15:21:40 +0200 Subject: [PATCH 08/36] speed up fuse slightly with sortperm --- src/spaces/gradedspace.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 5ca39f912..0a79293c2 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -180,9 +180,10 @@ function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorD end end end - ks = sort!(collect(keys(acc))) #TODO: sortperm? - vs = [acc[k] for k in ks] - return typeof(V₁)(SectorDict{I, Int}(ks, vs), false) + ks0 = collect(keys(acc)) + vs0 = collect(values(acc)) + perm = sortperm(ks0) + return typeof(V₁)(SectorDict{I, Int}(ks0[perm], vs0[perm]), false) end function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} vals = values(I) From 411dc4472b20d214bb851e53b1468a5938a7c29d Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 11 Aug 2026 15:21:55 +0200 Subject: [PATCH 09/36] import thing --- src/factorizations/factorizations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index fbe87a63a..80d96ff9a 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -6,7 +6,7 @@ module Factorizations export copy_oftype, factorisation_scalartype, one!, truncspace using ..TensorKit -using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, +using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, findindex, blocktype, foreachblock, one!, similar_diagonal, similarstoragetype From 0adec3ba4357537db2857dac9c25125673302a0d Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 12 Aug 2026 11:47:24 +0200 Subject: [PATCH 10/36] splat with type annotation above Val --- src/factorizations/truncation.jl | 2 +- src/spaces/gradedspace.jl | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 67dac8e56..f3c64ed7c 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -46,7 +46,7 @@ function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sec n_write = findindex(vals, c) # output is never dual, so c is fine as-is newdims[n_write] = _blocklength(V.dims[n_read], ind) # dim(c) = dim(dual(c)) end - return typeof(V)(ntuple(i -> newdims[i], Val(N)), false) + return typeof(V)((newdims...,)::NTuple{N, Int}, false) end function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector} dualV = isdual(V) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 0a79293c2..3e363d19c 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -30,17 +30,17 @@ end sectortype(::Type{<:GradedSpace{I}}) where {I <: Sector} = I function GradedSpace{I, NTuple{N, Int}}(dims; dual::Bool = false) where {I, N} - d = ntuple(n -> 0, N) - isset = ntuple(n -> false, N) + d = zeros(Int, N) + isset = falses(N) for (c, dc) in dims k = convert(I, c) i = findindex(values(I), k) - k = dc < 0 && throw(ArgumentError(lazy"Sector $k has negative dimension $dc")) + dc < 0 && throw(ArgumentError(lazy"Sector $k has negative dimension $dc")) isset[i] && throw(ArgumentError(lazy"Sector $c appears multiple times")) - isset = TupleTools.setindex(isset, true, i) - d = TupleTools.setindex(d, dc, i) + isset[i] = true + d[i] = dc end - return GradedSpace{I, NTuple{N, Int}}(d, dual) + return GradedSpace{I, NTuple{N, Int}}((d...,)::NTuple{N, Int}, dual) end function GradedSpace{I, NTuple{N, Int}}(dims::Pair; dual::Bool = false) where {I, N} return GradedSpace{I, NTuple{N, Int}}((dims,); dual = dual) @@ -206,7 +206,7 @@ function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{ end end end - return typeof(V₁)(ntuple(i -> newdims[i], Val(N)), false) + return typeof(V₁)((newdims...,)::NTuple{N, Int}, false) end function infimum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} From 13677f3129b19757c7ce4fd72a0129b577c64e88 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 21 Aug 2026 15:33:03 +0200 Subject: [PATCH 11/36] actually don't splat, but construct directly where previously a vector was made --- src/factorizations/truncation.jl | 2 +- src/spaces/gradedspace.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index f3c64ed7c..3ec0457a9 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -46,7 +46,7 @@ function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sec n_write = findindex(vals, c) # output is never dual, so c is fine as-is newdims[n_write] = _blocklength(V.dims[n_read], ind) # dim(c) = dim(dual(c)) end - return typeof(V)((newdims...,)::NTuple{N, Int}, false) + return typeof(V)(NTuple{N, Int}(newdims), false) end function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector} dualV = isdual(V) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 3e363d19c..ede620e62 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -40,7 +40,7 @@ function GradedSpace{I, NTuple{N, Int}}(dims; dual::Bool = false) where {I, N} isset[i] = true d[i] = dc end - return GradedSpace{I, NTuple{N, Int}}((d...,)::NTuple{N, Int}, dual) + return GradedSpace{I, NTuple{N, Int}}(NTuple{N, Int}(d), dual) end function GradedSpace{I, NTuple{N, Int}}(dims::Pair; dual::Bool = false) where {I, N} return GradedSpace{I, NTuple{N, Int}}((dims,); dual = dual) @@ -206,7 +206,7 @@ function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{ end end end - return typeof(V₁)((newdims...,)::NTuple{N, Int}, false) + return typeof(V₁)(NTuple{N, Int}(newdims), false) end function infimum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} From 399e121f677afbd330df396feaed54b97918c2c4 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 21 Aug 2026 16:17:06 +0200 Subject: [PATCH 12/36] make slightly more readable maybe perhaps --- src/spaces/gradedspace.jl | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index ede620e62..bdea064e1 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -168,13 +168,13 @@ function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorD acc = Dict{I, Int}() # SectorDict `get` within the double for loop accumulates O(N^2) ` findindex` calls -> sort afterwards k1, k2 = V₁.dims.keys, V₂.dims.keys v1, v2 = V₁.dims.values, V₂.dims.values - @inbounds for n1 in eachindex(k1) - a0 = k1[n1]; d1 = v1[n1] - a = dual1 ? dual(a0) : a0 - for n2 in eachindex(k2) - b0 = k2[n2]; d2 = v2[n2] - b = dual2 ? dual(b0) : b0 - dab = d1 * d2 + @inbounds for na in eachindex(k1) + a₀, da = k1[na], v1[na] + a = dual1 ? dual(a₀) : a₀ + for nb in eachindex(k2) + b₀, db = k2[nb], v2[nb] + b = dual2 ? dual(b₀) : b₀ + dab = da * db for c in a ⊗ b acc[c] = get(acc, c, 0) + Nsymbol(a, b, c) * dab end @@ -189,17 +189,17 @@ function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{ vals = values(I) dual1, dual2 = isdual(V₁), isdual(V₂) newdims = zeros(Int, N) - @inbounds for n1 in 1:N - d1 = V₁.dims[n1] - iszero(d1) && continue - a0 = vals[n1] # avoid call to sectors(V₁) - a = dual1 ? dual(a0) : a0 - for n2 in 1:N - d2 = V₂.dims[n2] - iszero(d2) && continue - b0 = vals[n2] # idem for V₂ - b = dual2 ? dual(b0) : b0 - dab = d1 * d2 + @inbounds for na in 1:N + da = V₁.dims[na] + iszero(da) && continue + a₀ = vals[na] # avoid call to sectors(V₁) + a = dual1 ? dual(a₀) : a₀ + for nb in 1:N + db = V₂.dims[nb] + iszero(db) && continue + b₀ = vals[nb] # idem for V₂ + b = dual2 ? dual(b₀) : b₀ + dab = da * db for c in a ⊗ b nc = findindex(vals, c) newdims[nc] += Nsymbol(a, b, c) * dab From d2a41f16fff0ff10108dda3966433a68b1269ef0 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 25 Aug 2026 17:46:59 +0200 Subject: [PATCH 13/36] truncate_space always has non-dual entry spaces --- src/factorizations/truncation.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 3ec0457a9..87971ca88 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -39,20 +39,18 @@ function truncate_space(V::ElementarySpace, inds) end function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N} vals = values(I) - dualV = isdual(V) newdims = zeros(Int, N) for (c, ind) in pairs(inds) - n_read = findindex(vals, dualV ? dual(c) : c) # dual-adjusted index for reading V.dims - n_write = findindex(vals, c) # output is never dual, so c is fine as-is - newdims[n_write] = _blocklength(V.dims[n_read], ind) # dim(c) = dim(dual(c)) + d = dim(V, c) + n_write = findindex(vals, c) + newdims[n_write] = _blocklength(d, ind) end return typeof(V)(NTuple{N, Int}(newdims), false) end function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector} - dualV = isdual(V) ks, vs = Vector{I}(), Vector{Int}() # accumulate and sort once at the end for (c, ind) in pairs(inds) - d = get(V.dims, dualV ? dual(c) : c, 0) + d = dim(V, c) len = _blocklength(d, ind) if !iszero(len) push!(ks, c) From d7bdb07189022986ebff1be5a6846ff1fad20775 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 25 Aug 2026 18:04:44 +0200 Subject: [PATCH 14/36] overkill iszero check in dim --- src/spaces/gradedspace.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index bdea064e1..cde32d709 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -97,10 +97,8 @@ function dim(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} init = zero(dimscalartype(I)) D = init vals = values(I) - @inbounds for n in 1:N - d = V.dims[n] - iszero(d) && continue - D += dim(vals[n]) * d # dim(c) = dim(dual(c)) + @inbounds for (c, d) in zip(values(I), V.dims) + D += dim(c) * d end return D end From 271e96b4bfe1004661712c30652b934e7288c357 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 26 Aug 2026 09:48:22 +0200 Subject: [PATCH 15/36] assert truncate_space spaces being non-dual --- src/factorizations/truncation.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 87971ca88..03e636298 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -35,9 +35,11 @@ _blocklength(ax::Base.OneTo, ind::AbstractVector{<:Integer}) = length(ind) _blocklength(ax::Base.OneTo, ind::AbstractVector{Bool}) = count(ind) function truncate_space(V::ElementarySpace, inds) + @assert !isdual(V) return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in pairs(inds)) end function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N} + @assert !isdual(V) vals = values(I) newdims = zeros(Int, N) for (c, ind) in pairs(inds) @@ -48,6 +50,7 @@ function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sec return typeof(V)(NTuple{N, Int}(newdims), false) end function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector} + @assert !isdual(V) ks, vs = Vector{I}(), Vector{Int}() # accumulate and sort once at the end for (c, ind) in pairs(inds) d = dim(V, c) From 335785f282cf6ee21a856c04e5c3f479d06c13a3 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 28 Aug 2026 17:18:18 +0200 Subject: [PATCH 16/36] introduce `sectorstoragetype` --- src/spaces/gradedspace.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index cde32d709..b7b35c8f9 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -323,12 +323,20 @@ specify `D`. const Vect = SpaceTable() Base.getindex(::SpaceTable) = ComplexSpace Base.getindex(::SpaceTable, ::Type{Trivial}) = ComplexSpace -function Base.getindex(::SpaceTable, I::Type{<:Sector}) +Base.getindex(::SpaceTable, I::Type{<:Sector}) = GradedSpace{I, sectorstoragetype(I)} + +""" + sectorstoragetype(I::Type{<:Sector}) -> Type + +The storage type `D` used for the `dims` field of `GradedSpace{I, D}`. +This is `NTuple{N,Int}` with `N = length(values(I))` if `I` has a finite, known length, +or `SectorDict{I,Int}` otherwise. +""" +Base.@assume_effects :foldable function sectorstoragetype(::Type{I}) where {I <: Sector} if Base.IteratorSize(values(I)) isa Union{HasLength, HasShape} - N = length(values(I)) - return GradedSpace{I, NTuple{N, Int}} + return NTuple{length(values(I)), Int} else - return GradedSpace{I, SectorDict{I, Int}} + return SectorDict{I, Int} end end From d79931ec59734c0dbefb589d83010b295b586677 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 28 Aug 2026 18:12:59 +0200 Subject: [PATCH 17/36] specialise truncation code to sectorstoragetype --- src/factorizations/factorizations.jl | 2 +- src/factorizations/pullbacks.jl | 12 ++- src/factorizations/truncation.jl | 152 ++++++++++++++++++++------- 3 files changed, 123 insertions(+), 43 deletions(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 80d96ff9a..b16a28e9b 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -8,7 +8,7 @@ export copy_oftype, factorisation_scalartype, one!, truncspace using ..TensorKit using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, findindex, blocktype, foreachblock, one!, - similar_diagonal, similarstoragetype + similar_diagonal, similarstoragetype, sectorstoragetype using LinearAlgebra: LinearAlgebra, BlasFloat, Diagonal, svdvals, svdvals!, eigen, eigen!, diff --git a/src/factorizations/pullbacks.jl b/src/factorizations/pullbacks.jl index b910a184c..68ac57bb4 100644 --- a/src/factorizations/pullbacks.jl +++ b/src/factorizations/pullbacks.jl @@ -24,16 +24,22 @@ for pullback! in (:qr_null_pullback!, :lq_null_pullback!) return Δt end end -_notrunc_ind(t) = SectorDict(c => Colon() for c in blocksectors(t)) +function _notrunc_ind(t) + I = sectortype(t) + return _builddensemap(sectorstoragetype(I), I, blocks(t), Colon) do _, _ + Colon() + end +end for pullback! in (:svd_pullback!, :eig_pullback!, :eigh_pullback!) @eval function MAK.$pullback!( Δt::AbstractTensorMap, t::AbstractTensorMap, F, ΔF, inds = _notrunc_ind(t); kwargs... ) + Isec = sectortype(t) foreachblock(Δt, t) do c, (Δb, b) - haskey(inds, c) || return nothing - ind = inds[c] + ind = _denseget(inds, Isec, c) + isnothing(ind) && return nothing Fc = block.(F, Ref(c)) ΔFc = block.(ΔF, Ref(c)) MAK.$pullback!(Δb, b, Fc, ΔFc, ind; kwargs...) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 03e636298..d2e455ebd 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -34,15 +34,55 @@ _blocklength(ax, ind) = length(ax[ind]) _blocklength(ax::Base.OneTo, ind::AbstractVector{<:Integer}) = length(ind) _blocklength(ax::Base.OneTo, ind::AbstractVector{Bool}) = count(ind) +# TODO: it quacks like a duck, just define a subtype of AbstractDict? +# represent the sector-index mapping as Vector{Union{Nothing, V}} where V is the type of the index +# mapping is indexed through findindex +# the type V is needed because the concrete type of ind depends on the strategy (except for intersect/union) +_densenew(::Type{I}, ::Type{V}) where {I <: Sector, V} = + Vector{Union{Nothing, V}}(nothing, length(values(I))) + +function _denseset!(v::Vector, ::Type{I}, c::I, val) where {I <: Sector} + v[findindex(values(I), c)] = val + return v +end +_denseget(v::Vector, ::Type{I}, c::I) where {I <: Sector} = v[findindex(values(I), c)] +function _densepairs(v::Vector, ::Type{I}) where {I <: Sector} + vals = values(I) + return (vals[i] => x for (i, x) in enumerate(v) if !isnothing(x)) +end +_densekeys(v::Vector, ::Type{I}) where {I <: Sector} = (c for (c, _) in _densepairs(v, I)) + +# fallbacks to catch SectorVector/SectorDict, even for NTuple sectorstoragetype +_denseget(v, ::Type{I}, c::I) where {I <: Sector} = get(v, c, nothing) +_densekeys(v, ::Type{I}) where {I <: Sector} = keys(v) +_densepairs(v, ::Type{I}) where {I <: Sector} = pairs(v) + +# builds either a dense Vector or SectorDict based on sectorstoragetype +# mapping each (c, v) pair's sector c to f(c, v) +# so every `findtruncated` method shares one output-construction path +# `pairsiter` are c => v pairs, can be c => nothing for NoTruncation/TruncationIntersection/TruncationUnion +function _builddensemap(f, ::Type{D}, ::Type{I}, pairsiter, ::Type{V}) where {D <: Tuple, I <: Sector, V} + d = _densenew(I, V) + for (c, v) in pairsiter + _denseset!(d, I, c, f(c, v)) + end + return d +end +function _builddensemap(f, ::Type{D}, ::Type{I}, pairsiter, ::Type{V}) where {D <: SectorDict, I <: Sector, V} + return SectorDict(c => f(c, v) for (c, v) in pairsiter) # V unused +end + function truncate_space(V::ElementarySpace, inds) @assert !isdual(V) - return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in pairs(inds)) + I = sectortype(V) + @assert I == Trivial + return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in _densepairs(inds, I)) end function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N} @assert !isdual(V) vals = values(I) newdims = zeros(Int, N) - for (c, ind) in pairs(inds) + for (c, ind) in _densepairs(inds, I) d = dim(V, c) n_write = findindex(vals, c) newdims[n_write] = _blocklength(d, ind) @@ -65,27 +105,30 @@ function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Secto end function truncate_domain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds) + Isec = sectortype(tdst) for (c, b) in blocks(tdst) - I = get(inds, c, nothing) - @assert !isnothing(I) + I = _denseget(inds, Isec, c) + @assert !isnothing(I) # kept for safety, but should be guaranteed by _densepairs b′ = block(tsrc, c) b .= view(b′, :, I) end return tdst end function truncate_codomain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds) + Isec = sectortype(tdst) for (c, b) in blocks(tdst) - I = get(inds, c, nothing) - @assert !isnothing(I) + I = _denseget(inds, Isec, c) + @assert !isnothing(I) # kept for safety, but should be guaranteed by _densepairs b′ = block(tsrc, c) b .= view(b′, I, :) end return tdst end function truncate_diagonal!(Ddst::DiagonalTensorMap, Dsrc::DiagonalTensorMap, inds) + Isec = sectortype(Ddst) for (c, b) in blocks(Ddst) - I = get(inds, c, nothing) - @assert !isnothing(I) + I = _denseget(inds, Isec, c) + @assert !isnothing(I) # kept for safety, but should be guaranteed by _densepairs diagview(b) .= view(diagview(block(Dsrc, c)), I) end return Ddst @@ -140,7 +183,10 @@ end function MAK.truncate( ::typeof(left_null!), (U, S)::NTuple{2, AbstractTensorMap}, strategy::NoTruncation ) - ind = SectorDict(c => (size(b, 2) + 1):size(b, 1) for (c, b) in blocks(S)) + I = sectortype(S) + ind = _builddensemap(sectorstoragetype(I), I, blocks(S), UnitRange{Int}) do _, b + (size(b, 2) + 1):size(b, 1) + end V_truncated = truncate_space(space(S, 1), ind) Ũ = similar(U, codomain(U) ← V_truncated) truncate_domain!(Ũ, U, ind) @@ -149,7 +195,10 @@ end function MAK.truncate( ::typeof(right_null!), (S, Vᴴ)::NTuple{2, AbstractTensorMap}, strategy::NoTruncation ) - ind = SectorDict(c => (size(b, 1) + 1):size(b, 2) for (c, b) in blocks(S)) + I = sectortype(S) + ind = _builddensemap(sectorstoragetype(I), I, blocks(S), UnitRange{Int}) do _, b + (size(b, 1) + 1):size(b, 2) + end V_truncated = truncate_space(dual(space(S, 2)), ind) Ṽᴴ = similar(Vᴴ, V_truncated ← domain(Vᴴ)) truncate_codomain!(Ṽᴴ, Vᴴ, ind) @@ -186,7 +235,10 @@ function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationStrateg end function MAK.findtruncated(values::SectorVector, ::NoTruncation) - return SectorDict(c => Colon() for c in keys(values)) + I = sectortype(values) + return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in keys(values)), Colon) do _, _ + Colon() + end end # Need to select the first k values here after sorting across blocks, weighted by quantum dimension @@ -232,18 +284,29 @@ MAK.findtruncated_svd(values::SectorVector, strategy::TruncationByOrder) = MAK.findtruncated(values, strategy) function MAK.findtruncated(values::SectorVector, strategy::TruncationByFilter) - return SectorDict(c => findall(strategy.filter, d) for (c, d) in pairs(values)) + I = sectortype(values) + return _builddensemap(sectorstoragetype(I), I, pairs(values), Vector{Int}) do _, v + findall(strategy.filter, v) + end end function MAK.findtruncated(values::SectorVector, strategy::TruncationByValue) + I = sectortype(values) atol = rtol_to_atol(values, strategy.p, strategy.atol, strategy.rtol) strategy′ = trunctol(; atol, strategy.by, strategy.keep_below) - return SectorDict(c => MAK.findtruncated(d, strategy′) for (c, d) in pairs(values)) + V = Base.promote_op(MAK.findtruncated, valtype(values), typeof(strategy′)) + return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do _, v + MAK.findtruncated(v, strategy′) + end end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationByValue) + I = sectortype(values) atol = rtol_to_atol(values, strategy.p, strategy.atol, strategy.rtol) strategy′ = trunctol(; atol, strategy.by, strategy.keep_below) - return SectorDict(c => MAK.findtruncated_svd(d, strategy′) for (c, d) in pairs(values)) + V = Base.promote_op(MAK.findtruncated_svd, valtype(values), typeof(strategy′)) + return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do _, v + MAK.findtruncated_svd(v, strategy′) + end end # Need to select the first k values here after sorting by error across blocks, @@ -285,14 +348,24 @@ MAK.findtruncated_svd(values::SectorVector, strategy::TruncationByError) = MAK.findtruncated(values, strategy) function MAK.findtruncated(values::SectorVector, strategy::TruncationSpace) - sectortype(values) == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values")) + I = sectortype(values) + I == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values")) blockstrategy(c) = truncrank(dim(strategy.space, c); strategy.by, strategy.rev) - return SectorDict(c => MAK.findtruncated(d, blockstrategy(c)) for (c, d) in pairs(values)) + Vstrategy = Base.promote_op(blockstrategy, I) + V = Base.promote_op(MAK.findtruncated, valtype(values), Vstrategy) + return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do c, v + MAK.findtruncated(v, blockstrategy(c)) + end end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationSpace) - sectortype(values) == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values")) + I = sectortype(values) + I == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values")) blockstrategy(c) = truncrank(dim(strategy.space, c); strategy.by, strategy.rev) - return SectorDict(c => MAK.findtruncated_svd(d, blockstrategy(c)) for (c, d) in pairs(values)) + Vstrategy = Base.promote_op(blockstrategy, I) + V = Base.promote_op(MAK.findtruncated_svd, valtype(values), Vstrategy) + return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do c, v + MAK.findtruncated_svd(v, blockstrategy(c)) + end end # The implementations below assume that the `SectorDict` always contains an entry for every block sector @@ -300,40 +373,40 @@ end # This is always the case in the implementations above. function MAK.findtruncated(values::SectorVector, strategy::TruncationIntersection) + I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated, values), strategy.components) - @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now" - sectors = keys(first(inds)) - vals = map(keys(first(inds))) do c - mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_intersect, inds) + @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now" + sectors = collect(_densekeys(first(inds), I)) + return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ + mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_intersect, inds) end - return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals) end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationIntersection) + I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated_svd, values), strategy.components) - @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now" - sectors = keys(first(inds)) - vals = map(keys(first(inds))) do c - mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_intersect, inds) + @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now" + sectors = collect(_densekeys(first(inds), I)) + return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ + mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_intersect, inds) end - return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals) end function MAK.findtruncated(values::SectorVector, strategy::TruncationUnion) + I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated, values), strategy.components) - @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now" - sectors = keys(first(inds)) - vals = map(keys(first(inds))) do c - mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_union, inds) + @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now" + sectors = collect(_densekeys(first(inds), I)) + return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ + mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_union, inds) end - return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals) end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationUnion) + I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated_svd, values), strategy.components) - @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now" - sectors = keys(first(inds)) - vals = map(keys(first(inds))) do c - mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_union, inds) + @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now" + sectors = collect(_densekeys(first(inds), I)) + return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ + mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_union, inds) end - return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals) end # Truncation error @@ -341,7 +414,8 @@ end MAK.truncation_error(values::SectorVector, ind) = MAK.truncation_error!(copy(values), ind) function MAK.truncation_error!(values::SectorVector, ind) - for (c, ind_c) in pairs(ind) + Isec = sectortype(values) + for (c, ind_c) in _densepairs(ind, Isec) v = values[c] v[ind_c] .= zero(eltype(v)) end From d42aa28bea2e3a3b15c339112e21011d92d3c768 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 28 Aug 2026 18:31:07 +0200 Subject: [PATCH 18/36] put `_sortedmerge` in `Base.mergewith` and use where possible --- src/auxiliary/dicts.jl | 67 ++++++++++++++++++++++----------------- src/spaces/gradedspace.jl | 12 +++---- 2 files changed, 42 insertions(+), 37 deletions(-) diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl index 17297177e..c84a9e1e5 100644 --- a/src/auxiliary/dicts.jl +++ b/src/auxiliary/dicts.jl @@ -186,17 +186,19 @@ function Base.:(==)(d1::SortedVectorDict, d2::SortedVectorDict) return true end -# merge over two SORTED vector pairs representing keys and values -# - combine(v1,v2): value for a key present in both operands -# - only1(v1) / only2(v2): value for a key present in only one operand; -# pass `nothing` to drop such keys entirely (e.g. for an intersection) -# zero results are dropped (either from `combine` or `only1`/`only2`), matching how GradedSpace never stores an explicit zero dimension -# k1 and k2 originate from GradedSpace.dims.keys, which are guaranteed to be sorted -function _sortedmerge(k1::Vector{I}, v1::Vector{Int}, k2::Vector{I}, v2::Vector{Int}, combine, only1, only2) where {I} +# merge over two SectorDicts +# the intersect case for infimum is kind of tricky, so there's an extra bool +# to indicate keeping keys that are only present in one of the two dicts +# zero results are dropped, matching how GradedSpace never stores an explicit zero dimension +function _sortedmerge( + combine, ::Val{keepunique}, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V} + ) where {keepunique, K, V} + k1, v1 = d1.keys, d1.values + k2, v2 = d2.keys, d2.values n1, n2 = length(k1), length(k2) - ks, vs = Vector{I}(), Vector{Int}() - sizehint!(ks, n1 + n2) - sizehint!(vs, n1 + n2) + ks, vs = Vector{K}(), Vector{V}() + sizehint!(ks, keepunique ? n1 + n2 : min(n1, n2)) + sizehint!(vs, keepunique ? n1 + n2 : min(n1, n2)) i, j = 1, 1 @inbounds while i <= n1 && j <= n2 if k1[i] == k2[j] @@ -208,33 +210,40 @@ function _sortedmerge(k1::Vector{I}, v1::Vector{Int}, k2::Vector{I}, v2::Vector{ i += 1 j += 1 elseif k1[i] < k2[j] - _mergeonly!(ks, vs, k1[i], v1[i], only1) + if keepunique + push!(ks, k1[i]) + push!(vs, v1[i]) + end i += 1 else - _mergeonly!(ks, vs, k2[j], v2[j], only2) + if keepunique + push!(ks, k2[j]) + push!(vs, v2[j]) + end j += 1 end end - @inbounds while i <= n1 - _mergeonly!(ks, vs, k1[i], v1[i], only1) - i += 1 - end - @inbounds while j <= n2 - _mergeonly!(ks, vs, k2[j], v2[j], only2) - j += 1 - end - return ks, vs -end -@inline _mergeonly!(ks, vs, k, v, ::Nothing) = nothing -@inline function _mergeonly!(ks, vs, k, v, f) - d = f(v) - if !iszero(d) - push!(ks, k) - push!(vs, d) + if keepunique + @inbounds while i <= n1 + push!(ks, k1[i]) + push!(vs, v1[i]) + i += 1 + end + @inbounds while j <= n2 + push!(ks, k2[j]) + push!(vs, v2[j]) + j += 1 + end end - return nothing + return SortedVectorDict{K, V}(ks, vs) end +Base.mergewith(combine, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {K, V} = + _sortedmerge(combine, Val(true), d1, d2) + +_sortedintersect(combine, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {K, V} = + _sortedmerge(combine, Val(false), d1, d2) + """ Hashed(value, hashfunction = Base.hash, isequal = Base.isequal) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index b7b35c8f9..25bce0ec5 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -138,8 +138,7 @@ zerospace(S::Type{<:GradedSpace}) = S() function ⊕(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} dual1 = isdual(V₁) dual1 == isdual(V₂) || throw(SpaceMismatch("Direct sum of a vector space and a dual space does not exist")) - ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, +, identity, identity) - return typeof(V₁)(SectorDict{I, Int}(ks, vs), dual1) + return typeof(V₁)(mergewith(+, V₁.dims, V₂.dims), dual1) end function ⊕(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} dual1 = isdual(V₁) @@ -157,8 +156,7 @@ end function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) where {I <: Sector} dualV = isdual(V) V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)")) - ks, vs = _sortedmerge(V.dims.keys, V.dims.values, W.dims.keys, W.dims.values, -, identity, nothing) - return typeof(V)(SectorDict{I, Int}(ks, vs), dualV) + return typeof(V)(mergewith(-, V.dims, W.dims), dualV) end function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} @@ -216,8 +214,7 @@ end function infimum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} Visdual = isdual(V₁) Visdual == isdual(V₂) || throw(SpaceMismatch("Infimum of space and dual space does not exist")) - ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, min, nothing, nothing) - return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual) + return typeof(V₁)(_sortedintersect(min, V₁.dims, V₂.dims), Visdual) end function supremum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} Visdual = isdual(V₁) @@ -228,8 +225,7 @@ end function supremum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} Visdual = isdual(V₁) Visdual == isdual(V₂) || throw(SpaceMismatch("Supremum of space and dual space does not exist")) - ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, max, identity, identity) - return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual) + return typeof(V₁)(mergewith(max, V₁.dims, V₂.dims), Visdual) end hassector(V::GradedSpace{I}, s::I) where {I <: Sector} = dim(V, s) != 0 From f8b0288cfd4d7bb11e51348f9f7a998530b62631 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Tue, 1 Sep 2026 18:04:01 +0200 Subject: [PATCH 19/36] introduce and use `blockdims` --- src/spaces/gradedspace.jl | 51 +++++++++++++++++++++++---------------- 1 file changed, 30 insertions(+), 21 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 25bce0ec5..e4bb53b21 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -31,7 +31,7 @@ sectortype(::Type{<:GradedSpace{I}}) where {I <: Sector} = I function GradedSpace{I, NTuple{N, Int}}(dims; dual::Bool = false) where {I, N} d = zeros(Int, N) - isset = falses(N) + isset = falses(N) # see if this is still needed if we're restricting to small N for (c, dc) in dims k = convert(I, c) i = findindex(values(I), k) @@ -111,9 +111,9 @@ end Base.axes(V::GradedSpace) = Base.OneTo(dim(V)) function Base.axes(V::GradedSpace{I}, c::I) where {I <: Sector} offset = 0 - for c′ in sectors(V) + for (c′, d′) in blockdims(V) c′ == c && break - offset += dim(c′) * dim(V, c′) + offset += dim(c′) * d′ end return (offset + 1):(offset + dim(c) * dim(V, c)) end @@ -124,9 +124,9 @@ isdual(V::GradedSpace) = V.dual isconj(V::GradedSpace) = isdual(V) function flip(V::GradedSpace{I}) where {I <: Sector} return if isdual(V) - typeof(V)(c => dim(V, c) for c in sectors(V)) + typeof(V)(blockdims(V)) else - typeof(V)(dual(c) => dim(V, c) for c in sectors(V))' + typeof(V)(dual(c) => d for (c, d) in blockdims(V))' end end @@ -160,20 +160,11 @@ function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) w end function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} - dual1, dual2 = isdual(V₁), isdual(V₂) - acc = Dict{I, Int}() # SectorDict `get` within the double for loop accumulates O(N^2) ` findindex` calls -> sort afterwards - k1, k2 = V₁.dims.keys, V₂.dims.keys - v1, v2 = V₁.dims.values, V₂.dims.values - @inbounds for na in eachindex(k1) - a₀, da = k1[na], v1[na] - a = dual1 ? dual(a₀) : a₀ - for nb in eachindex(k2) - b₀, db = k2[nb], v2[nb] - b = dual2 ? dual(b₀) : b₀ - dab = da * db - for c in a ⊗ b - acc[c] = get(acc, c, 0) + Nsymbol(a, b, c) * dab - end + acc = Dict{I, Int}() # SectorDict `get` within the double for loop accumulates O(N^2) `findindex` calls -> sort afterwards + for (a, da) in blockdims(V₁), (b, db) in blockdims(V₂) + dab = da * db + for c in a ⊗ b + acc[c] = get(acc, c, 0) + Nsymbol(a, b, c) * dab end end ks0 = collect(keys(acc)) @@ -238,6 +229,24 @@ function sectors(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} end end +""" + blockdims(V::GradedSpace) + +Return an iterator over the non-zero blocks of the graded space `V`. +These blocks contain the `Sector`s and their corresponding degeneracy, +i.e. the number of times the sector appears in the direct sum decomposition of `V`. +""" +function blockdims(V::GradedSpace{I, <:AbstractDict}) where {I <: Sector} + return ((isdual(V) ? dual(c) : c) => d for (c, d) in V.dims) +end +function blockdims(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} + vals = values(I) + return ( + (isdual(V) ? dual(vals[n]) : vals[n]) => V.dims[n] + for n in 1:N if !iszero(V.dims[n]) + ) +end + Base.hash(V::GradedSpace, h::UInt) = hash(V.dual, hash(V.dims, h)) function Base.:(==)(V₁::GradedSpace, V₂::GradedSpace) return sectortype(V₁) == sectortype(V₂) && (V₁.dims == V₂.dims) && V₁.dual == V₂.dual @@ -270,7 +279,7 @@ function Base.show(io::IO, V::GradedSpace) cls = ")" end - v = [c => dim(V, c) for c in sectors(V)] + v = collect(blockdims(V)) # logic stolen from Base.show_vector limited = get(io, :limit, false)::Bool @@ -301,7 +310,7 @@ function Base.show(io::IO, ::MIME"text/plain", V::GradedSpace) # print detailed sector information - hijack Base.Vector printing print(io, ":\n") isdual(V) && (V = dual(V)) - print_data = [c => dim(V, c) for c in sectors(V)] + print_data = collect(blockdims(V)) ioc = IOContext(io, :typeinfo => eltype(print_data)) Base.print_matrix(ioc, print_data) From 2fbf6feb2b8d7bd294ef80d1700b2591686f5be5 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 2 Sep 2026 17:06:56 +0200 Subject: [PATCH 20/36] introduce the ntuple cutoff --- src/spaces/gradedspace.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index e4bb53b21..c2f6758f9 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -330,19 +330,22 @@ Base.getindex(::SpaceTable) = ComplexSpace Base.getindex(::SpaceTable, ::Type{Trivial}) = ComplexSpace Base.getindex(::SpaceTable, I::Type{<:Sector}) = GradedSpace{I, sectorstoragetype(I)} +# based on Julia tuple unrolling range +const _ntuple_storage_threshold = 32 + """ sectorstoragetype(I::Type{<:Sector}) -> Type The storage type `D` used for the `dims` field of `GradedSpace{I, D}`. -This is `NTuple{N,Int}` with `N = length(values(I))` if `I` has a finite, known length, -or `SectorDict{I,Int}` otherwise. +This is `NTuple{N,Int}` with `N = length(values(I))` if `I` has a finite, known length +of at most `$_ntuple_storage_threshold`, or `SectorDict{I,Int}` otherwise. """ Base.@assume_effects :foldable function sectorstoragetype(::Type{I}) where {I <: Sector} if Base.IteratorSize(values(I)) isa Union{HasLength, HasShape} - return NTuple{length(values(I)), Int} - else - return SectorDict{I, Int} + N = length(values(I)) + N <= _ntuple_storage_threshold && return NTuple{N, Int} end + return SectorDict{I, Int} end Base.getindex(::ComplexNumbers, I::Type{<:Sector}) = Vect[I] From e5877dce1be4fd3023b95c2d236e0c4ca42d3eb1 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 4 Sep 2026 10:31:22 +0200 Subject: [PATCH 21/36] apply code suggestions --- src/auxiliary/dicts.jl | 21 ++++++--------------- src/spaces/gradedspace.jl | 13 ++----------- 2 files changed, 8 insertions(+), 26 deletions(-) diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl index c84a9e1e5..dbce97087 100644 --- a/src/auxiliary/dicts.jl +++ b/src/auxiliary/dicts.jl @@ -89,21 +89,12 @@ end Base.empty(::SortedVectorDict, ::Type{K}, ::Type{V}) where {K, V} = SortedVectorDict{K, V}() Base.empty!(d::SortedVectorDict) = (empty!(d.keys); empty!(d.values); return d) -_searchsortedfirst(v::Vector, k) = searchsortedfirst(v, k) -# function _searchsortedfirst(v::Vector, k) -# i = 1 -# @inbounds while i <= length(v) && isless(v[i], k) -# i += 1 -# end -# return i -# end - function Base.delete!(d::SortedVectorDict{K}, k) where {K} key = convert(K, k) if !isequal(k, key) return d end - i = _searchsortedfirst(d.keys, key) + i = searchsortedfirst(d.keys, key) if i <= length(d) && isequal(d.keys[i], key) deleteat!(d.keys, i) deleteat!(d.values, i) @@ -118,7 +109,7 @@ function Base.haskey(d::SortedVectorDict{K}, k) where {K} if !isequal(k, key) return false end - i = _searchsortedfirst(d.keys, key) + i = searchsortedfirst(d.keys, key) return (i <= length(d) && isequal(d.keys[i], key)) end function Base.getindex(d::SortedVectorDict{K}, k) where {K} @@ -126,7 +117,7 @@ function Base.getindex(d::SortedVectorDict{K}, k) where {K} if !isequal(k, key) throw(KeyError(k)) end - i = _searchsortedfirst(d.keys, key) + i = searchsortedfirst(d.keys, key) @inbounds if (i <= length(d) && isequal(d.keys[i], key)) return d.values[i] else @@ -138,7 +129,7 @@ function Base.setindex!(d::SortedVectorDict{K}, v, k) where {K} if !isequal(k, key) throw(ArgumentError("$k is not a valid key for type $K")) end - i = _searchsortedfirst(d.keys, key) + i = searchsortedfirst(d.keys, key) if i <= length(d) && isequal(d.keys[i], key) d.values[i] = v else @@ -153,7 +144,7 @@ function Base.get(d::SortedVectorDict{K}, k, default) where {K} if !isequal(k, key) return default end - i = _searchsortedfirst(d.keys, key) + i = searchsortedfirst(d.keys, key) @inbounds begin return (i <= length(d) && isequal(d.keys[i], key)) ? d.values[i] : default end @@ -163,7 +154,7 @@ function Base.get(f::Union{Function, Type}, d::SortedVectorDict{K}, k) where {K} if !isequal(k, key) return f() end - i = _searchsortedfirst(d.keys, key) + i = searchsortedfirst(d.keys, key) @inbounds begin return (i <= length(d) && isequal(d.keys[i], key)) ? d.values[i] : f() end diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index c2f6758f9..facb63811 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -89,18 +89,9 @@ GradedSpace(g::AbstractDict; dual::Bool = false) = GradedSpace(g...; dual = dual field(::Type{<:GradedSpace}) = ℂ InnerProductStyle(::Type{<:GradedSpace}) = EuclideanInnerProduct() -function dim(V::GradedSpace{I, <:AbstractDict}) where {I <: Sector} +function dim(V::GradedSpace{I}) where {I <: Sector} init = zero(dimscalartype(I)) - return sum(((c, d),) -> dim(c) * d, V.dims; init) -end -function dim(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} - init = zero(dimscalartype(I)) - D = init - vals = values(I) - @inbounds for (c, d) in zip(values(I), V.dims) - D += dim(c) * d - end - return D + return sum(((c, d),) -> dim(c) * d, blockdims(V); init) end function dim(V::GradedSpace{I, <:AbstractDict}, c::I) where {I <: Sector} return get(V.dims, isdual(V) ? dual(c) : c, 0) From 820d2229b75a7c1b5b5aa80fb3506f7388b4204d Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 8 Sep 2026 11:21:28 -0400 Subject: [PATCH 22/36] slight refactor of sortmerge implementation --- src/auxiliary/dicts.jl | 69 ++++++++++++++++++--------------------- src/spaces/gradedspace.jl | 2 +- 2 files changed, 33 insertions(+), 38 deletions(-) diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl index dbce97087..52d4d7695 100644 --- a/src/auxiliary/dicts.jl +++ b/src/auxiliary/dicts.jl @@ -177,63 +177,58 @@ function Base.:(==)(d1::SortedVectorDict, d2::SortedVectorDict) return true end -# merge over two SectorDicts -# the intersect case for infimum is kind of tricky, so there's an extra bool -# to indicate keeping keys that are only present in one of the two dicts -# zero results are dropped, matching how GradedSpace never stores an explicit zero dimension -function _sortedmerge( - combine, ::Val{keepunique}, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V} - ) where {keepunique, K, V} +# merge two SortedVectorDicts of `GradedSpace` dimensions, applying `combine` to keys present in +# both; keys present in only one dict are kept as is or dropped according to `_keepunmatched(combine)` +# zero results are dropped since `GradedSpace` never stores an explicit zero dimension +_keepunmatched(::Any) = true +_keepunmatched(::typeof(min)) = false # infimum: a missing sector has dimension zero, so min drops it + +function _sortedmerge(combine::F, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {F, K, V <: Integer} + keep = _keepunmatched(combine) k1, v1 = d1.keys, d1.values k2, v2 = d2.keys, d2.values n1, n2 = length(k1), length(k2) - ks, vs = Vector{K}(), Vector{V}() - sizehint!(ks, keepunique ? n1 + n2 : min(n1, n2)) - sizehint!(vs, keepunique ? n1 + n2 : min(n1, n2)) - i, j = 1, 1 + len = keep ? n1 + n2 : min(n1, n2) + ks = Vector{K}(undef, len) + vs = Vector{V}(undef, len) + i, j, n = 1, 1, 0 @inbounds while i <= n1 && j <= n2 - if k1[i] == k2[j] - d = combine(v1[i], v2[j]) - if !iszero(d) - push!(ks, k1[i]) - push!(vs, d) - end + a, b = k1[i], k2[j] + if isless(a, b) + keep && (n = _mergestore!(ks, vs, n, a, v1[i])) i += 1 + elseif isless(b, a) + keep && (n = _mergestore!(ks, vs, n, b, v2[j])) j += 1 - elseif k1[i] < k2[j] - if keepunique - push!(ks, k1[i]) - push!(vs, v1[i]) - end - i += 1 else - if keepunique - push!(ks, k2[j]) - push!(vs, v2[j]) - end + n = _mergestore!(ks, vs, n, a, combine(v1[i], v2[j])) + i += 1 j += 1 end end - if keepunique + if keep @inbounds while i <= n1 - push!(ks, k1[i]) - push!(vs, v1[i]) + n = _mergestore!(ks, vs, n, k1[i], v1[i]) i += 1 end @inbounds while j <= n2 - push!(ks, k2[j]) - push!(vs, v2[j]) + n = _mergestore!(ks, vs, n, k2[j], v2[j]) j += 1 end end + resize!(ks, n) + resize!(vs, n) return SortedVectorDict{K, V}(ks, vs) end +# write into slot `n + 1` and only advance the length when the value is nonzero +@inline function _mergestore!(ks, vs, n, k, d) + @inbounds ks[n + 1] = k + @inbounds vs[n + 1] = d + return n + !iszero(d) +end -Base.mergewith(combine, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {K, V} = - _sortedmerge(combine, Val(true), d1, d2) - -_sortedintersect(combine, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {K, V} = - _sortedmerge(combine, Val(false), d1, d2) +Base.mergewith(combine, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {K, V <: Integer} = + _sortedmerge(combine, d1, d2) """ Hashed(value, hashfunction = Base.hash, isequal = Base.isequal) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index facb63811..98403d5cf 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -196,7 +196,7 @@ end function infimum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} Visdual = isdual(V₁) Visdual == isdual(V₂) || throw(SpaceMismatch("Infimum of space and dual space does not exist")) - return typeof(V₁)(_sortedintersect(min, V₁.dims, V₂.dims), Visdual) + return typeof(V₁)(mergewith(min, V₁.dims, V₂.dims), Visdual) end function supremum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} Visdual = isdual(V₁) From 09e0d6f4cfae12bc13b3c6c52703eb492777a2b7 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 9 Sep 2026 11:24:21 +0200 Subject: [PATCH 23/36] introduce `FullVectorDict` and use it --- src/auxiliary/dicts.jl | 79 ++++++++++++++++++++++++++ src/factorizations/factorizations.jl | 2 +- src/factorizations/pullbacks.jl | 3 +- src/factorizations/truncation.jl | 85 +++++++++------------------- 4 files changed, 109 insertions(+), 60 deletions(-) diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl index 52d4d7695..de636b5cd 100644 --- a/src/auxiliary/dicts.jl +++ b/src/auxiliary/dicts.jl @@ -230,6 +230,85 @@ end Base.mergewith(combine, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {K, V <: Integer} = _sortedmerge(combine, d1, d2) +""" + FullVectorDict{K<:Sector,V} <: AbstractDict{K,V} + +Dictionary-like type that reserves one slot for every possible sector `c::K`, indexed +through `findindex(values(K), c)`. Absent entries are stored as `nothing`. +This is the counterpart of `SectorDict` for `NTuple{N,Int}`-backed storage. +Both are needed to build a "sector => value" map over some value type `V`. +""" +struct FullVectorDict{K <: Sector, V} <: AbstractDict{K, V} + slots::Vector{Union{Nothing, V}} + function FullVectorDict{K, V}(slots::Vector{Union{Nothing, V}}) where {K <: Sector, V} + @assert length(slots) == length(values(K)) + return new{K, V}(slots) + end +end +FullVectorDict{K, V}() where {K <: Sector, V} = + FullVectorDict{K, V}(Vector{Union{Nothing, V}}(nothing, length(values(K)))) +function FullVectorDict{K, V}(kv) where {K <: Sector, V} + d = FullVectorDict{K, V}() + for (k, v) in kv + d[k] = v + end + return d +end +FullVectorDict{K, V}(kv::Pair{K, V}...) where {K <: Sector, V} = FullVectorDict{K, V}(kv) + +Base.length(d::FullVectorDict) = count(!isnothing, d.slots) + +Base.copy(d::FullVectorDict{K, V}) where {K, V} = FullVectorDict{K, V}(copy(d.slots)) +Base.empty(::FullVectorDict{K}, ::Type{K}, ::Type{V}) where {K <: Sector, V} = FullVectorDict{K, V}() +Base.empty!(d::FullVectorDict) = (fill!(d.slots, nothing); return d) + +function Base.delete!(d::FullVectorDict{K}, k) where {K} + key = convert(K, k) + isequal(k, key) && (d.slots[findindex(values(K), key)] = nothing) + return d +end + +function Base.haskey(d::FullVectorDict{K}, k) where {K} + key = convert(K, k) + return isequal(k, key) && !isnothing(d.slots[findindex(values(K), key)]) +end +function Base.getindex(d::FullVectorDict{K}, k) where {K} + key = convert(K, k) + isequal(k, key) || throw(KeyError(k)) + v = d.slots[findindex(values(K), key)] + return isnothing(v) ? throw(KeyError(key)) : v +end +function Base.setindex!(d::FullVectorDict{K}, v, k) where {K} + key = convert(K, k) + isequal(k, key) || throw(ArgumentError("$k is not a valid key for type $K")) + d.slots[findindex(values(K), key)] = v + return d +end + +function Base.get(d::FullVectorDict{K}, k, default) where {K} + key = convert(K, k) + isequal(k, key) || return default + v = d.slots[findindex(values(K), key)] + return isnothing(v) ? default : v +end +function Base.get(f::Union{Function, Type}, d::FullVectorDict{K}, k) where {K} + key = convert(K, k) + isequal(k, key) || return f() + v = d.slots[findindex(values(K), key)] + return isnothing(v) ? f() : v +end + +function Base.iterate(d::FullVectorDict{K}, i = 1) where {K} + vals = values(K) + n = length(d.slots) + @inbounds while i <= n + v = d.slots[i] + isnothing(v) || return (vals[i] => v), i + 1 + i += 1 + end + return nothing +end + """ Hashed(value, hashfunction = Base.hash, isequal = Base.isequal) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index b16a28e9b..281070a77 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -6,7 +6,7 @@ module Factorizations export copy_oftype, factorisation_scalartype, one!, truncspace using ..TensorKit -using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, findindex, +using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, FullVectorDict, findindex, blocktype, foreachblock, one!, similar_diagonal, similarstoragetype, sectorstoragetype diff --git a/src/factorizations/pullbacks.jl b/src/factorizations/pullbacks.jl index 68ac57bb4..c8385107f 100644 --- a/src/factorizations/pullbacks.jl +++ b/src/factorizations/pullbacks.jl @@ -36,9 +36,8 @@ for pullback! in (:svd_pullback!, :eig_pullback!, :eigh_pullback!) Δt::AbstractTensorMap, t::AbstractTensorMap, F, ΔF, inds = _notrunc_ind(t); kwargs... ) - Isec = sectortype(t) foreachblock(Δt, t) do c, (Δb, b) - ind = _denseget(inds, Isec, c) + ind = get(inds, c, nothing) isnothing(ind) && return nothing Fc = block.(F, Ref(c)) ΔFc = block.(ΔF, Ref(c)) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index d2e455ebd..9d747532a 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -34,55 +34,30 @@ _blocklength(ax, ind) = length(ax[ind]) _blocklength(ax::Base.OneTo, ind::AbstractVector{<:Integer}) = length(ind) _blocklength(ax::Base.OneTo, ind::AbstractVector{Bool}) = count(ind) -# TODO: it quacks like a duck, just define a subtype of AbstractDict? -# represent the sector-index mapping as Vector{Union{Nothing, V}} where V is the type of the index -# mapping is indexed through findindex +# represent the sector-index mapping as either a FullVectorDict or SectorDict based on sectorstoragetype(I) # the type V is needed because the concrete type of ind depends on the strategy (except for intersect/union) -_densenew(::Type{I}, ::Type{V}) where {I <: Sector, V} = - Vector{Union{Nothing, V}}(nothing, length(values(I))) +_densemaptype(::Type{<:Tuple}, ::Type{I}, ::Type{V}) where {I <: Sector, V} = FullVectorDict{I, V} +_densemaptype(::Type{<:SectorDict}, ::Type{I}, ::Type{V}) where {I <: Sector, V} = SectorDict{I, V} -function _denseset!(v::Vector, ::Type{I}, c::I, val) where {I <: Sector} - v[findindex(values(I), c)] = val - return v -end -_denseget(v::Vector, ::Type{I}, c::I) where {I <: Sector} = v[findindex(values(I), c)] -function _densepairs(v::Vector, ::Type{I}) where {I <: Sector} - vals = values(I) - return (vals[i] => x for (i, x) in enumerate(v) if !isnothing(x)) -end -_densekeys(v::Vector, ::Type{I}) where {I <: Sector} = (c for (c, _) in _densepairs(v, I)) - -# fallbacks to catch SectorVector/SectorDict, even for NTuple sectorstoragetype -_denseget(v, ::Type{I}, c::I) where {I <: Sector} = get(v, c, nothing) -_densekeys(v, ::Type{I}) where {I <: Sector} = keys(v) -_densepairs(v, ::Type{I}) where {I <: Sector} = pairs(v) - -# builds either a dense Vector or SectorDict based on sectorstoragetype -# mapping each (c, v) pair's sector c to f(c, v) +# builds a sector => f(c, v) map through either FullVectorDict or SectorDict storage # so every `findtruncated` method shares one output-construction path # `pairsiter` are c => v pairs, can be c => nothing for NoTruncation/TruncationIntersection/TruncationUnion -function _builddensemap(f, ::Type{D}, ::Type{I}, pairsiter, ::Type{V}) where {D <: Tuple, I <: Sector, V} - d = _densenew(I, V) - for (c, v) in pairsiter - _denseset!(d, I, c, f(c, v)) - end - return d -end -function _builddensemap(f, ::Type{D}, ::Type{I}, pairsiter, ::Type{V}) where {D <: SectorDict, I <: Sector, V} - return SectorDict(c => f(c, v) for (c, v) in pairsiter) # V unused +function _builddensemap(f, ::Type{D}, ::Type{I}, pairsiter, ::Type{V}) where {D, I <: Sector, V} + Dout = _densemaptype(D, I, V) + return Dout(c => f(c, v) for (c, v) in pairsiter) end function truncate_space(V::ElementarySpace, inds) @assert !isdual(V) I = sectortype(V) @assert I == Trivial - return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in _densepairs(inds, I)) + return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in pairs(inds)) end function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N} @assert !isdual(V) vals = values(I) newdims = zeros(Int, N) - for (c, ind) in _densepairs(inds, I) + for (c, ind) in pairs(inds) d = dim(V, c) n_write = findindex(vals, c) newdims[n_write] = _blocklength(d, ind) @@ -105,30 +80,27 @@ function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Secto end function truncate_domain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds) - Isec = sectortype(tdst) for (c, b) in blocks(tdst) - I = _denseget(inds, Isec, c) - @assert !isnothing(I) # kept for safety, but should be guaranteed by _densepairs + I = get(inds, c, nothing) + @assert !isnothing(I) # kept for safety, but should be guaranteed by pairs(inds) covering every block sector b′ = block(tsrc, c) b .= view(b′, :, I) end return tdst end function truncate_codomain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds) - Isec = sectortype(tdst) for (c, b) in blocks(tdst) - I = _denseget(inds, Isec, c) - @assert !isnothing(I) # kept for safety, but should be guaranteed by _densepairs + I = get(inds, c, nothing) + @assert !isnothing(I) # kept for safety, but should be guaranteed by pairs(inds) covering every block sector b′ = block(tsrc, c) b .= view(b′, I, :) end return tdst end function truncate_diagonal!(Ddst::DiagonalTensorMap, Dsrc::DiagonalTensorMap, inds) - Isec = sectortype(Ddst) for (c, b) in blocks(Ddst) - I = _denseget(inds, Isec, c) - @assert !isnothing(I) # kept for safety, but should be guaranteed by _densepairs + I = get(inds, c, nothing) + @assert !isnothing(I) # kept for safety, but should be guaranteed by pairs(inds) covering every block sector diagview(b) .= view(diagview(block(Dsrc, c)), I) end return Ddst @@ -375,37 +347,37 @@ end function MAK.findtruncated(values::SectorVector, strategy::TruncationIntersection) I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated, values), strategy.components) - @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now" - sectors = collect(_densekeys(first(inds), I)) + @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" + sectors = collect(keys(first(inds))) return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ - mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_intersect, inds) + mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_intersect, inds) end end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationIntersection) I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated_svd, values), strategy.components) - @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now" - sectors = collect(_densekeys(first(inds), I)) + @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" + sectors = collect(keys(first(inds))) return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ - mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_intersect, inds) + mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_intersect, inds) end end function MAK.findtruncated(values::SectorVector, strategy::TruncationUnion) I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated, values), strategy.components) - @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now" - sectors = collect(_densekeys(first(inds), I)) + @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" + sectors = collect(keys(first(inds))) return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ - mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_union, inds) + mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_union, inds) end end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationUnion) I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated_svd, values), strategy.components) - @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now" - sectors = collect(_densekeys(first(inds), I)) + @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" + sectors = collect(keys(first(inds))) return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ - mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_union, inds) + mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_union, inds) end end @@ -414,8 +386,7 @@ end MAK.truncation_error(values::SectorVector, ind) = MAK.truncation_error!(copy(values), ind) function MAK.truncation_error!(values::SectorVector, ind) - Isec = sectortype(values) - for (c, ind_c) in _densepairs(ind, Isec) + for (c, ind_c) in pairs(ind) v = values[c] v[ind_c] .= zero(eltype(v)) end From 040b8c3259f1fc8da2ab34260441d0c22ce46775 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Wed, 9 Sep 2026 11:31:59 +0200 Subject: [PATCH 24/36] remove some reminder comments --- src/spaces/gradedspace.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 98403d5cf..3c02726e1 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -31,7 +31,7 @@ sectortype(::Type{<:GradedSpace{I}}) where {I <: Sector} = I function GradedSpace{I, NTuple{N, Int}}(dims; dual::Bool = false) where {I, N} d = zeros(Int, N) - isset = falses(N) # see if this is still needed if we're restricting to small N + isset = falses(N) for (c, dc) in dims k = convert(I, c) i = findindex(values(I), k) @@ -151,7 +151,7 @@ function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) w end function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} - acc = Dict{I, Int}() # SectorDict `get` within the double for loop accumulates O(N^2) `findindex` calls -> sort afterwards + acc = Dict{I, Int}() for (a, da) in blockdims(V₁), (b, db) in blockdims(V₂) dab = da * db for c in a ⊗ b @@ -170,12 +170,12 @@ function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{ @inbounds for na in 1:N da = V₁.dims[na] iszero(da) && continue - a₀ = vals[na] # avoid call to sectors(V₁) + a₀ = vals[na] a = dual1 ? dual(a₀) : a₀ for nb in 1:N db = V₂.dims[nb] iszero(db) && continue - b₀ = vals[nb] # idem for V₂ + b₀ = vals[nb] b = dual2 ? dual(b₀) : b₀ dab = da * db for c in a ⊗ b From b653bbf617cd52fbaa533056068dbef23880773c Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 9 Sep 2026 14:51:18 -0400 Subject: [PATCH 25/36] avoid `NTuple{N,Int}(::Vector)` conversion cliff at N == 32 `NTuple{N,Int}(v)` falls off Julia's unrolled fast path at exactly N = 32: 9.7 ns at N = 31 versus 825 ns at N = 32. That is precisely `_ntuple_storage_threshold`, so the largest sector sets that tuple storage is meant to cover were by far the slowest to build. `ntuple(i -> @inbounds(v[i]), Val(N))` stays on the unrolled path and is faster at every N (3.1 ns at N = 32), so the threshold can stay at 32. Measured end-to-end, before -> after: Vect[ZNIrrep{32}] constructor 1067 ns -> 137 ns fuse, Z32, full 1329 ns -> 342 ns, 864 B/4 allocs -> 0 fuse, Z32, sparse 1088 ns -> 134 ns, 864 B/4 allocs -> 0 Co-Authored-By: Claude Opus 5 (1M context) --- src/factorizations/truncation.jl | 2 +- src/spaces/gradedspace.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 9d747532a..0a88307b1 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -62,7 +62,7 @@ function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sec n_write = findindex(vals, c) newdims[n_write] = _blocklength(d, ind) end - return typeof(V)(NTuple{N, Int}(newdims), false) + return typeof(V)(ntuple(i -> @inbounds(newdims[i]), Val(N)), false) end function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector} @assert !isdual(V) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 3c02726e1..d67e016ec 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -40,7 +40,7 @@ function GradedSpace{I, NTuple{N, Int}}(dims; dual::Bool = false) where {I, N} isset[i] = true d[i] = dc end - return GradedSpace{I, NTuple{N, Int}}(NTuple{N, Int}(d), dual) + return GradedSpace{I, NTuple{N, Int}}(ntuple(i -> @inbounds(d[i]), Val(N)), dual) end function GradedSpace{I, NTuple{N, Int}}(dims::Pair; dual::Bool = false) where {I, N} return GradedSpace{I, NTuple{N, Int}}((dims,); dual = dual) @@ -184,7 +184,7 @@ function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{ end end end - return typeof(V₁)(NTuple{N, Int}(newdims), false) + return typeof(V₁)(ntuple(i -> @inbounds(newdims[i]), Val(N)), false) end function infimum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} From 8808e7662381911b339b4ba3f4f74a9ac36da249 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 9 Sep 2026 14:51:38 -0400 Subject: [PATCH 26/36] build tuple storage in a `MutableNTuple` with a bitmask The scatter-write paths that build `NTuple{N,Int}` storage used two heap buffers: a `zeros(Int, N)` accumulator and a `falses(N)` flag vector, both only there because `Base.setindex` on a tuple in a loop is slow. `TupleTools.MutableNTuple` replaces the accumulator: the optimizer elides the mutable box in these call patterns, so nothing escapes to the heap. Tracking which sectors have already been set needs one bit per sector, and tuple storage is only used for `N <= _ntuple_storage_threshold`, so a `UInt64` mask does that with no allocation at all. `Base.setindex` was measured as the third option and ruled out: 2-3x behind at N = 16 and 11x behind on the `fuse` accumulation pattern at N = 31. Constructor, before -> after: N = 4 121 ns, 240 B, 6 allocs -> 80 ns, 48 B, 1 alloc N = 16 115 ns, 432 B, 6 allocs -> 86 ns, 144 B, 1 alloc N = 32 137 ns, 560 B, 6 allocs -> 102 ns, 272 B, 1 alloc `fuse` deliberately keeps its plain `Vector` accumulator: it does not escape either, and `MutableNTuple` is ~9% slower there at N = 16 because the pointer round-trip blocks register promotion of the accumulator. Co-Authored-By: Claude Opus 5 (1M context) --- src/factorizations/factorizations.jl | 2 ++ src/factorizations/truncation.jl | 6 +++--- src/spaces/gradedspace.jl | 12 +++++++----- 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index dac48a6f7..7a58e6d58 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -18,6 +18,8 @@ using LinearAlgebra: LinearAlgebra, BlasFloat, Diagonal, using TensorOperations: Index2Tuple +using TupleTools: MutableNTuple, StaticLength + using MatrixAlgebraKit import MatrixAlgebraKit as MAK using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, DiagonalAlgorithm diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 0a88307b1..9b0d3836f 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -56,13 +56,13 @@ end function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N} @assert !isdual(V) vals = values(I) - newdims = zeros(Int, N) + newdims = MutableNTuple(ntuple(_ -> 0, StaticLength(N))) for (c, ind) in pairs(inds) d = dim(V, c) n_write = findindex(vals, c) - newdims[n_write] = _blocklength(d, ind) + @inbounds newdims[n_write] = _blocklength(d, ind) end - return typeof(V)(ntuple(i -> @inbounds(newdims[i]), Val(N)), false) + return typeof(V)(Tuple(newdims), false) end function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector} @assert !isdual(V) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index d67e016ec..584b52f70 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -30,17 +30,19 @@ end sectortype(::Type{<:GradedSpace{I}}) where {I <: Sector} = I function GradedSpace{I, NTuple{N, Int}}(dims; dual::Bool = false) where {I, N} - d = zeros(Int, N) - isset = falses(N) + @assert N <= 64 "the `UInt64` bitmask tracking which sectors have been set holds 64 bits" + d = TupleTools.MutableNTuple(ntuple(Returns(0), StaticLength(N))) + mask = zero(UInt64) for (c, dc) in dims k = convert(I, c) i = findindex(values(I), k) dc < 0 && throw(ArgumentError(lazy"Sector $k has negative dimension $dc")) - isset[i] && throw(ArgumentError(lazy"Sector $c appears multiple times")) - isset[i] = true + bit = one(UInt64) << (i - 1) + iszero(mask & bit) || throw(ArgumentError(lazy"Sector $c appears multiple times")) + mask |= bit d[i] = dc end - return GradedSpace{I, NTuple{N, Int}}(ntuple(i -> @inbounds(d[i]), Val(N)), dual) + return GradedSpace{I, NTuple{N, Int}}(Tuple(d), dual) end function GradedSpace{I, NTuple{N, Int}}(dims::Pair; dual::Bool = false) where {I, N} return GradedSpace{I, NTuple{N, Int}}((dims,); dual = dual) From dd3da4e3781e8d72d79a8b2b2ec5c16b86ebca64 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 9 Sep 2026 14:52:01 -0400 Subject: [PATCH 27/36] specialise `flip` to tuple storage The generic `flip` rebuilds the space through `GradedSpace{I,NTuple{N,Int}}(dims)`, paying a `findindex` per sector plus the whole constructor. For tuple storage the result is just a fixed permutation of the dims with `isdual` negated, and that permutation depends only on `I`, so `Base.@assume_effects :foldable` lets it fold away entirely. flip, Z32 1083 ns, 1280 B, 9 allocs -> 5.2 ns, 0 B, 0 allocs flip, Z4 x Z8 1950 ns, 1280 B, 9 allocs -> 4.9 ns, 0 B, 0 allocs Verified to agree with the generic method on both `isdual` branches for Z4, Z8, Z32, a product sector, and empty spaces. Co-Authored-By: Claude Opus 5 (1M context) --- src/spaces/gradedspace.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 584b52f70..a3bd3308e 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -123,6 +123,18 @@ function flip(V::GradedSpace{I}) where {I <: Sector} end end +# the permutation of `values(I)` induced by `dual`; only depends on the type, so it folds away +Base.@assume_effects :foldable function _dualpermutation(::Type{I}, ::Val{N}) where {I <: Sector, N} + vals = values(I) + return ntuple(n -> findindex(vals, dual(vals[n])), Val(N)) +end +function flip(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} + # `flip` maps `c => d` to `dual(c) => d` and negates `isdual`, which for tuple storage is + # just a fixed permutation of the dims, so the generic constructor can be skipped + newdims = TupleTools.getindices(V.dims, _dualpermutation(I, Val(N))) + return GradedSpace{I, NTuple{N, Int}}(newdims, !isdual(V)) +end + function unitspace(S::Type{<:GradedSpace{I}}) where {I <: Sector} return S(unit => 1 for unit in allunits(I)) end From 4855f8fa0e5d07cf93c37b5e11f768028aac3d51 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 9 Sep 2026 14:52:48 -0400 Subject: [PATCH 28/36] make `ominus` a single pass MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `⊖` validated with `V ≿ W` and only then subtracted the dimensions, so every call traversed the sectors twice. `isepimorphic` also `collect`s `blocksectors(W)` and does a `dim` lookup per sector, which is what made `⊖` allocate on tuple storage where the subtraction itself is allocation-free. Both methods now subtract and validate in one pass: `map` over the two tuples, or one walk over the two sorted key vectors reusing `_mergestore!` so zero results are dropped. This also closes a latent trap. `_keepunmatched(-)` fell through to the `true` default, so `mergewith(-, ...)` gave a sector present only in `W` a *positive* dimension in the result; only the `≿` pre-check kept that unreachable. `⊖` was the sole `mergewith(-, ...)` caller, and there are none left. Rep[Z4] (tuple, N=4) 18.8 ns, 64 B, 2 allocs -> 2.1 ns, 0 B, 0 allocs ZNIrrep{8} (tuple, N=8) 42.8 ns, 112 B, 3 allocs -> 2.4 ns, 0 B, 0 allocs Rep[U1] (dict, 5 sect.) 76.3 ns, 384 B, 6 allocs -> 31.1 ns, 224 B, 5 allocs Rep[SU2] (dict, 6 sect.) 59.3 ns, 432 B, 6 allocs -> 27.9 ns, 256 B, 5 allocs Co-Authored-By: Claude Opus 5 (1M context) --- src/spaces/gradedspace.jl | 36 ++++++++++++++++++++++++++++++++---- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index a3bd3308e..550bff781 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -152,16 +152,44 @@ function ⊕(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where newdims = map(+, V₁.dims, V₂.dims) return typeof(V₁)(newdims, dual1) end +@noinline _throw_not_subspace(V, W) = throw(SpaceMismatch(lazy"$(W) is not a subspace of $(V)")) function ⊖(V::GradedSpace{I, <:Tuple}, W::GradedSpace{I, <:Tuple}) where {I <: Sector} dualV = isdual(V) - V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)")) - newdims = map(-, V.dims, W.dims) + dualV == isdual(W) || _throw_not_subspace(V, W) + # single unrolled pass: subtract and validate non-negativity at once + newdims = map((dV, dW) -> dV < dW ? _throw_not_subspace(V, W) : dV - dW, V.dims, W.dims) return typeof(V)(newdims, dualV) end function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) where {I <: Sector} dualV = isdual(V) - V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)")) - return typeof(V)(mergewith(-, V.dims, W.dims), dualV) + dualV == isdual(W) || _throw_not_subspace(V, W) + k1, v1 = keys(V.dims), values(V.dims) + k2, v2 = keys(W.dims), values(W.dims) + n1, n2 = length(k1), length(k2) + ks, vs = Vector{I}(undef, n1), Vector{Int}(undef, n1) + i, j, n = 1, 1, 0 + @inbounds while i <= n1 && j <= n2 + a, b = k1[i], k2[j] + if isless(a, b) + n = _mergestore!(ks, vs, n, a, v1[i]) + i += 1 + elseif isless(b, a) + _throw_not_subspace(V, W) # sector of `W` absent from `V` + else + v1[i] < v2[j] && _throw_not_subspace(V, W) + n = _mergestore!(ks, vs, n, a, v1[i] - v2[j]) + i += 1 + j += 1 + end + end + j <= n2 && _throw_not_subspace(V, W) # leftover sectors of `W` absent from `V` + @inbounds while i <= n1 + n = _mergestore!(ks, vs, n, k1[i], v1[i]) + i += 1 + end + resize!(ks, n) + resize!(vs, n) + return typeof(V)(SectorDict{I, Int}(ks, vs), dualV) end function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} From 7d474d70bf8dea897215ff2622ad6f209bfa051d Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 9 Sep 2026 14:53:02 -0400 Subject: [PATCH 29/36] replace `_builddensemap` by `sectormap` in dicts.jl `_builddensemap(f, D, I, pairsiter, V)` took five arguments, two of which (`sectorstoragetype(I)` and `I`) were derivable and were spelled out at all fourteen call sites. It is really a `map` over `c => v` pairs into whichever dense sector map matches the storage choice for `GradedSpace{I}`, so it belongs next to the dict types rather than in the truncation code. - `sectormaptype(I, V)` folds the two-step storage dispatch into one function, keeping `_densemaptype` as the internal dispatch on the storage type. - `sectormap(f, [V,] pairsiter)` derives `I` from `eltype(pairsiter)` and, when `V` is omitted, the value type from `f` via `promote_op` - which is what four call sites were computing by hand. Only the intersection/union sites still pass `V` explicitly, to pin it to `Any`. - The separate keys-only builder is gone; callers that need just one half of the pair take a dummy argument. `pairs` normalises the intersection/union inputs, whose components may be `SectorVector` masks rather than dicts. `pairs(::SectorVector)` materialises a `SectorDict` that `sectormap` immediately re-iterates, but making it lazy would return a `Base.Generator` with `eltype == Any` and break the sector-type deduction at every call site, so it is left alone. Co-Authored-By: Claude Opus 5 (1M context) --- src/auxiliary/dicts.jl | 32 +++++++++++++++ src/factorizations/factorizations.jl | 4 +- src/factorizations/pullbacks.jl | 7 +--- src/factorizations/truncation.jl | 59 ++++++---------------------- 4 files changed, 47 insertions(+), 55 deletions(-) diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl index de636b5cd..61ae8c99f 100644 --- a/src/auxiliary/dicts.jl +++ b/src/auxiliary/dicts.jl @@ -309,6 +309,38 @@ function Base.iterate(d::FullVectorDict{K}, i = 1) where {K} return nothing end +# the dense sector map is backed by a `FullVectorDict` or a `SortedVectorDict`, mirroring the +# `NTuple`/`SortedVectorDict` storage choice that `sectorstoragetype` makes for `GradedSpace` +_densemaptype(::Type{<:Tuple}, ::Type{I}, ::Type{V}) where {I <: Sector, V} = FullVectorDict{I, V} +_densemaptype(::Type{<:SortedVectorDict}, ::Type{I}, ::Type{V}) where {I <: Sector, V} = + SortedVectorDict{I, V} + +""" + sectormaptype(::Type{I}, ::Type{V}) where {I <: Sector, V} -> Type + +The dense `AbstractDict{I, V}` type used to map sectors of type `I` onto values of type `V`, +chosen to match the storage type of `GradedSpace{I}`, see [`sectorstoragetype`](@ref). +""" +sectormaptype(::Type{I}, ::Type{V}) where {I <: Sector, V} = + _densemaptype(sectorstoragetype(I), I, V) + +""" + sectormap(f, [V::Type, ] pairsiter) + +Map `f(c, v)` over an iterator `pairsiter` of `c => v` pairs, collecting the results into the +dense sector map returned by [`sectormaptype`](@ref). The value type `V` of the result is +inferred from `f` unless it is given explicitly. +""" +sectormap(f, pairsiter) = _sectormap(f, pairsiter, eltype(pairsiter)) +sectormap(f, ::Type{V}, pairsiter) where {V} = _sectormap(f, V, pairsiter, eltype(pairsiter)) + +function _sectormap(f, pairsiter, ::Type{<:Pair{I, W}}) where {I <: Sector, W} + return _sectormap(f, Base.promote_op(f, I, W), pairsiter, Pair{I, W}) +end +function _sectormap(f, ::Type{V}, pairsiter, ::Type{<:Pair{I}}) where {V, I <: Sector} + return sectormaptype(I, V)(c => f(c, v) for (c, v) in pairsiter) +end + """ Hashed(value, hashfunction = Base.hash, isequal = Base.isequal) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 7a58e6d58..9b214b776 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -6,9 +6,9 @@ module Factorizations export copy_oftype, factorisation_scalartype, one!, truncspace using ..TensorKit -using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, FullVectorDict, findindex, +using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, findindex, blocktype, foreachblock, one!, - similar_diagonal, similarstoragetype, sectorstoragetype + similar_diagonal, similarstoragetype, sectormap using ..TensorKit: GLOBAL_TIMER using TimerOutputs: @timeit_debug diff --git a/src/factorizations/pullbacks.jl b/src/factorizations/pullbacks.jl index c8385107f..b63bc9a94 100644 --- a/src/factorizations/pullbacks.jl +++ b/src/factorizations/pullbacks.jl @@ -24,12 +24,7 @@ for pullback! in (:qr_null_pullback!, :lq_null_pullback!) return Δt end end -function _notrunc_ind(t) - I = sectortype(t) - return _builddensemap(sectorstoragetype(I), I, blocks(t), Colon) do _, _ - Colon() - end -end +_notrunc_ind(t) = sectormap(Returns(Colon()), blocks(t)) for pullback! in (:svd_pullback!, :eig_pullback!, :eigh_pullback!) @eval function MAK.$pullback!( diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 9b0d3836f..71bb54b68 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -34,19 +34,6 @@ _blocklength(ax, ind) = length(ax[ind]) _blocklength(ax::Base.OneTo, ind::AbstractVector{<:Integer}) = length(ind) _blocklength(ax::Base.OneTo, ind::AbstractVector{Bool}) = count(ind) -# represent the sector-index mapping as either a FullVectorDict or SectorDict based on sectorstoragetype(I) -# the type V is needed because the concrete type of ind depends on the strategy (except for intersect/union) -_densemaptype(::Type{<:Tuple}, ::Type{I}, ::Type{V}) where {I <: Sector, V} = FullVectorDict{I, V} -_densemaptype(::Type{<:SectorDict}, ::Type{I}, ::Type{V}) where {I <: Sector, V} = SectorDict{I, V} - -# builds a sector => f(c, v) map through either FullVectorDict or SectorDict storage -# so every `findtruncated` method shares one output-construction path -# `pairsiter` are c => v pairs, can be c => nothing for NoTruncation/TruncationIntersection/TruncationUnion -function _builddensemap(f, ::Type{D}, ::Type{I}, pairsiter, ::Type{V}) where {D, I <: Sector, V} - Dout = _densemaptype(D, I, V) - return Dout(c => f(c, v) for (c, v) in pairsiter) -end - function truncate_space(V::ElementarySpace, inds) @assert !isdual(V) I = sectortype(V) @@ -155,8 +142,7 @@ end function MAK.truncate( ::typeof(left_null!), (U, S)::NTuple{2, AbstractTensorMap}, strategy::NoTruncation ) - I = sectortype(S) - ind = _builddensemap(sectorstoragetype(I), I, blocks(S), UnitRange{Int}) do _, b + ind = sectormap(blocks(S)) do _, b (size(b, 2) + 1):size(b, 1) end V_truncated = truncate_space(space(S, 1), ind) @@ -167,8 +153,7 @@ end function MAK.truncate( ::typeof(right_null!), (S, Vᴴ)::NTuple{2, AbstractTensorMap}, strategy::NoTruncation ) - I = sectortype(S) - ind = _builddensemap(sectorstoragetype(I), I, blocks(S), UnitRange{Int}) do _, b + ind = sectormap(blocks(S)) do _, b (size(b, 1) + 1):size(b, 2) end V_truncated = truncate_space(dual(space(S, 2)), ind) @@ -207,10 +192,7 @@ function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationStrateg end function MAK.findtruncated(values::SectorVector, ::NoTruncation) - I = sectortype(values) - return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in keys(values)), Colon) do _, _ - Colon() - end + return sectormap(Returns(Colon()), pairs(values)) end # Need to select the first k values here after sorting across blocks, weighted by quantum dimension @@ -256,27 +238,22 @@ MAK.findtruncated_svd(values::SectorVector, strategy::TruncationByOrder) = MAK.findtruncated(values, strategy) function MAK.findtruncated(values::SectorVector, strategy::TruncationByFilter) - I = sectortype(values) - return _builddensemap(sectorstoragetype(I), I, pairs(values), Vector{Int}) do _, v + return sectormap(pairs(values)) do _, v findall(strategy.filter, v) end end function MAK.findtruncated(values::SectorVector, strategy::TruncationByValue) - I = sectortype(values) atol = rtol_to_atol(values, strategy.p, strategy.atol, strategy.rtol) strategy′ = trunctol(; atol, strategy.by, strategy.keep_below) - V = Base.promote_op(MAK.findtruncated, valtype(values), typeof(strategy′)) - return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do _, v + return sectormap(pairs(values)) do _, v MAK.findtruncated(v, strategy′) end end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationByValue) - I = sectortype(values) atol = rtol_to_atol(values, strategy.p, strategy.atol, strategy.rtol) strategy′ = trunctol(; atol, strategy.by, strategy.keep_below) - V = Base.promote_op(MAK.findtruncated_svd, valtype(values), typeof(strategy′)) - return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do _, v + return sectormap(pairs(values)) do _, v MAK.findtruncated_svd(v, strategy′) end end @@ -323,9 +300,7 @@ function MAK.findtruncated(values::SectorVector, strategy::TruncationSpace) I = sectortype(values) I == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values")) blockstrategy(c) = truncrank(dim(strategy.space, c); strategy.by, strategy.rev) - Vstrategy = Base.promote_op(blockstrategy, I) - V = Base.promote_op(MAK.findtruncated, valtype(values), Vstrategy) - return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do c, v + return sectormap(pairs(values)) do c, v MAK.findtruncated(v, blockstrategy(c)) end end @@ -333,9 +308,7 @@ function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationSpace) I = sectortype(values) I == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values")) blockstrategy(c) = truncrank(dim(strategy.space, c); strategy.by, strategy.rev) - Vstrategy = Base.promote_op(blockstrategy, I) - V = Base.promote_op(MAK.findtruncated_svd, valtype(values), Vstrategy) - return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do c, v + return sectormap(pairs(values)) do c, v MAK.findtruncated_svd(v, blockstrategy(c)) end end @@ -345,38 +318,30 @@ end # This is always the case in the implementations above. function MAK.findtruncated(values::SectorVector, strategy::TruncationIntersection) - I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated, values), strategy.components) @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" - sectors = collect(keys(first(inds))) - return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ + return sectormap(Any, pairs(first(inds))) do c, _ mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_intersect, inds) end end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationIntersection) - I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated_svd, values), strategy.components) @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" - sectors = collect(keys(first(inds))) - return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ + return sectormap(Any, pairs(first(inds))) do c, _ mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_intersect, inds) end end function MAK.findtruncated(values::SectorVector, strategy::TruncationUnion) - I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated, values), strategy.components) @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" - sectors = collect(keys(first(inds))) - return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ + return sectormap(Any, pairs(first(inds))) do c, _ mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_union, inds) end end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationUnion) - I = sectortype(values) inds = map(Base.Fix1(MAK.findtruncated_svd, values), strategy.components) @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" - sectors = collect(keys(first(inds))) - return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _ + return sectormap(Any, pairs(first(inds))) do c, _ mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_union, inds) end end From 2779e2764bb2fa0adb1f68d4686f78b9d3c5f326 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 9 Sep 2026 14:54:26 -0400 Subject: [PATCH 30/36] document why `fuse` accumulates into a `Dict` Inserting straight into the sorted result was measured as the alternative. It wins below a few hundred output sectors but is O(k^2) in them, and product sectors reach that regime quickly: U1 x U1 with 196 input sectors fuses to 729 output sectors at 1.03 ms against 0.57 ms for the `Dict`, and the gap grows. The `Dict` has no such cliff, so it stays. Also drops a stray blank line left before `_dualpermutation`. Co-Authored-By: Claude Opus 5 (1M context) --- src/spaces/gradedspace.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 550bff781..2745439d4 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -122,7 +122,6 @@ function flip(V::GradedSpace{I}) where {I <: Sector} typeof(V)(dual(c) => d for (c, d) in blockdims(V))' end end - # the permutation of `values(I)` induced by `dual`; only depends on the type, so it folds away Base.@assume_effects :foldable function _dualpermutation(::Type{I}, ::Val{N}) where {I <: Sector, N} vals = values(I) @@ -193,7 +192,7 @@ function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) w end function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} - acc = Dict{I, Int}() + acc = Dict{I, Int}() # Accumulation into Dict is more efficient than repeated insertion in sorted vector for (a, da) in blockdims(V₁), (b, db) in blockdims(V₂) dab = da * db for c in a ⊗ b From aa04e15373f9fa5549925be07c8e7ed1e1cbfce6 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 14 Sep 2026 20:46:44 -0400 Subject: [PATCH 31/36] address review: drop `FullVectorDict` and the `sectormap` paradigm MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `FullVectorDict`/`sectormap` existed only to give the truncation code a dense "sector => index" map for tuple-backed spaces. Their call sites all go through `pairs(::SectorVector)`, which materialises a `SectorDict` anyway, so the dense map bought little for a lot of `Sector`-specific machinery in `auxiliary/`. Every truncation index map is a plain `SectorDict` again; instead `pairs(::SectorVector)` is built directly from the already-sorted structure rather than by repeated sorted insertion. The storage-specialised `truncate_space` methods are kept, since they never used the dense map. Also from the review: - `⊖` for dict storage reuses `_sortedmerge`, whose `_keepunmatched` trait is generalised into per-side `_unmatched1`/`_unmatched2` hooks - `_ntuple_storage_threshold` -> `_NTUPLE_STORAGE_THRESHOLD`, lowered to 8 - `ZNSpace{N}` deprecated in favour of `Vect[ZNIrrep{N}]`, which the alias can no longer track once `N` exceeds the threshold - drop the over-strong `I == Trivial` assert in `truncate_space`, and name the index variables `ind`/`inds` instead of shadowing the sector type `I` Co-Authored-By: Claude Opus 5 (1M context) --- docs/src/Changelog.md | 2 + docs/src/man/gradedspaces.md | 7 +- src/auxiliary/dicts.jl | 150 ++++----------------------- src/factorizations/factorizations.jl | 2 +- src/factorizations/pullbacks.jl | 2 +- src/factorizations/truncation.jl | 138 ++++++++++++------------ src/spaces/gradedspace.jl | 93 +++++++++-------- src/tensors/sectorvector.jl | 6 +- 8 files changed, 146 insertions(+), 254 deletions(-) diff --git a/docs/src/Changelog.md b/docs/src/Changelog.md index 9a6b231b9..fbbab7cb3 100644 --- a/docs/src/Changelog.md +++ b/docs/src/Changelog.md @@ -26,12 +26,14 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - For sector types with `GenericUnit` such that colorings are not unique, `GradedSpace`, `ProductSpace` and `HomSpace` now check for this compatibility. In particular, this prevents the construction of `TensorMap`s with incompatible colorings, which previously either errored or produced empty tensors inconsistently. ([#515](https://github.com/QuantumKitHub/TensorKit.jl/pull/515)) ### Deprecated +- The type alias `ZNSpace{N}` is deprecated in favour of `Vect[ZNIrrep{N}]`: a type alias cannot compute the storage type from `N`, so the two only agree for small `N`. ([#511](https://github.com/QuantumKitHub/TensorKit.jl/pull/511)) ### Removed ### Fixed ### Performance +- `GradedSpace` operations (`dim`, `flip`, `⊕`, `⊖`, `fuse`, `infimum`, `supremum`, truncation) are now specialised on the storage type of the degeneracy dimensions, and tuple storage is used only for sector types with at most `TensorKit._NTUPLE_STORAGE_THRESHOLD` sectors so that sector types with many sectors no longer burden the compiler. ([#511](https://github.com/QuantumKitHub/TensorKit.jl/pull/511)) ## [0.17.1](https://github.com/QuantumKitHub/TensorKit.jl/compare/v0.17.0...v0.17.1) - 2026-07-13 diff --git a/docs/src/man/gradedspaces.md b/docs/src/man/gradedspaces.md index 69342dcba..529586cf3 100644 --- a/docs/src/man/gradedspaces.md +++ b/docs/src/man/gradedspaces.md @@ -25,16 +25,17 @@ However, this is mostly to lower the barrier, as really the instances of `Graded ## Implementation details As mentioned, the way in which the degeneracy dimensions ``n_a`` are stored depends on the specific sector type `I`, more specifically on the `IteratorSize` of `values(I)`. -If `IteratorSize(values(I)) isa Union{IsInfinite, SizeUnknown}`, the different sectors ``a`` and their corresponding degeneracy ``n_a`` are stored as key value pairs in an `Associative` array, i.e. a dictionary `dims::SectorDict`. +If `IteratorSize(values(I)) isa Union{IsInfinite, SizeUnknown}`, or if `values(I)` has a known length that exceeds `TensorKit._NTUPLE_STORAGE_THRESHOLD`, the different sectors ``a`` and their corresponding degeneracy ``n_a`` are stored as key value pairs in an `Associative` array, i.e. a dictionary `dims::SectorDict`. As the total number of sectors in `values(I)` can be infinite, only sectors ``a`` for which ``n_a`` are stored. Here, `SectorDict` is a constant type alias for a specific dictionary implementation, which currently resorts to `SortedVectorDict` implemented in TensorKit.jl. Hence, the sectors and their corresponding dimensions are stored as two matching lists (`Vector` instances), which are ordered based on the property `isless(a::I, b::I)`. This ensures that the space ``V = ⨁_a ℂ^{n_a} ⊗ R_{a}`` has some unique canonical order in the direct sum decomposition, such that two different but equal instances created independently always match. -If `IteratorSize(values(I)) isa Union{HasLength, HasShape}`, the degeneracy dimensions `n_a` are stored for all sectors `a ∈ values(I)` (also if `n_a == 0`) in a tuple, more specifically a `NTuple{N, Int}` with `N = length(values(I))`. +If `IteratorSize(values(I)) isa Union{HasLength, HasShape}` and `N = length(values(I))` is at most `TensorKit._NTUPLE_STORAGE_THRESHOLD`, the degeneracy dimensions `n_a` are stored for all sectors `a ∈ values(I)` (also if `n_a == 0`) in a tuple, more specifically a `NTuple{N, Int}`. The methods `getindex(values(I), i)` and `findindex(values(I), a)` are used to map between a sector `a ∈ values(I)` and a corresponding index `i ∈ 1:N`. As `N` is a compile time constant, these types can be created in a type stable manner. -Note however that this implies that for large values of `N`, it can be beneficial to define `IteratorSize(values(a)) = SizeUnknown()` to not overly burden the compiler. +For larger `N` this would overly burden the compiler, which is precisely why the dictionary storage takes over above the threshold. +The storage type of a given sector type is reported by `TensorKit.sectorstoragetype(I)`, and the canonical space type is always obtained as `Vect[I]`. ## Constructing instances diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl index 61ae8c99f..cf6f6fd2e 100644 --- a/src/auxiliary/dicts.jl +++ b/src/auxiliary/dicts.jl @@ -178,27 +178,33 @@ function Base.:(==)(d1::SortedVectorDict, d2::SortedVectorDict) end # merge two SortedVectorDicts of `GradedSpace` dimensions, applying `combine` to keys present in -# both; keys present in only one dict are kept as is or dropped according to `_keepunmatched(combine)` -# zero results are dropped since `GradedSpace` never stores an explicit zero dimension -_keepunmatched(::Any) = true -_keepunmatched(::typeof(min)) = false # infimum: a missing sector has dimension zero, so min drops it +# both, and `_unmatched1`/`_unmatched2` to keys present in only the first/second dict, which may +# keep the value, drop the entry by returning `nothing`, or throw. Zero results are always dropped, +# since `GradedSpace` never stores an explicit zero dimension. +_unmatched1(::Any, v) = v +_unmatched2(::Any, v) = v +_unmatched1(::typeof(min), v) = nothing # infimum: a missing sector has dimension zero +_unmatched2(::typeof(min), v) = nothing + +# upper bound on the number of entries the merge can produce +_mergelength(::Any, n1, n2) = n1 + n2 +_mergelength(::typeof(min), n1, n2) = min(n1, n2) function _sortedmerge(combine::F, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {F, K, V <: Integer} - keep = _keepunmatched(combine) k1, v1 = d1.keys, d1.values k2, v2 = d2.keys, d2.values n1, n2 = length(k1), length(k2) - len = keep ? n1 + n2 : min(n1, n2) + len = _mergelength(combine, n1, n2) ks = Vector{K}(undef, len) vs = Vector{V}(undef, len) i, j, n = 1, 1, 0 @inbounds while i <= n1 && j <= n2 a, b = k1[i], k2[j] if isless(a, b) - keep && (n = _mergestore!(ks, vs, n, a, v1[i])) + n = _mergestore!(ks, vs, n, a, _unmatched1(combine, v1[i])) i += 1 elseif isless(b, a) - keep && (n = _mergestore!(ks, vs, n, b, v2[j])) + n = _mergestore!(ks, vs, n, b, _unmatched2(combine, v2[j])) j += 1 else n = _mergestore!(ks, vs, n, a, combine(v1[i], v2[j])) @@ -206,15 +212,13 @@ function _sortedmerge(combine::F, d1::SortedVectorDict{K, V}, d2::SortedVectorDi j += 1 end end - if keep - @inbounds while i <= n1 - n = _mergestore!(ks, vs, n, k1[i], v1[i]) - i += 1 - end - @inbounds while j <= n2 - n = _mergestore!(ks, vs, n, k2[j], v2[j]) - j += 1 - end + @inbounds while i <= n1 + n = _mergestore!(ks, vs, n, k1[i], _unmatched1(combine, v1[i])) + i += 1 + end + @inbounds while j <= n2 + n = _mergestore!(ks, vs, n, k2[j], _unmatched2(combine, v2[j])) + j += 1 end resize!(ks, n) resize!(vs, n) @@ -226,121 +230,11 @@ end @inbounds vs[n + 1] = d return n + !iszero(d) end +@inline _mergestore!(ks, vs, n, k, ::Nothing) = n Base.mergewith(combine, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {K, V <: Integer} = _sortedmerge(combine, d1, d2) -""" - FullVectorDict{K<:Sector,V} <: AbstractDict{K,V} - -Dictionary-like type that reserves one slot for every possible sector `c::K`, indexed -through `findindex(values(K), c)`. Absent entries are stored as `nothing`. -This is the counterpart of `SectorDict` for `NTuple{N,Int}`-backed storage. -Both are needed to build a "sector => value" map over some value type `V`. -""" -struct FullVectorDict{K <: Sector, V} <: AbstractDict{K, V} - slots::Vector{Union{Nothing, V}} - function FullVectorDict{K, V}(slots::Vector{Union{Nothing, V}}) where {K <: Sector, V} - @assert length(slots) == length(values(K)) - return new{K, V}(slots) - end -end -FullVectorDict{K, V}() where {K <: Sector, V} = - FullVectorDict{K, V}(Vector{Union{Nothing, V}}(nothing, length(values(K)))) -function FullVectorDict{K, V}(kv) where {K <: Sector, V} - d = FullVectorDict{K, V}() - for (k, v) in kv - d[k] = v - end - return d -end -FullVectorDict{K, V}(kv::Pair{K, V}...) where {K <: Sector, V} = FullVectorDict{K, V}(kv) - -Base.length(d::FullVectorDict) = count(!isnothing, d.slots) - -Base.copy(d::FullVectorDict{K, V}) where {K, V} = FullVectorDict{K, V}(copy(d.slots)) -Base.empty(::FullVectorDict{K}, ::Type{K}, ::Type{V}) where {K <: Sector, V} = FullVectorDict{K, V}() -Base.empty!(d::FullVectorDict) = (fill!(d.slots, nothing); return d) - -function Base.delete!(d::FullVectorDict{K}, k) where {K} - key = convert(K, k) - isequal(k, key) && (d.slots[findindex(values(K), key)] = nothing) - return d -end - -function Base.haskey(d::FullVectorDict{K}, k) where {K} - key = convert(K, k) - return isequal(k, key) && !isnothing(d.slots[findindex(values(K), key)]) -end -function Base.getindex(d::FullVectorDict{K}, k) where {K} - key = convert(K, k) - isequal(k, key) || throw(KeyError(k)) - v = d.slots[findindex(values(K), key)] - return isnothing(v) ? throw(KeyError(key)) : v -end -function Base.setindex!(d::FullVectorDict{K}, v, k) where {K} - key = convert(K, k) - isequal(k, key) || throw(ArgumentError("$k is not a valid key for type $K")) - d.slots[findindex(values(K), key)] = v - return d -end - -function Base.get(d::FullVectorDict{K}, k, default) where {K} - key = convert(K, k) - isequal(k, key) || return default - v = d.slots[findindex(values(K), key)] - return isnothing(v) ? default : v -end -function Base.get(f::Union{Function, Type}, d::FullVectorDict{K}, k) where {K} - key = convert(K, k) - isequal(k, key) || return f() - v = d.slots[findindex(values(K), key)] - return isnothing(v) ? f() : v -end - -function Base.iterate(d::FullVectorDict{K}, i = 1) where {K} - vals = values(K) - n = length(d.slots) - @inbounds while i <= n - v = d.slots[i] - isnothing(v) || return (vals[i] => v), i + 1 - i += 1 - end - return nothing -end - -# the dense sector map is backed by a `FullVectorDict` or a `SortedVectorDict`, mirroring the -# `NTuple`/`SortedVectorDict` storage choice that `sectorstoragetype` makes for `GradedSpace` -_densemaptype(::Type{<:Tuple}, ::Type{I}, ::Type{V}) where {I <: Sector, V} = FullVectorDict{I, V} -_densemaptype(::Type{<:SortedVectorDict}, ::Type{I}, ::Type{V}) where {I <: Sector, V} = - SortedVectorDict{I, V} - -""" - sectormaptype(::Type{I}, ::Type{V}) where {I <: Sector, V} -> Type - -The dense `AbstractDict{I, V}` type used to map sectors of type `I` onto values of type `V`, -chosen to match the storage type of `GradedSpace{I}`, see [`sectorstoragetype`](@ref). -""" -sectormaptype(::Type{I}, ::Type{V}) where {I <: Sector, V} = - _densemaptype(sectorstoragetype(I), I, V) - -""" - sectormap(f, [V::Type, ] pairsiter) - -Map `f(c, v)` over an iterator `pairsiter` of `c => v` pairs, collecting the results into the -dense sector map returned by [`sectormaptype`](@ref). The value type `V` of the result is -inferred from `f` unless it is given explicitly. -""" -sectormap(f, pairsiter) = _sectormap(f, pairsiter, eltype(pairsiter)) -sectormap(f, ::Type{V}, pairsiter) where {V} = _sectormap(f, V, pairsiter, eltype(pairsiter)) - -function _sectormap(f, pairsiter, ::Type{<:Pair{I, W}}) where {I <: Sector, W} - return _sectormap(f, Base.promote_op(f, I, W), pairsiter, Pair{I, W}) -end -function _sectormap(f, ::Type{V}, pairsiter, ::Type{<:Pair{I}}) where {V, I <: Sector} - return sectormaptype(I, V)(c => f(c, v) for (c, v) in pairsiter) -end - """ Hashed(value, hashfunction = Base.hash, isequal = Base.isequal) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 9b214b776..0c41648ce 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -8,7 +8,7 @@ export copy_oftype, factorisation_scalartype, one!, truncspace using ..TensorKit using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, findindex, blocktype, foreachblock, one!, - similar_diagonal, similarstoragetype, sectormap + similar_diagonal, similarstoragetype using ..TensorKit: GLOBAL_TIMER using TimerOutputs: @timeit_debug diff --git a/src/factorizations/pullbacks.jl b/src/factorizations/pullbacks.jl index b63bc9a94..a74acf8df 100644 --- a/src/factorizations/pullbacks.jl +++ b/src/factorizations/pullbacks.jl @@ -24,7 +24,7 @@ for pullback! in (:qr_null_pullback!, :lq_null_pullback!) return Δt end end -_notrunc_ind(t) = sectormap(Returns(Colon()), blocks(t)) +_notrunc_ind(t) = SectorDict(c => Colon() for c in blocksectors(t)) for pullback! in (:svd_pullback!, :eig_pullback!, :eigh_pullback!) @eval function MAK.$pullback!( diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 71bb54b68..ed595b3c3 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -36,8 +36,6 @@ _blocklength(ax::Base.OneTo, ind::AbstractVector{Bool}) = count(ind) function truncate_space(V::ElementarySpace, inds) @assert !isdual(V) - I = sectortype(V) - @assert I == Trivial return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in pairs(inds)) end function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N} @@ -68,27 +66,27 @@ end function truncate_domain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds) for (c, b) in blocks(tdst) - I = get(inds, c, nothing) - @assert !isnothing(I) # kept for safety, but should be guaranteed by pairs(inds) covering every block sector + ind = get(inds, c, nothing) + @assert !isnothing(ind) b′ = block(tsrc, c) - b .= view(b′, :, I) + b .= view(b′, :, ind) end return tdst end function truncate_codomain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds) for (c, b) in blocks(tdst) - I = get(inds, c, nothing) - @assert !isnothing(I) # kept for safety, but should be guaranteed by pairs(inds) covering every block sector + ind = get(inds, c, nothing) + @assert !isnothing(ind) b′ = block(tsrc, c) - b .= view(b′, I, :) + b .= view(b′, ind, :) end return tdst end function truncate_diagonal!(Ddst::DiagonalTensorMap, Dsrc::DiagonalTensorMap, inds) for (c, b) in blocks(Ddst) - I = get(inds, c, nothing) - @assert !isnothing(I) # kept for safety, but should be guaranteed by pairs(inds) covering every block sector - diagview(b) .= view(diagview(block(Dsrc, c)), I) + ind = get(inds, c, nothing) + @assert !isnothing(ind) + diagview(b) .= view(diagview(block(Dsrc, c)), ind) end return Ddst end @@ -97,17 +95,17 @@ function MAK.truncate( ::typeof(svd_trunc!), (U, S, Vᴴ)::NTuple{3, AbstractTensorMap}, strategy::TruncationStrategy ) - ind = MAK.findtruncated_svd(diagview(S), strategy) - V_truncated = truncate_space(space(S, 1), ind) + inds = MAK.findtruncated_svd(diagview(S), strategy) + V_truncated = truncate_space(space(S, 1), inds) Ũ = similar(U, codomain(U) ← V_truncated) - truncate_domain!(Ũ, U, ind) + truncate_domain!(Ũ, U, inds) S̃ = similar_diagonal(S, V_truncated) - truncate_diagonal!(S̃, S, ind) + truncate_diagonal!(S̃, S, inds) Ṽᴴ = similar(Vᴴ, V_truncated ← domain(Vᴴ)) - truncate_codomain!(Ṽᴴ, Vᴴ, ind) + truncate_codomain!(Ṽᴴ, Vᴴ, inds) - return (Ũ, S̃, Ṽᴴ), ind + return (Ũ, S̃, Ṽᴴ), inds end function MAK.truncate( @@ -117,11 +115,11 @@ function MAK.truncate( for (c, b) in blocks(S) copyto!(extended_S[c], diagview(b)) # copyto! since `b` might be shorter end - ind = MAK.findtruncated(extended_S, strategy) - V_truncated = truncate_space(space(S, 1), ind) + inds = MAK.findtruncated(extended_S, strategy) + V_truncated = truncate_space(space(S, 1), inds) Ũ = similar(U, codomain(U) ← V_truncated) - truncate_domain!(Ũ, U, ind) - return Ũ, ind + truncate_domain!(Ũ, U, inds) + return Ũ, inds end function MAK.truncate( ::typeof(right_null!), (S, Vᴴ)::NTuple{2, AbstractTensorMap}, strategy::TruncationStrategy @@ -130,11 +128,11 @@ function MAK.truncate( for (c, b) in blocks(S) copyto!(extended_S[c], diagview(b)) # copyto! since `b` might be shorter end - ind = MAK.findtruncated(extended_S, strategy) - V_truncated = truncate_space(dual(space(S, 2)), ind) + inds = MAK.findtruncated(extended_S, strategy) + V_truncated = truncate_space(dual(space(S, 2)), inds) Ṽᴴ = similar(Vᴴ, V_truncated ← domain(Vᴴ)) - truncate_codomain!(Ṽᴴ, Vᴴ, ind) - return Ṽᴴ, ind + truncate_codomain!(Ṽᴴ, Vᴴ, inds) + return Ṽᴴ, inds end # special case `NoTruncation` for null: should keep exact zeros due to rectangularity @@ -142,24 +140,20 @@ end function MAK.truncate( ::typeof(left_null!), (U, S)::NTuple{2, AbstractTensorMap}, strategy::NoTruncation ) - ind = sectormap(blocks(S)) do _, b - (size(b, 2) + 1):size(b, 1) - end - V_truncated = truncate_space(space(S, 1), ind) + inds = SectorDict(c => (size(b, 2) + 1):size(b, 1) for (c, b) in blocks(S)) + V_truncated = truncate_space(space(S, 1), inds) Ũ = similar(U, codomain(U) ← V_truncated) - truncate_domain!(Ũ, U, ind) - return Ũ, ind + truncate_domain!(Ũ, U, inds) + return Ũ, inds end function MAK.truncate( ::typeof(right_null!), (S, Vᴴ)::NTuple{2, AbstractTensorMap}, strategy::NoTruncation ) - ind = sectormap(blocks(S)) do _, b - (size(b, 1) + 1):size(b, 2) - end - V_truncated = truncate_space(dual(space(S, 2)), ind) + inds = SectorDict(c => (size(b, 1) + 1):size(b, 2) for (c, b) in blocks(S)) + V_truncated = truncate_space(dual(space(S, 2)), inds) Ṽᴴ = similar(Vᴴ, V_truncated ← domain(Vᴴ)) - truncate_codomain!(Ṽᴴ, Vᴴ, ind) - return Ṽᴴ, ind + truncate_codomain!(Ṽᴴ, Vᴴ, inds) + return Ṽᴴ, inds end for f! in (:eig_trunc!, :eigh_trunc!) @@ -168,16 +162,16 @@ for f! in (:eig_trunc!, :eigh_trunc!) (D, V)::Tuple{DiagonalTensorMap, AbstractTensorMap}, strategy::TruncationStrategy ) - ind = MAK.findtruncated(diagview(D), strategy) - V_truncated = truncate_space(space(D, 1), ind) + inds = MAK.findtruncated(diagview(D), strategy) + V_truncated = truncate_space(space(D, 1), inds) D̃ = similar_diagonal(D, V_truncated) - truncate_diagonal!(D̃, D, ind) + truncate_diagonal!(D̃, D, inds) Ṽ = similar(V, codomain(V) ← V_truncated) - truncate_domain!(Ṽ, V, ind) + truncate_domain!(Ṽ, V, inds) - return (D̃, Ṽ), ind + return (D̃, Ṽ), inds end end @@ -192,7 +186,7 @@ function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationStrateg end function MAK.findtruncated(values::SectorVector, ::NoTruncation) - return sectormap(Returns(Colon()), pairs(values)) + return SectorDict(c => Colon() for c in keys(values)) end # Need to select the first k values here after sorting across blocks, weighted by quantum dimension @@ -238,24 +232,18 @@ MAK.findtruncated_svd(values::SectorVector, strategy::TruncationByOrder) = MAK.findtruncated(values, strategy) function MAK.findtruncated(values::SectorVector, strategy::TruncationByFilter) - return sectormap(pairs(values)) do _, v - findall(strategy.filter, v) - end + return SectorDict(c => findall(strategy.filter, d) for (c, d) in pairs(values)) end function MAK.findtruncated(values::SectorVector, strategy::TruncationByValue) atol = rtol_to_atol(values, strategy.p, strategy.atol, strategy.rtol) strategy′ = trunctol(; atol, strategy.by, strategy.keep_below) - return sectormap(pairs(values)) do _, v - MAK.findtruncated(v, strategy′) - end + return SectorDict(c => MAK.findtruncated(d, strategy′) for (c, d) in pairs(values)) end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationByValue) atol = rtol_to_atol(values, strategy.p, strategy.atol, strategy.rtol) strategy′ = trunctol(; atol, strategy.by, strategy.keep_below) - return sectormap(pairs(values)) do _, v - MAK.findtruncated_svd(v, strategy′) - end + return SectorDict(c => MAK.findtruncated_svd(d, strategy′) for (c, d) in pairs(values)) end # Need to select the first k values here after sorting by error across blocks, @@ -297,20 +285,14 @@ MAK.findtruncated_svd(values::SectorVector, strategy::TruncationByError) = MAK.findtruncated(values, strategy) function MAK.findtruncated(values::SectorVector, strategy::TruncationSpace) - I = sectortype(values) - I == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values")) + sectortype(values) == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values")) blockstrategy(c) = truncrank(dim(strategy.space, c); strategy.by, strategy.rev) - return sectormap(pairs(values)) do c, v - MAK.findtruncated(v, blockstrategy(c)) - end + return SectorDict(c => MAK.findtruncated(d, blockstrategy(c)) for (c, d) in pairs(values)) end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationSpace) - I = sectortype(values) - I == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values")) + sectortype(values) == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values")) blockstrategy(c) = truncrank(dim(strategy.space, c); strategy.by, strategy.rev) - return sectormap(pairs(values)) do c, v - MAK.findtruncated_svd(v, blockstrategy(c)) - end + return SectorDict(c => MAK.findtruncated_svd(d, blockstrategy(c)) for (c, d) in pairs(values)) end # The implementations below assume that the `SectorDict` always contains an entry for every block sector @@ -319,31 +301,39 @@ end function MAK.findtruncated(values::SectorVector, strategy::TruncationIntersection) inds = map(Base.Fix1(MAK.findtruncated, values), strategy.components) - @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" - return sectormap(Any, pairs(first(inds))) do c, _ - mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_intersect, inds) + @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now" + sectors = keys(first(inds)) + vals = map(keys(first(inds))) do c + mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_intersect, inds) end + return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals) end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationIntersection) inds = map(Base.Fix1(MAK.findtruncated_svd, values), strategy.components) - @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" - return sectormap(Any, pairs(first(inds))) do c, _ - mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_intersect, inds) + @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now" + sectors = keys(first(inds)) + vals = map(keys(first(inds))) do c + mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_intersect, inds) end + return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals) end function MAK.findtruncated(values::SectorVector, strategy::TruncationUnion) inds = map(Base.Fix1(MAK.findtruncated, values), strategy.components) - @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" - return sectormap(Any, pairs(first(inds))) do c, _ - mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_union, inds) + @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now" + sectors = keys(first(inds)) + vals = map(keys(first(inds))) do c + mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_union, inds) end + return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals) end function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationUnion) inds = map(Base.Fix1(MAK.findtruncated_svd, values), strategy.components) - @assert TensorKit._allequal(v -> collect(keys(v)), inds) "missing blocks are not supported right now" - return sectormap(Any, pairs(first(inds))) do c, _ - mapreduce(v -> get(v, c, nothing), MatrixAlgebraKit._ind_union, inds) + @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now" + sectors = keys(first(inds)) + vals = map(keys(first(inds))) do c + mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_union, inds) end + return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals) end # Truncation error diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 61b407b17..cd197c68d 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -1,3 +1,6 @@ +# budget on how many sectors are worth specializing the compiler on +const _NTUPLE_STORAGE_THRESHOLD = 8 + """ struct GradedSpace{I<:Sector, D} <: ElementarySpace GradedSpace{I,D}(dims; dual::Bool = false) where {I<:Sector, D} @@ -12,13 +15,13 @@ and pivotal (pre-, multi-) fusion category. Here `dims` represents the degeneracy or multiplicity of every sector. The data structure `D` of `dims` will depend on the result `Base.IteratorSize(values(I))`. -If the result is of type `HasLength` or `HasShape`, `dims` will be stored in a -`NTuple{N,Int}` with `N = length(values(I))`. This requires that a sector `s::I` can be -transformed into an index via `s == getindex(values(I), i)` and -`i == findindex(values(I), s)`. If `Base.IteratorElsize(values(I))` results `IsInfinite()` -or `SizeUnknown()`, a `SectorDict{I,Int}` is used to store the non-zero degeneracy -dimensions with the corresponding sector as key. The parameter `D` is hidden from the user -and should typically be of no concern. +If the result is of type `HasLength` or `HasShape` and `N = length(values(I))` does not exceed +`$_NTUPLE_STORAGE_THRESHOLD`, `dims` will be stored in a `NTuple{N,Int}`. This requires that a +sector `s::I` can be transformed into an index via `s == getindex(values(I), i)` and +`i == findindex(values(I), s)`. Otherwise, i.e. for more sectors than that or if +`Base.IteratorSize(values(I))` results `IsInfinite()` or `SizeUnknown()`, a `SectorDict{I,Int}` +is used to store the non-zero degeneracy dimensions with the corresponding sector as key. +The parameter `D` is hidden from the user and should typically be of no concern. The concrete type `GradedSpace{I,D}` with correct `D` can be obtained as `Vect[I]`, or if `I == Irrep[G]` for some `G<:Group`, as `Rep[G]`. @@ -172,36 +175,21 @@ function ⊖(V::GradedSpace{I, <:Tuple}, W::GradedSpace{I, <:Tuple}) where {I <: newdims = map((dV, dW) -> dV < dW ? _throw_not_subspace(V, W) : dV - dW, V.dims, W.dims) return typeof(V)(newdims, dualV) end + +# combiner for `⊖`, carrying the spaces to report which subspace condition was violated +struct SubtractDims{S} + V::S + W::S +end +(f::SubtractDims)(dV, dW) = dV < dW ? _throw_not_subspace(f.V, f.W) : dV - dW +_unmatched1(::SubtractDims, d) = d # sector only in `V`: keep +_unmatched2(f::SubtractDims, d) = _throw_not_subspace(f.V, f.W) # sector only in `W`: not a subspace +_mergelength(::SubtractDims, n1, n2) = n1 + function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) where {I <: Sector} dualV = isdual(V) dualV == isdual(W) || _throw_not_subspace(V, W) - k1, v1 = keys(V.dims), values(V.dims) - k2, v2 = keys(W.dims), values(W.dims) - n1, n2 = length(k1), length(k2) - ks, vs = Vector{I}(undef, n1), Vector{Int}(undef, n1) - i, j, n = 1, 1, 0 - @inbounds while i <= n1 && j <= n2 - a, b = k1[i], k2[j] - if isless(a, b) - n = _mergestore!(ks, vs, n, a, v1[i]) - i += 1 - elseif isless(b, a) - _throw_not_subspace(V, W) # sector of `W` absent from `V` - else - v1[i] < v2[j] && _throw_not_subspace(V, W) - n = _mergestore!(ks, vs, n, a, v1[i] - v2[j]) - i += 1 - j += 1 - end - end - j <= n2 && _throw_not_subspace(V, W) # leftover sectors of `W` absent from `V` - @inbounds while i <= n1 - n = _mergestore!(ks, vs, n, k1[i], v1[i]) - i += 1 - end - resize!(ks, n) - resize!(vs, n) - return typeof(V)(SectorDict{I, Int}(ks, vs), dualV) + return typeof(V)(mergewith(SubtractDims(V, W), V.dims, W.dims), dualV) end function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} @@ -375,20 +363,17 @@ Base.getindex(::SpaceTable) = ComplexSpace Base.getindex(::SpaceTable, ::Type{Trivial}) = ComplexSpace Base.getindex(::SpaceTable, I::Type{<:Sector}) = GradedSpace{I, sectorstoragetype(I)} -# based on Julia tuple unrolling range -const _ntuple_storage_threshold = 32 - """ sectorstoragetype(I::Type{<:Sector}) -> Type The storage type `D` used for the `dims` field of `GradedSpace{I, D}`. This is `NTuple{N,Int}` with `N = length(values(I))` if `I` has a finite, known length -of at most `$_ntuple_storage_threshold`, or `SectorDict{I,Int}` otherwise. +of at most `$_NTUPLE_STORAGE_THRESHOLD`, or `SectorDict{I,Int}` otherwise. """ Base.@assume_effects :foldable function sectorstoragetype(::Type{I}) where {I <: Sector} if Base.IteratorSize(values(I)) isa Union{HasLength, HasShape} N = length(values(I)) - N <= _ntuple_storage_threshold && return NTuple{N, Int} + N <= _NTUPLE_STORAGE_THRESHOLD && return NTuple{N, Int} end return SectorDict{I, Int} end @@ -426,17 +411,33 @@ function type_repr(::Type{<:GradedSpace{ProductSector{T}}}) where end # Specific constructors for Z_N +""" + const ZNSpace{N} + +Type alias for the tuple-backed `GradedSpace{ZNIrrep{N}, NTuple{N,Int}}`. + +!!! warning "Deprecated" + A type alias cannot compute the storage type from `N`, so this only coincides with + `Vect[ZNIrrep{N}]` while `N <= $_NTUPLE_STORAGE_THRESHOLD`. Use `Vect[ZNIrrep{N}]` instead. +""" const ZNSpace{N} = GradedSpace{ZNIrrep{N}, NTuple{N, Int}} -ZNSpace{N}(dims::NTuple{N, Int}; dual::Bool = false) where {N} = ZNSpace{N}(dims, dual) -ZNSpace{N}(dims::Vararg{Int, N}; dual::Bool = false) where {N} = ZNSpace{N}(dims, dual) -ZNSpace(dims::NTuple{N, Int}; dual::Bool = false) where {N} = ZNSpace{N}(dims, dual) -ZNSpace(dims::Vararg{Int, N}; dual::Bool = false) where {N} = ZNSpace{N}(dims, dual) +@noinline function _throw_znspace_storage(N) + msg = lazy"`ZNSpace{$N}` is not the canonical space type for `ZNIrrep{$N}`, which stores more than $(_NTUPLE_STORAGE_THRESHOLD) sectors in a `SectorDict`; use `Vect[ZNIrrep{$N}]` instead" + return throw(ArgumentError(msg)) +end +function ZNSpace{N}(dims::NTuple{N, Int}; dual::Bool = false) where {N} + N <= _NTUPLE_STORAGE_THRESHOLD || _throw_znspace_storage(N) + return ZNSpace{N}(dims, dual) +end +ZNSpace{N}(dims::Vararg{Int, N}; dual::Bool = false) where {N} = ZNSpace{N}(dims; dual) +ZNSpace(dims::NTuple{N, Int}; dual::Bool = false) where {N} = ZNSpace{N}(dims; dual) +ZNSpace(dims::Vararg{Int, N}; dual::Bool = false) where {N} = ZNSpace{N}(dims; dual) # TODO: Do we still need all of those # ASCII type aliases -const Z2Space = ZNSpace{2} -const Z3Space = ZNSpace{3} -const Z4Space = ZNSpace{4} +const Z2Space = Vect[ZNIrrep{2}] +const Z3Space = Vect[ZNIrrep{3}] +const Z4Space = Vect[ZNIrrep{4}] const U1Space = Rep[U₁] const CU1Space = Rep[CU₁] const SU2Space = Rep[SU₂] diff --git a/src/tensors/sectorvector.jl b/src/tensors/sectorvector.jl index f586d4aef..5744c23f5 100644 --- a/src/tensors/sectorvector.jl +++ b/src/tensors/sectorvector.jl @@ -53,7 +53,11 @@ Base.valtype(::Type{SectorVector{T, I, A}}) where {T, I, A} = SubArray{T, 1, A, Base.keys(v::SectorVector) = keys(v.structure) Base.values(v::SectorVector) = (v[c] for c in keys(v)) -Base.pairs(v::SectorVector) = SectorDict(c => v[c] for c in keys(v)) +function Base.pairs(v::SectorVector) + # `structure` is already sorted, so the dict can be built without repeated insertion + vals = valtype(v)[view(parent(v), r) for r in values(v.structure)] + return SectorDict{keytype(v), valtype(v)}(copy(v.structure.keys), vals) +end Base.get(v::SectorVector{<:Any, I}, key::I, default) where {I} = haskey(v, key) ? v[key] : default Base.haskey(v::SectorVector{<:Any, I}, key::I) where {I} = key in keys(v) From 192cf486144bdea71208724480ed8dddc1a12e63 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 15 Sep 2026 00:48:18 -0400 Subject: [PATCH 32/36] fix `pairs(::SectorVector)` for non-`SubArray` views `valtype(::SectorVector)` claims a `SubArray`, but `view` of a GPU array is itself a `CuArray`/`ROCArray`, so pinning the element type to `valtype` broke every GPU factorization. Take the element type from the views instead. Co-Authored-By: Claude Opus 5 (1M context) --- src/tensors/sectorvector.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/tensors/sectorvector.jl b/src/tensors/sectorvector.jl index 5744c23f5..11be45826 100644 --- a/src/tensors/sectorvector.jl +++ b/src/tensors/sectorvector.jl @@ -55,8 +55,9 @@ Base.keys(v::SectorVector) = keys(v.structure) Base.values(v::SectorVector) = (v[c] for c in keys(v)) function Base.pairs(v::SectorVector) # `structure` is already sorted, so the dict can be built without repeated insertion - vals = valtype(v)[view(parent(v), r) for r in values(v.structure)] - return SectorDict{keytype(v), valtype(v)}(copy(v.structure.keys), vals) + # note: `view` need not produce `valtype(v)`, e.g. for GPU arrays + vals = map(Base.Fix1(view, parent(v)), values(v.structure)) + return SectorDict{keytype(v), eltype(vals)}(copy(v.structure.keys), vals) end Base.get(v::SectorVector{<:Any, I}, key::I, default) where {I} = haskey(v, key) ? v[key] : default From d1802c547322a397fea1034b76ade8952b22c316 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 15 Sep 2026 11:14:07 -0400 Subject: [PATCH 33/36] review: storage-variant aliases, shared merge handlers, lazy `pairs` MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - add `TupleGradedSpace{I,N}` / `DictGradedSpace{I}` aliases for the two variants `sectorstoragetype` selects between, and use them for dispatch - `SubtractDims` doubles as the unmatched handler for `⊖`, and is used for the tuple variant as well, so both paths share one callable - `_sortedmerge` takes the unmatched handlers as arguments; `mergewith` selects them inline - `pairs(::SectorVector)` is lazy, like `blocks(::AbstractTensorMap)`: every consumer only iterates, and lookups go through the vector itself - document `sectorstoragetype`, whose docstring interpolates the threshold Co-Authored-By: Claude Opus 5 (1M context) --- docs/src/lib/spaces.md | 13 ++++-- docs/src/man/gradedspaces.md | 2 +- src/auxiliary/dicts.jl | 50 +++++++++++++--------- src/factorizations/factorizations.jl | 3 +- src/factorizations/truncation.jl | 4 +- src/spaces/gradedspace.jl | 64 +++++++++++++++------------- src/tensors/sectorvector.jl | 9 ++-- 7 files changed, 80 insertions(+), 65 deletions(-) diff --git a/docs/src/lib/spaces.md b/docs/src/lib/spaces.md index fbeaacb6d..228c98726 100644 --- a/docs/src/lib/spaces.md +++ b/docs/src/lib/spaces.md @@ -37,13 +37,18 @@ Vect Rep ``` +The storage that `Vect[I]` selects for the degeneracy dimensions depends on how many sectors `I` has: + +```@docs +TensorKit.sectorstoragetype +``` + In this respect, there are also a number of type aliases for the `GradedSpace` types associated with the most common sectors, namely ```julia -const ZNSpace{N} = Vect[ZNIrrep{N}] -const Z2Space = ZNSpace{2} -const Z3Space = ZNSpace{3} -const Z4Space = ZNSpace{4} +const Z2Space = Vect[ZNIrrep{2}] +const Z3Space = Vect[ZNIrrep{3}] +const Z4Space = Vect[ZNIrrep{4}] const U1Space = Rep[U₁] const CU1Space = Rep[CU₁] const SU2Space = Rep[SU₂] diff --git a/docs/src/man/gradedspaces.md b/docs/src/man/gradedspaces.md index 529586cf3..2d0f986c2 100644 --- a/docs/src/man/gradedspaces.md +++ b/docs/src/man/gradedspaces.md @@ -35,7 +35,7 @@ If `IteratorSize(values(I)) isa Union{HasLength, HasShape}` and `N = length(valu The methods `getindex(values(I), i)` and `findindex(values(I), a)` are used to map between a sector `a ∈ values(I)` and a corresponding index `i ∈ 1:N`. As `N` is a compile time constant, these types can be created in a type stable manner. For larger `N` this would overly burden the compiler, which is precisely why the dictionary storage takes over above the threshold. -The storage type of a given sector type is reported by `TensorKit.sectorstoragetype(I)`, and the canonical space type is always obtained as `Vect[I]`. +The exact threshold is documented with [`TensorKit.sectorstoragetype`](@ref), which reports the storage type of a given sector type; the canonical space type is always obtained as `Vect[I]`. ## Constructing instances diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl index cf6f6fd2e..642373e19 100644 --- a/src/auxiliary/dicts.jl +++ b/src/auxiliary/dicts.jl @@ -177,34 +177,28 @@ function Base.:(==)(d1::SortedVectorDict, d2::SortedVectorDict) return true end -# merge two SortedVectorDicts of `GradedSpace` dimensions, applying `combine` to keys present in -# both, and `_unmatched1`/`_unmatched2` to keys present in only the first/second dict, which may -# keep the value, drop the entry by returning `nothing`, or throw. Zero results are always dropped, -# since `GradedSpace` never stores an explicit zero dimension. -_unmatched1(::Any, v) = v -_unmatched2(::Any, v) = v -_unmatched1(::typeof(min), v) = nothing # infimum: a missing sector has dimension zero -_unmatched2(::typeof(min), v) = nothing - -# upper bound on the number of entries the merge can produce -_mergelength(::Any, n1, n2) = n1 + n2 -_mergelength(::typeof(min), n1, n2) = min(n1, n2) - -function _sortedmerge(combine::F, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {F, K, V <: Integer} +# merge two SortedVectorDicts of `GradedSpace` dimensions: `combine(v1, v2)` is applied to keys +# present in both, `unmatched1(v)`/`unmatched2(v)` to keys present in only the first/second dict, +# which may keep the value, drop the entry by returning `nothing`, or throw. Zero results are +# always dropped, since `GradedSpace` never stores an explicit zero dimension. +function _sortedmerge( + combine::F, unmatched1::F1, unmatched2::F2, + d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V} + ) where {F, F1, F2, K, V <: Integer} k1, v1 = d1.keys, d1.values k2, v2 = d2.keys, d2.values n1, n2 = length(k1), length(k2) - len = _mergelength(combine, n1, n2) + len = _mergelength(unmatched1, unmatched2, n1, n2) ks = Vector{K}(undef, len) vs = Vector{V}(undef, len) i, j, n = 1, 1, 0 @inbounds while i <= n1 && j <= n2 a, b = k1[i], k2[j] if isless(a, b) - n = _mergestore!(ks, vs, n, a, _unmatched1(combine, v1[i])) + n = _mergestore!(ks, vs, n, a, unmatched1(v1[i])) i += 1 elseif isless(b, a) - n = _mergestore!(ks, vs, n, b, _unmatched2(combine, v2[j])) + n = _mergestore!(ks, vs, n, b, unmatched2(v2[j])) j += 1 else n = _mergestore!(ks, vs, n, a, combine(v1[i], v2[j])) @@ -213,11 +207,11 @@ function _sortedmerge(combine::F, d1::SortedVectorDict{K, V}, d2::SortedVectorDi end end @inbounds while i <= n1 - n = _mergestore!(ks, vs, n, k1[i], _unmatched1(combine, v1[i])) + n = _mergestore!(ks, vs, n, k1[i], unmatched1(v1[i])) i += 1 end @inbounds while j <= n2 - n = _mergestore!(ks, vs, n, k2[j], _unmatched2(combine, v2[j])) + n = _mergestore!(ks, vs, n, k2[j], unmatched2(v2[j])) j += 1 end resize!(ks, n) @@ -232,8 +226,22 @@ end end @inline _mergestore!(ks, vs, n, k, ::Nothing) = n -Base.mergewith(combine, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {K, V <: Integer} = - _sortedmerge(combine, d1, d2) +# upper bound on the number of entries the merge can produce: a dropping handler contributes none +function _mergelength(unmatched1, unmatched2, n1, n2) + drop = Returns(nothing) + return if unmatched1 == drop + unmatched2 == drop ? min(n1, n2) : n2 + else + unmatched2 == drop ? n1 : n1 + n2 + end +end + +function Base.mergewith(combine, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {K, V <: Integer} + # keys occurring in only one of the dicts are kept, except for `min`, where a missing sector + # has dimension zero and thus drops out of the result + unmatched = combine == min ? Returns(nothing) : identity + return _sortedmerge(combine, unmatched, unmatched, d1, d2) +end """ Hashed(value, hashfunction = Base.hash, isequal = Base.isequal) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 0c41648ce..12b8b4437 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -6,7 +6,8 @@ module Factorizations export copy_oftype, factorisation_scalartype, one!, truncspace using ..TensorKit -using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, findindex, +using ..TensorKit: AdjointTensorMap, DictGradedSpace, SectorDict, SectorVector, + TupleGradedSpace, findindex, blocktype, foreachblock, one!, similar_diagonal, similarstoragetype using ..TensorKit: GLOBAL_TIMER diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index ed595b3c3..3e6ef8ee8 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -38,7 +38,7 @@ function truncate_space(V::ElementarySpace, inds) @assert !isdual(V) return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in pairs(inds)) end -function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N} +function truncate_space(V::TupleGradedSpace{I, N}, inds) where {I <: Sector, N} @assert !isdual(V) vals = values(I) newdims = MutableNTuple(ntuple(_ -> 0, StaticLength(N))) @@ -49,7 +49,7 @@ function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sec end return typeof(V)(Tuple(newdims), false) end -function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector} +function truncate_space(V::DictGradedSpace{I}, inds) where {I <: Sector} @assert !isdual(V) ks, vs = Vector{I}(), Vector{Int}() # accumulate and sort once at the end for (c, ind) in pairs(inds) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index cd197c68d..d46137043 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -32,6 +32,10 @@ struct GradedSpace{I <: Sector, D} <: ElementarySpace end sectortype(::Type{<:GradedSpace{I}}) where {I <: Sector} = I +# the two storage variants that `sectorstoragetype` selects between, used for dispatch +const TupleGradedSpace{I, N} = GradedSpace{I, NTuple{N, Int}} +const DictGradedSpace{I} = GradedSpace{I, SectorDict{I, Int}} + # elementary spaces are homogeneously colored: all sectors share a left and a right unit. function _check_unit_homogeneity(::Type{I}, sectors) where {I <: Sector} (UnitStyle(I) isa SimpleUnit || isempty(sectors)) && return nothing @@ -109,10 +113,10 @@ function dim(V::GradedSpace{I}) where {I <: Sector} init = zero(dimscalartype(I)) return sum(((c, d),) -> dim(c) * d, blockdims(V); init) end -function dim(V::GradedSpace{I, <:AbstractDict}, c::I) where {I <: Sector} +function dim(V::DictGradedSpace{I}, c::I) where {I <: Sector} return get(V.dims, isdual(V) ? dual(c) : c, 0) end -function dim(V::GradedSpace{I, <:Tuple}, c::I) where {I <: Sector} +function dim(V::TupleGradedSpace{I}, c::I) where {I <: Sector} return V.dims[findindex(values(I), isdual(V) ? dual(c) : c)] end Base.axes(V::GradedSpace) = Base.OneTo(dim(V)) @@ -141,7 +145,7 @@ Base.@assume_effects :foldable function _dualpermutation(::Type{I}, ::Val{N}) wh vals = values(I) return ntuple(n -> findindex(vals, dual(vals[n])), Val(N)) end -function flip(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} +function flip(V::TupleGradedSpace{I, N}) where {I <: Sector, N} # `flip` maps `c => d` to `dual(c) => d` and negates `isdual`, which for tuple storage is # just a fixed permutation of the dims, so the generic constructor can be skipped newdims = TupleTools.getindices(V.dims, _dualpermutation(I, Val(N))) @@ -155,12 +159,12 @@ function unitspace(S::Type{<:GradedSpace{I}}) where {I <: Sector} end zerospace(S::Type{<:GradedSpace}) = S() -function ⊕(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} +function ⊕(V₁::DictGradedSpace{I}, V₂::DictGradedSpace{I}) where {I <: Sector} dual1 = isdual(V₁) dual1 == isdual(V₂) || throw(SpaceMismatch("Direct sum of a vector space and a dual space does not exist")) return typeof(V₁)(mergewith(+, V₁.dims, V₂.dims), dual1) end -function ⊕(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} +function ⊕(V₁::TupleGradedSpace{I}, V₂::TupleGradedSpace{I}) where {I <: Sector} dual1 = isdual(V₁) dual1 == isdual(V₂) || throw(SpaceMismatch("Direct sum of a vector space and a dual space does not exist")) @@ -168,31 +172,31 @@ function ⊕(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where return typeof(V₁)(newdims, dual1) end @noinline _throw_not_subspace(V, W) = throw(SpaceMismatch(lazy"$(W) is not a subspace of $(V)")) -function ⊖(V::GradedSpace{I, <:Tuple}, W::GradedSpace{I, <:Tuple}) where {I <: Sector} - dualV = isdual(V) - dualV == isdual(W) || _throw_not_subspace(V, W) - # single unrolled pass: subtract and validate non-negativity at once - newdims = map((dV, dW) -> dV < dW ? _throw_not_subspace(V, W) : dV - dW, V.dims, W.dims) - return typeof(V)(newdims, dualV) -end -# combiner for `⊖`, carrying the spaces to report which subspace condition was violated +# `⊖` as a callable, carrying the spaces to report which subspace condition was violated: +# as a combiner it subtracts and validates in one pass, and as an unmatched handler it rejects +# a sector of `W` that does not occur in `V` struct SubtractDims{S} V::S W::S end (f::SubtractDims)(dV, dW) = dV < dW ? _throw_not_subspace(f.V, f.W) : dV - dW -_unmatched1(::SubtractDims, d) = d # sector only in `V`: keep -_unmatched2(f::SubtractDims, d) = _throw_not_subspace(f.V, f.W) # sector only in `W`: not a subspace -_mergelength(::SubtractDims, n1, n2) = n1 +(f::SubtractDims)(d) = _throw_not_subspace(f.V, f.W) -function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) where {I <: Sector} +function ⊖(V::TupleGradedSpace{I}, W::TupleGradedSpace{I}) where {I <: Sector} + dualV = isdual(V) + dualV == isdual(W) || _throw_not_subspace(V, W) + newdims = map(SubtractDims(V, W), V.dims, W.dims) # single unrolled pass + return typeof(V)(newdims, dualV) +end +function ⊖(V::DictGradedSpace{I}, W::DictGradedSpace{I}) where {I <: Sector} dualV = isdual(V) dualV == isdual(W) || _throw_not_subspace(V, W) - return typeof(V)(mergewith(SubtractDims(V, W), V.dims, W.dims), dualV) + subtract = SubtractDims(V, W) + return typeof(V)(_sortedmerge(subtract, identity, subtract, V.dims, W.dims), dualV) end -function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} +function fuse(V₁::DictGradedSpace{I}, V₂::DictGradedSpace{I}) where {I <: Sector} acc = Dict{I, Int}() # Accumulation into Dict is more efficient than repeated insertion in sorted vector for (a, da) in blockdims(V₁), (b, db) in blockdims(V₂) dab = da * db @@ -205,7 +209,7 @@ function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorD perm = sortperm(ks0) return typeof(V₁)(SectorDict{I, Int}(ks0[perm], vs0[perm]), false) end -function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} +function fuse(V₁::TupleGradedSpace{I, N}, V₂::TupleGradedSpace{I, N}) where {I <: Sector, N} vals = values(I) dual1, dual2 = isdual(V₁), isdual(V₂) newdims = zeros(Int, N) @@ -229,34 +233,34 @@ function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{ return typeof(V₁)(ntuple(i -> @inbounds(newdims[i]), Val(N)), false) end -function infimum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} +function infimum(V₁::TupleGradedSpace{I}, V₂::TupleGradedSpace{I}) where {I <: Sector} Visdual = isdual(V₁) Visdual == isdual(V₂) || throw(SpaceMismatch("Infimum of space and dual space does not exist")) newdims = map(min, V₁.dims, V₂.dims) return typeof(V₁)(newdims, Visdual) end -function infimum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} +function infimum(V₁::DictGradedSpace{I}, V₂::DictGradedSpace{I}) where {I <: Sector} Visdual = isdual(V₁) Visdual == isdual(V₂) || throw(SpaceMismatch("Infimum of space and dual space does not exist")) return typeof(V₁)(mergewith(min, V₁.dims, V₂.dims), Visdual) end -function supremum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector} +function supremum(V₁::TupleGradedSpace{I}, V₂::TupleGradedSpace{I}) where {I <: Sector} Visdual = isdual(V₁) Visdual == isdual(V₂) || throw(SpaceMismatch("Supremum of space and dual space does not exist")) newdims = map(max, V₁.dims, V₂.dims) return typeof(V₁)(newdims, Visdual) end -function supremum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector} +function supremum(V₁::DictGradedSpace{I}, V₂::DictGradedSpace{I}) where {I <: Sector} Visdual = isdual(V₁) Visdual == isdual(V₂) || throw(SpaceMismatch("Supremum of space and dual space does not exist")) return typeof(V₁)(mergewith(max, V₁.dims, V₂.dims), Visdual) end hassector(V::GradedSpace{I}, s::I) where {I <: Sector} = dim(V, s) != 0 -function sectors(V::GradedSpace{I, <:AbstractDict}) where {I <: Sector} +function sectors(V::DictGradedSpace{I}) where {I <: Sector} return SectorSet{I}(s -> isdual(V) ? dual(s) : s, keys(V.dims)) end -function sectors(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} +function sectors(V::TupleGradedSpace{I, N}) where {I <: Sector, N} return SectorSet{I}(Iterators.filter(n -> V.dims[n] != 0, 1:N)) do n return isdual(V) ? dual(values(I)[n]) : values(I)[n] end @@ -269,10 +273,10 @@ Return an iterator over the non-zero blocks of the graded space `V`. These blocks contain the `Sector`s and their corresponding degeneracy, i.e. the number of times the sector appears in the direct sum decomposition of `V`. """ -function blockdims(V::GradedSpace{I, <:AbstractDict}) where {I <: Sector} +function blockdims(V::DictGradedSpace{I}) where {I <: Sector} return ((isdual(V) ? dual(c) : c) => d for (c, d) in V.dims) end -function blockdims(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N} +function blockdims(V::TupleGradedSpace{I, N}) where {I <: Sector, N} vals = values(I) return ( (isdual(V) ? dual(vals[n]) : vals[n]) => V.dims[n] @@ -285,7 +289,7 @@ function Base.:(==)(V₁::GradedSpace, V₂::GradedSpace) return sectortype(V₁) == sectortype(V₂) && (V₁.dims == V₂.dims) && V₁.dual == V₂.dual end -function sectorhash(V::GradedSpace{I, NTuple{N, Int}}, h::UInt) where {I, N} +function sectorhash(V::TupleGradedSpace{I, N}, h::UInt) where {I, N} return hash(iszero.(V.dims), hash(isdual(V), h)) end function sectorequal(V₁::GradedSpace{I, D}, V₂::GradedSpace{I, D}) where {I, N, D <: NTuple{N, Int}} @@ -293,7 +297,7 @@ function sectorequal(V₁::GradedSpace{I, D}, V₂::GradedSpace{I, D}) where {I, return iszero(d₁) == iszero(d₂) end end -function sectorhash(V::GradedSpace{I, <:SectorDict}, h::UInt) where {I} +function sectorhash(V::DictGradedSpace{I}, h::UInt) where {I} return hash(keys(V.dims), hash(isdual(V), h)) end function sectorequal(V₁::GradedSpace{I, D}, V₂::GradedSpace{I, D}) where {I, D <: SectorDict} diff --git a/src/tensors/sectorvector.jl b/src/tensors/sectorvector.jl index 11be45826..5c11c66e3 100644 --- a/src/tensors/sectorvector.jl +++ b/src/tensors/sectorvector.jl @@ -53,12 +53,9 @@ Base.valtype(::Type{SectorVector{T, I, A}}) where {T, I, A} = SubArray{T, 1, A, Base.keys(v::SectorVector) = keys(v.structure) Base.values(v::SectorVector) = (v[c] for c in keys(v)) -function Base.pairs(v::SectorVector) - # `structure` is already sorted, so the dict can be built without repeated insertion - # note: `view` need not produce `valtype(v)`, e.g. for GPU arrays - vals = map(Base.Fix1(view, parent(v)), values(v.structure)) - return SectorDict{keytype(v), eltype(vals)}(copy(v.structure.keys), vals) -end +# lazy, like `blocks(::AbstractTensorMap)`: the blocks are views, and lookups go through +# `getindex`/`get`/`block` on the `SectorVector` itself rather than through this iterator +Base.pairs(v::SectorVector) = Base.Iterators.map(((c, r),) -> c => view(parent(v), r), v.structure) Base.get(v::SectorVector{<:Any, I}, key::I, default) where {I} = haskey(v, key) ? v[key] : default Base.haskey(v::SectorVector{<:Any, I}, key::I) where {I} = key in keys(v) From 1ea0e8c6878e34f1ac7bed9e70075fbe53b3ffa3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 15 Sep 2026 18:38:02 -0400 Subject: [PATCH 34/36] Apply batched suggestions from code review [skip ci] Co-authored-by: Jutho --- docs/src/Changelog.md | 2 +- docs/src/lib/spaces.md | 6 +++--- src/spaces/gradedspace.jl | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/src/Changelog.md b/docs/src/Changelog.md index fbbab7cb3..282799375 100644 --- a/docs/src/Changelog.md +++ b/docs/src/Changelog.md @@ -26,7 +26,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - For sector types with `GenericUnit` such that colorings are not unique, `GradedSpace`, `ProductSpace` and `HomSpace` now check for this compatibility. In particular, this prevents the construction of `TensorMap`s with incompatible colorings, which previously either errored or produced empty tensors inconsistently. ([#515](https://github.com/QuantumKitHub/TensorKit.jl/pull/515)) ### Deprecated -- The type alias `ZNSpace{N}` is deprecated in favour of `Vect[ZNIrrep{N}]`: a type alias cannot compute the storage type from `N`, so the two only agree for small `N`. ([#511](https://github.com/QuantumKitHub/TensorKit.jl/pull/511)) +- The type alias `ZNSpace{N}` is deprecated in favour of `Vect[ZNIrrep{N}]` or `Rep[ℤ{N}]`: a type alias cannot compute the storage type from `N`, so the two only agree for small `N`. ([#511](https://github.com/QuantumKitHub/TensorKit.jl/pull/511)) ### Removed diff --git a/docs/src/lib/spaces.md b/docs/src/lib/spaces.md index 228c98726..d2bedc868 100644 --- a/docs/src/lib/spaces.md +++ b/docs/src/lib/spaces.md @@ -46,9 +46,9 @@ TensorKit.sectorstoragetype In this respect, there are also a number of type aliases for the `GradedSpace` types associated with the most common sectors, namely ```julia -const Z2Space = Vect[ZNIrrep{2}] -const Z3Space = Vect[ZNIrrep{3}] -const Z4Space = Vect[ZNIrrep{4}] +const Z2Space = Rep[ℤ{2}] +const Z3Space = Rep[ℤ{3}] +const Z4Space = Rep[ℤ{4}] const U1Space = Rep[U₁] const CU1Space = Rep[CU₁] const SU2Space = Rep[SU₂] diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index d46137043..7397e7be5 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -23,8 +23,8 @@ sector `s::I` can be transformed into an index via `s == getindex(values(I), i)` is used to store the non-zero degeneracy dimensions with the corresponding sector as key. The parameter `D` is hidden from the user and should typically be of no concern. -The concrete type `GradedSpace{I,D}` with correct `D` can be obtained as `Vect[I]`, or if -`I == Irrep[G]` for some `G<:Group`, as `Rep[G]`. +The concrete type `GradedSpace{I,D}` with correct `D` can be obtained as `Vect[I]`, +or as `Rep[G]` if `I == Irrep[G]` for some `G<:Group`. """ struct GradedSpace{I <: Sector, D} <: ElementarySpace dims::D From f26eccc05e312047c2e2cd602344716519414949 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 16 Sep 2026 08:39:40 -0400 Subject: [PATCH 35/36] review: drop StaticLength, simplify truncate_space, avoid re-sorting inds Address Jutho's final review round on #511: - replace `StaticLength(N)` with plain `N` (constant propagation makes the old TupleTools artefact unnecessary now that N is already a type parameter), and drop the now-dead `TupleTools: StaticLength` imports. - rewrite `truncate_space(::TupleGradedSpace, inds)` as a pure `ntuple` closure instead of mutating a `MutableNTuple`, guarding zero-dimension sectors before indexing into `inds` (which only holds keys for sectors with nonzero dimension). - `truncate_space(::DictGradedSpace, inds)` no longer collects and sorts: `inds` (whether a `SectorDict` or a `SectorVector`, depending on the truncation strategy) already iterates in sorted order by sector. - drop the now-unused `MutableNTuple`/`findindex` imports in the factorizations submodule. Co-Authored-By: Claude Sonnet 5 --- src/TensorKit.jl | 1 - src/factorizations/factorizations.jl | 4 +--- src/factorizations/truncation.jl | 23 ++++++++++++----------- src/spaces/gradedspace.jl | 2 +- 4 files changed, 14 insertions(+), 16 deletions(-) diff --git a/src/TensorKit.jl b/src/TensorKit.jl index a03e1af1a..41a7d8ae1 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -107,7 +107,6 @@ export empty_globalcaches! # Imports #--------- using TupleTools -using TupleTools: StaticLength using Strided diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 12b8b4437..e49cc8b29 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -7,7 +7,7 @@ export copy_oftype, factorisation_scalartype, one!, truncspace using ..TensorKit using ..TensorKit: AdjointTensorMap, DictGradedSpace, SectorDict, SectorVector, - TupleGradedSpace, findindex, + TupleGradedSpace, blocktype, foreachblock, one!, similar_diagonal, similarstoragetype using ..TensorKit: GLOBAL_TIMER @@ -19,8 +19,6 @@ using LinearAlgebra: LinearAlgebra, BlasFloat, Diagonal, using TensorOperations: Index2Tuple -using TupleTools: MutableNTuple, StaticLength - using MatrixAlgebraKit import MatrixAlgebraKit as MAK using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, DiagonalAlgorithm diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 3e6ef8ee8..bd6956440 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -41,27 +41,28 @@ end function truncate_space(V::TupleGradedSpace{I, N}, inds) where {I <: Sector, N} @assert !isdual(V) vals = values(I) - newdims = MutableNTuple(ntuple(_ -> 0, StaticLength(N))) - for (c, ind) in pairs(inds) - d = dim(V, c) - n_write = findindex(vals, c) - @inbounds newdims[n_write] = _blocklength(d, ind) + # `inds` has a key for every sector with nonzero `dim(V, c)`, so zero-dim sectors are skipped rather than looked up. + newdims = ntuple(N) do n + d = V.dims[n] + return iszero(d) ? 0 : _blocklength(d, inds[vals[n]]) end - return typeof(V)(Tuple(newdims), false) + return typeof(V)(newdims, false) end function truncate_space(V::DictGradedSpace{I}, inds) where {I <: Sector} @assert !isdual(V) - ks, vs = Vector{I}(), Vector{Int}() # accumulate and sort once at the end + # `inds` (a `SectorDict` or `SectorVector`, depending on the truncation strategy) always + # iterates in sorted order by sector, so no need to sort again here. + ks, vs = Vector{I}(undef, 0), Vector{Int}(undef, 0) + sizehint!(ks, length(inds)) + sizehint!(vs, length(inds)) for (c, ind) in pairs(inds) - d = dim(V, c) - len = _blocklength(d, ind) + len = _blocklength(dim(V, c), ind) if !iszero(len) push!(ks, c) push!(vs, len) end end - perm = sortperm(ks) - return typeof(V)(SectorDict{I, Int}(ks[perm], vs[perm]), false) + return typeof(V)(SectorDict{I, Int}(ks, vs), false) end function truncate_domain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 7397e7be5..67169a764 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -47,7 +47,7 @@ end function GradedSpace{I, NTuple{N, Int}}(dims; dual::Bool = false) where {I, N} @assert N <= 64 "the `UInt64` bitmask tracking which sectors have been set holds 64 bits" - d = TupleTools.MutableNTuple(ntuple(Returns(0), StaticLength(N))) + d = TupleTools.MutableNTuple(ntuple(Returns(0), N)) mask = zero(UInt64) for (c, dc) in dims k = convert(I, c) From 5521fac06a4a6b48062b7a4b0ce057d3c0ed30e0 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 16 Sep 2026 10:23:21 -0400 Subject: [PATCH 36/36] test: fix pairs(::SectorVector) comparison against Dict `pairs(::SectorVector)` is intentionally lazy (a `Base.Generator`), so comparing it directly to a `Dict` with `==` silently evaluates to `false` rather than erroring, since no `==` is defined between those unrelated iterator types. Collect it before comparing. Co-Authored-By: Claude Sonnet 5 --- test/tensors/sectorvector.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/tensors/sectorvector.jl b/test/tensors/sectorvector.jl index 1c8b268e6..22832404c 100644 --- a/test/tensors/sectorvector.jl +++ b/test/tensors/sectorvector.jl @@ -83,7 +83,7 @@ spacelist = ( end @test collect(@testinferred(keys(v))) == collect(sectors(V)) @test collect(values(v)) == [v[c] for c in keys(v)] - @test pairs(v) == Dict(c => v[c] for c in keys(v)) + @test Dict(pairs(v)) == Dict(c => v[c] for c in keys(v)) end @timedtestset "VectorInterface ($T)" begin