Derive named auxiliary legs in project

Project.toml

@@ -1,6 +1,6 @@
 name = "ITensorBase"
 uuid = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
-version = "0.13.6"
+version = "0.13.7"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]

src/tensoralgebra.jl

@@ -793,56 +793,116 @@ end
 # Projection into a symmetry-restricted named tensor.
 #
 
+# Attach `input_names` to the leading axes of the `projected` array and mint a fresh unique name for
+# a trailing surplus axis if the backend derived one (it appends at most one, as the last domain
+# axis). A surplus axis is an auxiliary leg the backend adds so the result is symmetry-allowed (for
+# example a flux-canceling charge leg for a charge-shifting operator), and naming it returns it as a
+# dimension the caller can read off the result. A `nothing` (from `tryproject`) passes straight
+# through.
+function name_projected(projected, input_names)
+    isnothing(projected) && return nothing
+    aux_names = ntuple(
+        _ -> uniquename(eltype(input_names)),
+        TA.ndims(projected) - length(input_names)
+    )
+    return nameddims(projected, (input_names..., aux_names...))
+end
+
+# Each `<verb>_nameddims` runs the named-index layer of a `TensorAlgebra` verb: strip the axes to
+# their unnamed ranges, lower to the dense verb, and reattach the names, minting a name for the aux
+# leg the backend may append (see `name_projected`). The one body also covers an empty codomain,
+# since `unnamed.(())` and `name.(())` are both `()`, so the all-domain (co-state) case needs no
+# separate path.
+function project_nameddims(a, codomain_inds, domain_inds; kwargs...)
+    projected = TA.project(a, unnamed.(codomain_inds), unnamed.(domain_inds); kwargs...)
+    return name_projected(projected, (name.(codomain_inds)..., name.(domain_inds)...))
+end
+function tryproject_nameddims(a, codomain_inds, domain_inds; kwargs...)
+    projected = TA.tryproject(a, unnamed.(codomain_inds), unnamed.(domain_inds); kwargs...)
+    return name_projected(projected, (name.(codomain_inds)..., name.(domain_inds)...))
+end
+function unchecked_project_nameddims(a, codomain_inds, domain_inds; kwargs...)
+    projected =
+        TA.unchecked_project(a, unnamed.(codomain_inds), unnamed.(domain_inds); kwargs...)
+    return name_projected(projected, (name.(codomain_inds)..., name.(domain_inds)...))
+end
+
+# Shared body for the named-index `project` family docstrings. Each function's summary states its
+# own verification behavior; this describes what all three have in common.
+const _project_named_body = """
+The three-argument form takes an explicit codomain/domain split (an operator), and the
+two-argument form a flat list of indices (a state, i.e. an empty domain). The index axes select
+the backend: dense ranges give an `Array`, graded ranges a block-sparse array, and TensorKit
+spaces a `TensorMap`. `a` is indexed positionally in the order `(codomain_inds..., domain_inds...)`.
+
+When `a` carries one more axis than `codomain_inds` and `domain_inds` account for, that trailing
+surplus axis is an auxiliary leg the backend derived to make the result symmetry-allowed (for
+example a flux-canceling charge leg for a charge-shifting operator). It is returned as a named
+dimension with a freshly generated name the caller can read off the result.
 """
+
+const _project_named_docstring = """
     TensorAlgebra.project(a::AbstractArray, codomain_inds, domain_inds; kwargs...) -> t
     TensorAlgebra.project(a::AbstractArray, inds; kwargs...) -> t
 
-Build a named tensor from the dense array `a` by projecting it into the
-symmetry-restricted space described by the indices, verifying that only a
-negligible component of `a` is discarded and throwing an `InexactError`
-otherwise (keyword arguments are forwarded to the `isapprox` tolerance
-check). The three-argument form takes an explicit codomain/domain split (an
-operator); the two-argument form takes a flat list of indices (a state, i.e.
-an empty domain). The index axes select the backend: dense ranges give an
-`Array`, graded ranges a block-sparse array, and TensorKit spaces a
-`TensorMap`. `a` is indexed positionally in the order
-`(codomain_inds..., domain_inds...)`.
-
-`TensorAlgebra.tryproject` returns `nothing` instead of throwing, and
-`TensorAlgebra.unchecked_project` skips the verification.
+Build a named tensor by projecting the dense array `a` into the symmetry-restricted space
+described by the indices, verifying that only a negligible component of `a` is discarded and
+throwing an `InexactError` otherwise (keyword arguments are forwarded to the `isapprox` tolerance
+check).
+
+$(_project_named_body)
+See also `TensorAlgebra.tryproject` and `TensorAlgebra.unchecked_project`.
+"""
+
+const _tryproject_named_docstring = """
+    TensorAlgebra.tryproject(a::AbstractArray, codomain_inds, domain_inds; kwargs...) -> Union{t, Nothing}
+    TensorAlgebra.tryproject(a::AbstractArray, inds; kwargs...) -> Union{t, Nothing}
+
+Like `TensorAlgebra.project`, but return `nothing` instead of throwing when a non-negligible
+component of `a` would be discarded (keyword arguments are forwarded to the `isapprox` tolerance
+check).
+
+$(_project_named_body)
+See also `TensorAlgebra.project` and `TensorAlgebra.unchecked_project`.
 """
