diff --git a/Project.toml b/Project.toml index b1f463b0..4caf5d30 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ITensorBase" uuid = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" -version = "0.13.8" +version = "0.13.9" authors = ["ITensor developers and contributors"] [workspace] diff --git a/src/ITensorBase.jl b/src/ITensorBase.jl index 95c71fd4..6e3de08f 100644 --- a/src/ITensorBase.jl +++ b/src/ITensorBase.jl @@ -2,8 +2,11 @@ module ITensorBase export AbstractNamedTensor, NamedTensor, AbstractITensor, ITensor, Index, NamedUnitRange, aligndims, aligneddims, apply, commonind, commoninds, - dimnames, dimnametype, hascommoninds, id, inds, inputnames, mapinds, named, nameddims, - noncommonind, noncommoninds, noprime, operator, outputnames, prime, replaceinds, sim, + dimnames, dimnametype, hascommoninds, id, inds, inputaxes, inputinds, inputnames, + mapinds, + named, nameddims, noncommonind, noncommoninds, noprime, operator, outputaxes, + outputinds, + outputnames, prime, replaceinds, sim, similar_operator, state, trycommonind, trynoncommonind, tryuniqueind, uniqueind, uniqueinds, unioninds, uniquename diff --git a/src/namedtensoroperator.jl b/src/namedtensoroperator.jl index 6759baef..d8ca4393 100644 --- a/src/namedtensoroperator.jl +++ b/src/namedtensoroperator.jl @@ -74,6 +74,79 @@ See also [`outputnames`](@ref), [`operator`](@ref), [`apply`](@ref). """ inputnames(a::AbstractNamedTensor{DimName}) where {DimName} = DimName[] +""" + outputinds(a) + +The output (codomain) indices of an operator `a` — the space it maps onto — in the same +order as [`outputnames`](@ref). Together with [`inputinds`](@ref) these are the codomain +and domain of `a` viewed as a map, so `id(eltype(a), outputinds(a), inputinds(a))` rebuilds +an operator of the same shape. Compare with `outputnames`, which returns just the names. +A plain tensor is a trivial operator with no pairing, so its output indices are empty. + +# Examples + +```jldoctest +julia> op = operator(zeros(2, 2), ("i",), ("j",)); + +julia> outputinds(op) +1-element Vector{NamedUnitRange{String, Int64, Base.OneTo{Int64}}}: + named(Base.OneTo(2), "i") +``` + +See also [`outputnames`](@ref), [`inputinds`](@ref), [`outputaxes`](@ref), [`operator`](@ref). +""" +outputinds(a::AbstractNamedTensor) = inds(a)[dims(a, outputnames(a))] + +""" + outputaxes(a) + +The output (codomain) indices of an operator `a` as a `Tuple`, the tuple form of +[`outputinds`](@ref) (mirroring how `axes` relates to [`inds`](@ref)). The `Tuple` feeds +directly into constructors that take a `(codomain, domain)` pair of index tuples, such as +[`id`](@ref): `id(eltype(a), outputaxes(a), inputaxes(a))` rebuilds an operator of `a`'s +shape. + +See also [`outputinds`](@ref), [`inputaxes`](@ref), [`operator`](@ref). +""" +outputaxes(a::AbstractNamedTensor) = map(Base.Fix1(inds, a), Tuple(dims(a, outputnames(a)))) + +""" + inputinds(a) + +The input (domain) indices of an operator `a` — the space of states it acts on — in the +same order as [`inputnames`](@ref). This is the domain in the sense of TensorKit's `domain`: +the non-dual space a state occupies, so `a * randn(Tuple(inputinds(a)))` contracts. It is the +conjugate of the operator's input legs as they appear fused in [`inds`](@ref), which carry +the dual. Compare with `inputnames`, which returns just the names. A plain tensor is a +trivial operator with no pairing, so its input indices are empty. + +# Examples + +```jldoctest +julia> op = operator(zeros(2, 2), ("i",), ("j",)); + +julia> inputinds(op) +1-element Vector{NamedUnitRange{String, Int64, Base.OneTo{Int64}}}: + named(Base.OneTo(2), "j") +``` + +See also [`inputnames`](@ref), [`outputinds`](@ref), [`inputaxes`](@ref), [`operator`](@ref). +""" +inputinds(a::AbstractNamedTensor) = conj.(inds(a)[dims(a, inputnames(a))]) + +""" + inputaxes(a) + +The input (domain) indices of an operator `a` as a `Tuple`, the tuple form of +[`inputinds`](@ref) (mirroring how `axes` relates to [`inds`](@ref)). Like `inputinds` it is +the non-dual domain space, so `a * randn(inputaxes(a))` contracts. + +See also [`inputinds`](@ref), [`outputaxes`](@ref), [`operator`](@ref). +""" +function inputaxes(a::AbstractNamedTensor) + return map(conj ∘ Base.Fix1(inds, a), Tuple(dims(a, inputnames(a)))) +end + # `outputname(a, i, default)` returns the output name paired with input name `i`, and # `inputname(a, i, default)` returns the input name paired with output name `i`, each # returning `default` when `i` is not one of the operator's paired names (mirroring diff --git a/test/test_exports.jl b/test/test_exports.jl index e575ae45..f0c5e100 100644 --- a/test/test_exports.jl +++ b/test/test_exports.jl @@ -6,8 +6,9 @@ using Test: @test, @testset :Index, :NamedUnitRange, :aligndims, :aligneddims, :apply, :commonind, :commoninds, :dimnames, :dimnametype, :hascommoninds, :id, - :inds, :inputnames, :mapinds, :named, :nameddims, :noncommonind, :noncommoninds, - :noprime, :operator, :outputnames, + :inds, :inputaxes, :inputinds, :inputnames, :mapinds, :named, :nameddims, + :noncommonind, :noncommoninds, + :noprime, :operator, :outputaxes, :outputinds, :outputnames, :prime, :replaceinds, :sim, :similar_operator, :state, :trycommonind, :trynoncommonind, :tryuniqueind, :uniqueind, :uniqueinds, :unioninds, :uniquename, diff --git a/test/test_operator.jl b/test/test_operator.jl index 8a98b2a4..e499a0bf 100644 --- a/test/test_operator.jl +++ b/test/test_operator.jl @@ -1,9 +1,11 @@ -using ITensorBase: ITensorBase as NDA, NamedTensor, NamedTensorOperator, apply, dimnames, - id, inputname, inputnames, nameddims, namedoneto, operator, outputname, outputnames, - product, replacedimnames, similar_operator, state, unname, unnamed +using GradedArrays: U1, gradedrange, isdual +using ITensorBase: ITensorBase as NDA, Index, NamedTensor, NamedTensorOperator, apply, + dimnames, id, inds, inputaxes, inputinds, inputname, inputnames, nameddims, namedoneto, + operator, outputaxes, outputinds, outputname, outputnames, product, replacedimnames, + similar_operator, state, unname, unnamed using LinearAlgebra: I, norm using MatrixAlgebraKit: project_hermitian -using Random: Random +using Random: Random, randn using StableRNGs: StableRNG using TensorAlgebra.MatrixAlgebra: gram_eigh_full, gram_eigh_full_with_pinv, invsqrth_safe, sqrth_invsqrth_safe, sqrth_safe @@ -61,6 +63,45 @@ end @test_throws ArgumentError inputname(NamedTensor(randn(2), ("i",)), "i") end +@testset "output/input indices" begin + o = operator(randn(2, 3, 2, 3), ("i'", "j'"), ("i", "j")) + # The index versions carry the paired names, in name order. Output indices are the + # codomain legs as-is; input indices are the domain, conjugated off the fused `inds`. + @test NDA.name.(outputinds(o)) == outputnames(o) + @test NDA.name.(inputinds(o)) == inputnames(o) + @test outputinds(o) == inds(o)[1:2] + @test inputinds(o) == conj.(inds(o)[3:4]) + # The `axes` forms are the tuple versions of the `inds` forms. + @test outputaxes(o) == Tuple(outputinds(o)) + @test inputaxes(o) == Tuple(inputinds(o)) + # A plain tensor is a trivial operator with no pairing, so all are empty. + v = NamedTensor(randn(2, 3), ("i", "j")) + @test isempty(outputinds(v)) + @test isempty(inputinds(v)) + @test outputaxes(v) == () + @test inputaxes(v) == () +end + +@testset "input indices carry the operator's domain space" begin + # `id` needs matching codomain/domain structure, so `i` and `j` share a graded range. + i = Index(gradedrange([U1(0) => 1, U1(1) => 1])) + j = Index(gradedrange([U1(0) => 1, U1(1) => 1])) + o = operator(id(Float64, (i,), (j,)), (i,), (j,)) + # `inds(o)` is the fused view: the domain (input) leg carries the dual. + @test isdual(inds(o)[2]) + # `outputinds`/`inputinds` are the codomain/domain spaces, both non-dual (the + # spaces states occupy), matching TensorKit's `codomain`/`domain`. + @test !isdual(only(outputinds(o))) + @test !isdual(only(inputinds(o))) + # The `axes` forms are the tuple versions, feeding constructors that take index tuples. + @test outputaxes(o) == Tuple(outputinds(o)) + @test inputaxes(o) == Tuple(inputinds(o)) + # A state built over the input space contracts, landing in the output space. + @test only(inds(o * randn(inputaxes(o)))) == only(outputinds(o)) + # `outputaxes`/`inputaxes` reconstruct an operator of the same shape. + @test isdual.(inds(id(Float64, outputaxes(o), inputaxes(o)))) == isdual.(inds(o)) +end + @testset "apply composition" begin # Operator applied to another operator: `x` lands on `y`'s output, and the result is an # operator on `x`'s input space keeping `y`'s input.