diff --git a/Project.toml b/Project.toml index f4aea5f..d2840c0 100644 --- a/Project.toml +++ b/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 and contributors"] [workspace] diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl index e708b69..0dfebf3 100644 --- a/src/tensoralgebra.jl +++ b/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 diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl index 835eabd..3998e50 100644 --- a/test/test_tensoralgebra.jl +++ b/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 diff --git a/test/test_tensorkitext.jl b/test/test_tensorkitext.jl index 26930ff..c1be719 100644 --- a/test/test_tensorkitext.jl +++ b/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