-TA.project
 
-# Strip to `a`'s axes, lower to the TensorAlgebra verb, and reattach the names (passing
-# through a `nothing` from `tryproject`). Two split entries per verb read the index type from
-# whichever side is non-empty (an empty codomain is the all-domain case, the mirror of the
-# empty-domain state), so neither an empty codomain nor an empty domain falls through to the
-# unnamed-axis generic; the flat (state) form forwards to the split form.
+const _unchecked_project_named_docstring = """
+    TensorAlgebra.unchecked_project(a::AbstractArray, codomain_inds, domain_inds; kwargs...) -> t
+    TensorAlgebra.unchecked_project(a::AbstractArray, inds; kwargs...) -> t
+
+Like `TensorAlgebra.project`, but skip the verification: components of `a` outside the
+symmetry-allowed structure are dropped without inspection.
+
+$(_project_named_body)
+See also `TensorAlgebra.project` and `TensorAlgebra.tryproject`.
+"""
+
+# Forward each named-index signature to its worker, attaching the family docstring to the split
+# form. Two split entries per verb so an empty codomain or an empty domain still selects this
+# overload instead of the unnamed-axis generic. The flat (state) form forwards with an empty domain.
 for f in (:project, :tryproject, :unchecked_project)
+    fnamed = Symbol(f, :_nameddims)
+    doc = Symbol("_", f, "_named_docstring")
     @eval begin
-        function TA.$f(
+        @doc $doc function TA.$f(
                 a::AbstractArray,
                 codomain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}},
                 domain_inds::Tuple{Vararg{NamedUnitRange}}; kwargs...
             )
-            raw = TA.$f(a, unnamed.(codomain_inds), unnamed.(domain_inds); kwargs...)
-            isnothing(raw) && return nothing
-            return nameddims(raw, (name.(codomain_inds)..., name.(domain_inds)...))
+            return $fnamed(a, codomain_inds, domain_inds; kwargs...)
         end
         function TA.$f(
                 a::AbstractArray,
                 codomain_inds::Tuple{},
                 domain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}; kwargs...
             )
-            raw = TA.$f(a, (), unnamed.(domain_inds); kwargs...)
-            isnothing(raw) && return nothing
-            return nameddims(raw, name.(domain_inds))
+            return $fnamed(a, codomain_inds, domain_inds; kwargs...)
         end
         function TA.$f(
                 a::AbstractArray, inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}};
                 kwargs...
             )
-            return TA.$f(a, inds, (); kwargs...)
+            return $fnamed(a, inds, (); kwargs...)
         end
     end
 end

test/test_gradedarraysext.jl

@@ -1,7 +1,7 @@
 using GradedArrays: U1, sectors
-using ITensorBase: ITensorBase, Index, inds, space
+using ITensorBase: ITensorBase, Index, inds, prime, space
 using StableRNGs: StableRNG
-using TensorAlgebra: TensorAlgebra, isdual
+using TensorAlgebra: TensorAlgebra, isdual, project, tryproject, unchecked_project
 using TensorKitSectors: FermionNumber
 using Test: @test, @testset
 
@@ -63,3 +63,30 @@ using Test: @test, @testset
     @test length(inds(fl)) == 3
     @test eltype(fl) == elt
 end
