Add named project for building symmetry-restricted tensors

Project.toml

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

src/tensoralgebra.jl

@@ -508,3 +508,74 @@ for f in MATRIX_FUNCTIONS
         end
     end
 end
+
+#
+# Projection into a symmetry-restricted named tensor.
+#
+
+# Shared implementation: strip to `a`'s axes, lower to the TensorAlgebra hook, and reattach the
+# names. The dispatch entries below feed this from either a codomain/domain split or a flat list
+# (an empty domain). Two split entries per function 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.
+function project_nameddims(a, codomain_inds, domain_inds)
+    raw = TA.project(a, unnamed.(codomain_inds), unnamed.(domain_inds))
+    return nameddims(raw, (name.(codomain_inds)..., name.(domain_inds)...))
+end
+function checked_project_nameddims(a, codomain_inds, domain_inds; kwargs...)
+    raw = TA.checked_project(a, unnamed.(codomain_inds), unnamed.(domain_inds); kwargs...)
+    return nameddims(raw, (name.(codomain_inds)..., name.(domain_inds)...))
+end
+
+"""
+    TensorAlgebra.project(a::AbstractArray, codomain_inds, domain_inds) -> t
+    TensorAlgebra.project(a::AbstractArray, inds) -> t
+
+Build a named tensor from the dense array `a` by projecting it into the
+symmetry-restricted space described by the indices. 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...)`.
+
+See `TensorAlgebra.checked_project` for a version that verifies nothing outside
+the symmetry-allowed blocks was discarded.
+"""
+function TA.project(
+        a::AbstractArray,
+        codomain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}},
+        domain_inds::Tuple{Vararg{NamedUnitRange}}
+    )
+    return project_nameddims(a, codomain_inds, domain_inds)
+end
+function TA.project(
+        a::AbstractArray,
+        codomain_inds::Tuple{},
+        domain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}
+    )
+    return project_nameddims(a, codomain_inds, domain_inds)
+end
+function TA.project(a::AbstractArray, inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}})
+    return project_nameddims(a, inds, ())
+end
+
+function TA.checked_project(
+        a::AbstractArray,
+        codomain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}},
+        domain_inds::Tuple{Vararg{NamedUnitRange}}; kwargs...
+    )
+    return checked_project_nameddims(a, codomain_inds, domain_inds; kwargs...)
+end
+function TA.checked_project(
+        a::AbstractArray,
+        codomain_inds::Tuple{},
+        domain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}; kwargs...
+    )
+    return checked_project_nameddims(a, codomain_inds, domain_inds; kwargs...)
+end
+function TA.checked_project(
+        a::AbstractArray, inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}; kwargs...
+    )
+    return checked_project_nameddims(a, inds, (); kwargs...)
+end

test/test_tensoralgebra.jl

@@ -1,11 +1,12 @@
-using ITensorBase:
-    ITensorBase, dimnames, inds, namedoneto, replacedimnames, uniquename, unname, unnamed
+using ITensorBase: ITensorBase, Index, dimnames, inds, name, namedoneto, prime,
+    replacedimnames, uniquename, unname, unnamed
 using LinearAlgebra: LinearAlgebra, norm
 using MatrixAlgebraKit: left_null, left_orth, left_polar, lq_compact, lq_full, qr_compact,
     qr_full, right_null, right_orth, right_polar, svd_compact, svd_trunc
 using StableRNGs: StableRNG
 using TensorAlgebra.MatrixAlgebra: gram_eigh_full, gram_eigh_full_with_pinv
-using TensorAlgebra: TensorAlgebra, contract, matricize, unmatricize
+using TensorAlgebra:
+    TensorAlgebra, checked_project, contract, matricize, project, unmatricize
 using Test: @test, @test_broken, @testset
 
 @testset "TensorAlgebra (eltype=$(elt))" for elt in
@@ -156,4 +157,24 @@ using Test: @test, @test_broken, @testset
             @test YXmat ≈ LinearAlgebra.I(size(YXmat, 1))
         end
     end
+    @testset "project" begin
+        i = Index(2)
+        Sz = elt[0.5 0; 0 -0.5]
+        # the three-argument form builds an operator from the codomain/domain split
+        top = project(Sz, (prime(i),), (i,))
+        @test eltype(top) === elt
+        @test Set(dimnames(top)) == Set(name.((prime(i), i)))
+        @test unname(top, (prime(i), i)) == Sz
+        # `checked_project` accepts the (for dense, always exact) projection
+        @test unname(checked_project(Sz, (prime(i),), (i,)), (prime(i), i)) == Sz
+        # the two-argument form builds a state (empty domain)
+        v = elt[1, 0]
+        s = project(v, (i,))
+        @test dimnames(s) == [name(i)]
+        @test unname(s, (i,)) == v
+        # the empty-codomain form builds an all-domain tensor (mirror of the state)
+        bra = project(v, (), (i,))
+        @test dimnames(bra) == [name(i)]
+        @test unname(bra, (i,)) == v
+    end
 end

test/test_tensorkitext.jl

@@ -1,9 +1,10 @@
-using ITensorBase: ITensorBase, Index, dimnames, name, unnamed
+using ITensorBase: ITensorBase, Index, dimnames, name, prime, unnamed
 using LinearAlgebra: norm
 using MatrixAlgebraKit: qr_compact, svd_compact
 using StableRNGs: StableRNG
+using TensorAlgebra: checked_project, project
 using TensorKit: TensorKit, @tensor, AbstractTensorMap, SU2Irrep, U1Irrep, Vect, dim, dual,
-    scalar, space, ⊗
+    scalar, space, ←, ⊗
 using Test: @test, @test_throws, @testset
 
 # A native TensorKit space flows into `Index`, so an `ITensor` wraps a `TensorMap` directly.
@@ -68,4 +69,36 @@ using Test: @test, @test_throws, @testset
         q, r = qr_compact(a3, (i,), (j, k))
         @test sca((q * r) * w) ≈ sca(a3 * w)
     end
+
+    # `project` builds a `TensorMap`-backed operator/state from a dense basis matrix: the index
+    # spaces select the backend, so the same call that yields an `Array` on dense indices yields a
+    # `TensorMap` here, keeping only the symmetry-allowed blocks.
+    @testset "project" begin
+        W = Vect[U1Irrep](0 => 1, 1 => 1)
+        w = Index(W)
+        Sz = elt[0.5 0; 0 -0.5]
+        Sx = elt[0 0.5; 0.5 0]
+
+        top = project(Sz, (prime(w),), (w,))
+        @test unnamed(top) isa AbstractTensorMap
+        @test space(unnamed(top)) == (W ← W)
+        @test Set(dimnames(top)) == Set(name.((prime(w), w)))
+
+        # a charge-breaking operator is projected to zero; `checked_project` rejects the discard
+        @test norm(unnamed(project(Sx, (prime(w),), (w,)))) == 0
+        @test_throws InexactError checked_project(Sx, (prime(w),), (w,); atol = 0, rtol = 0)
+
+        # the two-argument form builds an all-codomain state; only the trivial-charge component
+        # of the dense vector survives
+        pv = project(elt[1, 0], (w,))
+        @test unnamed(pv) isa AbstractTensorMap
+        @test norm(unnamed(pv)) ≈ 1
+        @test norm(unnamed(project(elt[0, 1], (w,)))) == 0
+
+        # 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 Set(dimnames(cobra)) == Set((name(w),))
+    end
 end