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

@@ -808,35 +808,56 @@ an empty domain). The index axes select the backend: dense ranges give an
 `TensorMap`. `a` is indexed positionally in the order
 `(codomain_inds..., domain_inds...)`.
 
+When `a` has more axes than `codomain_inds` and `domain_inds` account for, each
+surplus trailing axis is an auxiliary leg that the backend derived to make the
+result symmetry-allowed (for example a flux-canceling charge leg for a
+charge-shifting operator). Such axes are returned as named dimensions with
+freshly generated names, which the caller can read off the result.
+
 `TensorAlgebra.tryproject` returns `nothing` instead of throwing, and
 `TensorAlgebra.unchecked_project` skips the verification.
 """
 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.
+# Attach `input_names` to the leading axes of the `projected` array and mint a fresh unique name for
+# each surplus axis. A surplus axis is an auxiliary leg the backend derived 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 project_nameddims(projected, input_names)
+    isnothing(projected) && return nothing
+    aux_names = ntuple(
+        _ -> uniquename(first(input_names)),
+        TA.ndims(projected) - length(input_names)
+    )
+    return nameddims(projected, (input_names..., aux_names...))
+end
+
+# Strip to `a`'s axes, lower to the TensorAlgebra verb, and reattach the names. 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.
 for f in (:project, :tryproject, :unchecked_project)
     @eval begin
         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)...))
+            projected = TA.$f(a, unnamed.(codomain_inds), unnamed.(domain_inds); kwargs...)
+            return project_nameddims(
+                projected,
+                (name.(codomain_inds)..., name.(domain_inds)...)
+            )
         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))
+            projected = TA.$f(a, (), unnamed.(domain_inds); kwargs...)
+            return project_nameddims(projected, name.(domain_inds))
         end
         function TA.$f(
                 a::AbstractArray, inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}};

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, unchecked_project
 using TensorKitSectors: FermionNumber
 using Test: @test, @testset
 
@@ -63,3 +63,28 @@ using Test: @test, @testset
     @test length(inds(fl)) == 3
     @test eltype(fl) == elt
 end
+
+# `project` derives 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 `project` mint the flux-canceling aux leg
+    op = project(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

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