+
+# `project` and its siblings derive the same kind of auxiliary leg: a trailing surplus axis on the
+# dense array becomes a named aux dimension carrying the operator's flux, so a charge-shifting
+# operator stays symmetry-allowed instead of being projected away.
+@testset "project derives a named auxiliary leg (eltype = $elt)" for elt in
+    (
+        Float64,
+        ComplexF64,
+    )
+    s = Index([U1(0) => 1, U1(1) => 1]; tags = "s")
+    cdag = elt[0 0; 1 0]   # raising operator, flux +1
+
+    # without a surplus axis the charge-shifting operator has nothing to carry its flux
+    @test iszero(unchecked_project(cdag, (prime(s),), (s,)))
+
+    # reshaping to a trailing length-1 axis lets each verb mint the flux-canceling aux leg
+    @testset "$f" for f in (project, tryproject, unchecked_project)
+        op = f(reshape(cdag, (2, 2, 1)), (prime(s),), (s,))
+        @test length(inds(op)) == 3
+        @test !iszero(op)
+        @test eltype(op) == elt
+        aux = only(setdiff(collect(inds(op)), [prime(s), s]))
+        @test length(aux) == 1
+        @test isdual(aux)                          # dualized, in the domain
+        @test only(sectors(space(aux))) == U1(1)   # carries the operator's flux
+    end
+end

test/test_tensorkitext.jl

@@ -3,8 +3,7 @@ using LinearAlgebra: norm
 using MatrixAlgebraKit: qr_compact, svd_compact
 using StableRNGs: StableRNG
 using TensorAlgebra: TensorAlgebra, project, unchecked_project
-using TensorKit: TensorKit, @tensor, AbstractTensorMap, SU2Irrep, U1Irrep, Vect, dim, dual,
-    scalar, space, ←, ⊗
+using TensorKit: TensorKit as TK, @tensor, AbstractTensorMap, SU2Irrep, U1Irrep, Vect, ←, ⊗
 using Test: @test, @test_throws, @testset
 
 # A native TensorKit space flows into `Index`, so an `ITensor` wraps a `TensorMap` directly.
@@ -33,13 +32,13 @@ using Test: @test, @test_throws, @testset
 
         # `conj(index)` round-trips to an `Index` carrying the dual space, same name.
         @test conj(i) isa Index
-        @test unnamed(conj(i)) == dual(Vi)
+        @test unnamed(conj(i)) == TK.dual(Vi)
         @test name(conj(i)) == name(i)
 
         # `isdual`/`dual` forward to the underlying space.
         @test TensorAlgebra.isdual(i) == false
         @test TensorAlgebra.isdual(conj(i)) == true
-        @test unnamed(TensorAlgebra.dual(i)) == dual(Vi)
+        @test unnamed(TensorAlgebra.dual(i)) == TK.dual(Vi)
         @test name(TensorAlgebra.dual(i)) == name(i)
 
         # Equality is dual-insensitive: an index equals its dual (same name, same ungraded
@@ -54,7 +53,7 @@ using Test: @test, @test_throws, @testset
         # Cold-start construction wraps a `TensorMap`; size/eltype report dense values.
         a = randn(rng, elt, i, j)
         @test unnamed(a) isa AbstractTensorMap
-        @test size(a) == (dim(Vi), dim(Vj))
+        @test size(a) == (TK.dim(Vi), TK.dim(Vj))
         @test eltype(a) == elt
         @test norm(unnamed(zeros(elt, i, j))) == 0
 
@@ -64,7 +63,7 @@ using Test: @test, @test_throws, @testset
         @test Set(dimnames(c)) == Set(name.((i, k)))
         ta, tb, gc = unnamed(a), unnamed(b), unnamed(c)
         @tensor ref[vi; vk] := ta[vi, vj] * tb[vj, vk]
-        @test space(ref) == space(gc)
+        @test TK.space(ref) == TK.space(gc)
         @test ref ≈ gc
 
         # Linear-combination broadcast lowers to `bipermutedimsopadd!`; element-wise errors.
@@ -78,7 +77,7 @@ using Test: @test, @test_throws, @testset
         # Checked by full contraction to a scalar, which is bipartition-independent.
         a3 = randn(rng, elt, i, j, k)
         w = randn(rng, elt, conj(i), conj(j), conj(k))
-        sca(x) = scalar(unnamed(x))
+        sca(x) = TK.scalar(unnamed(x))
         u, s, v = svd_compact(a3, (i,), (j, k))
         @test sca((u * s * v) * w) ≈ sca(a3 * w)
         q, r = qr_compact(a3, (i,), (j, k))
@@ -104,42 +103,42 @@ using Test: @test, @test_throws, @testset
         # view, the same as a flat graded array would have axes `(Vi, dual(Vj))`.
         m = randn(rng, elt, (i,), (j,))
         @test unnamed(m) isa AbstractTensorMap
