diff --git a/Project.toml b/Project.toml index 192d9b22..1c571a84 100644 --- a/Project.toml +++ b/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 and contributors"] [workspace] diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl index b89f7674..2e9285c6 100644 --- a/src/tensoralgebra.jl +++ b/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 `_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 diff --git a/test/test_gradedarraysext.jl b/test/test_gradedarraysext.jl index c132e6d7..aac4facb 100644 --- a/test/test_gradedarraysext.jl +++ b/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 diff --git a/test/test_tensorkitext.jl b/test/test_tensorkitext.jl index 9ab69884..ebe6258d 100644 --- a/test/test_tensorkitext.jl +++ b/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,22 +103,22 @@ 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 @@ -127,19 +126,19 @@ using Test: @test, @test_throws, @testset # 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