-        @test space(unnamed(m)) == (Vi ← Vj)
-        @test space(unnamed(m), 1) == Vi
-        @test space(unnamed(m), 2) == dual(Vj)
+        @test TK.space(unnamed(m)) == (Vi ← Vj)
+        @test TK.space(unnamed(m), 1) == Vi
+        @test TK.space(unnamed(m), 2) == TK.dual(Vj)
         @test unnamed(rand(rng, elt, (i,), (j,))) isa AbstractTensorMap
         @test norm(unnamed(zeros(elt, (i,), (j,)))) == 0
         # The friendly forms agree with the underlying `TensorAlgebra` map hooks.
-        @test space(unnamed(TensorAlgebra.randn_map(elt, (i, j), (k,)))) ==
-            space(unnamed(randn(elt, (i, j), (k,))))
+        @test TK.space(unnamed(TensorAlgebra.randn_map(elt, (i, j), (k,)))) ==
+            TK.space(unnamed(randn(elt, (i, j), (k,))))
         # An empty codomain builds an all-domain `TensorMap`, the mirror of an empty domain. The
         # space type is read from the domain, since the empty codomain carries none. An all-empty
         # split has no map meaning and errors rather than recursing.
         cd = randn(rng, elt, (), (j,))
         @test unnamed(cd) isa AbstractTensorMap
-        @test space(unnamed(cd)) == (one(Vj) ← Vj)
+        @test TK.space(unnamed(cd)) == (one(Vj) ← Vj)
         @test dimnames(cd) == [name(j)]
-        @test space(unnamed(zeros(elt, (), (j,)))) == (one(Vj) ← Vj)
+        @test TK.space(unnamed(zeros(elt, (), (j,)))) == (one(Vj) ← Vj)
         @test_throws MethodError randn(rng, elt, (), ())
 
         # `aligndims` reorders a `TensorMap`-backed tensor. The flat form gives an all-codomain
         # result and the map form re-expresses the requested codomain/domain split, both
         # carrying each index with its arrow to the new position.
         mf = aligndims(m, (j, i))
         @test dimnames(mf) == [name(j), name(i)]
-        @test space(unnamed(mf), 1) == dual(Vj)
-        @test space(unnamed(mf), 2) == Vi
+        @test TK.space(unnamed(mf), 1) == TK.dual(Vj)
+        @test TK.space(unnamed(mf), 2) == Vi
         md = aligndims(m, (j,), (i,))
         @test dimnames(md) == [name(j), name(i)]
-        @test space(unnamed(md)) == (dual(Vj) ← dual(Vi))
-        @test space(unnamed(md), 1) == dual(Vj)
-        @test space(unnamed(md), 2) == Vi
+        @test TK.space(unnamed(md)) == (TK.dual(Vj) ← TK.dual(Vi))
+        @test TK.space(unnamed(md), 1) == TK.dual(Vj)
+        @test TK.space(unnamed(md), 2) == Vi
         # An empty codomain moves both indices into the domain, preserving the outward axes.
         me = aligndims(m, (), (i, j))
         @test dimnames(me) == [name(i), name(j)]
-        @test space(unnamed(me)) == (one(Vi) ← (dual(Vi) ⊗ Vj))
-        @test space(unnamed(me), 1) == Vi
-        @test space(unnamed(me), 2) == dual(Vj)
+        @test TK.space(unnamed(me)) == (one(Vi) ← (TK.dual(Vi) ⊗ Vj))
+        @test TK.space(unnamed(me), 1) == Vi
+        @test TK.space(unnamed(me), 2) == TK.dual(Vj)
     end
 
     # `project` builds a `TensorMap`-backed operator/state from a dense basis matrix: the index
@@ -153,7 +152,7 @@ using Test: @test, @test_throws, @testset
 
         top = project(Sz, (prime(w),), (w,))
         @test unnamed(top) isa AbstractTensorMap
-        @test space(unnamed(top)) == (W ← W)
+        @test TK.space(unnamed(top)) == (W ← W)
         @test Set(dimnames(top)) == Set(name.((prime(w), w)))
 
         # a charge-breaking operator is projected to zero by `unchecked_project`; the checked
@@ -171,7 +170,7 @@ using Test: @test, @test_throws, @testset
         # the empty-codomain form builds an all-domain `TensorMap` (the mirror case)
         cobra = project(elt[1, 0], (), (w,))
         @test unnamed(cobra) isa AbstractTensorMap
-        @test space(unnamed(cobra)) == (one(W) ← W)
+        @test TK.space(unnamed(cobra)) == (one(W) ← W)
         @test Set(dimnames(cobra)) == Set((name(w),))
     end
 end