Add ITensorNetworkOperator and QuadraticFormNetwork data structures

src/ITensorNetworksNext.jl

@@ -13,8 +13,10 @@ include("select_algorithm.jl")
 include("AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl")
 include("abstracttensornetwork.jl")
 include("tensornetwork.jl")
-include("normnetwork.jl")
-include("normnetworkview.jl")
+include("itensornetworkoperator.jl")
+include("bilinearforms/abstractbilinearformnetwork.jl")
+include("bilinearforms/normnetwork.jl")
+include("bilinearforms/quadraticformnetwork.jl")
 include("ITensorNetworkGenerators/ITensorNetworkGenerators.jl")
 include("contract_network.jl")
 

src/abstracttensornetwork.jl

@@ -133,9 +133,15 @@ function insertlink!(tn::AbstractGraph, e)
 end
 
 function operator_support(tn::AbstractGraph, op::ITensorOperator)
+    return operator_support_names(tn, inputnames(op))
+end
+
+# Shared with the `ITensorNetworkOperator` method, which is defined alongside that type
+# because it is not yet known at this point in the include order.
+function operator_support_names(tn::AbstractGraph, opnames)
     support = Indices{vertextype(tn)}()
 
-    for name in inputnames(op)
+    for name in opnames
         vertices = dimnamevertices(tn, name)
 
         if length(vertices) > 1

src/beliefpropagation/beliefpropagation.jl

@@ -237,24 +237,13 @@ end
     contraction_alg::ContractionAlg = Exact()
 end
 
-# The tensors making up the factor at `vertex`, as separate operands for `contract_network`. A
-# `NormNetwork`'s factor is a lazy `ket * conj(bra)` product, and the contraction order sees each
-# operand as one node carrying only its outer axes — which hides the physical index the two layers
-# share, forcing the doubled vertex to be formed before any message is absorbed (χ^(2 * degree)
-# rather than the χ^(degree + 1) an interleaved order reaches).
-factor_tensors(factors, vertex) = [factors[vertex]]
-function factor_tensors(factors::NormNetwork, vertex)
-    return [kettensor(factors, vertex), bratensor(factors, vertex)]
-end
-
 # Contract the incoming messages into the source factor to form the (unnormalized) new message on
 # `edge`.
 function updated_message(algorithm::SimpleMessageUpdate, cache, factors, edge)
     messages = collect(incoming_messages(cache, edge))
-    # TODO: Remove `factor_tensors` once `contract_network` handles lazy tensors in
-    # contraction sequences properly.
     return contract_network(
-        [messages; factor_tensors(factors, src(edge))]; alg = algorithm.contraction_alg
+        [messages; [factors[src(edge)]]];
+        alg = algorithm.contraction_alg
     )
 end
 
@@ -275,9 +264,9 @@ end
 # the message is positive semidefinite and its trace is a positive normalization.
 function message_update!(algorithm::SimpleMessageUpdate, cache, factors::NormNetwork, edge)
     new_tensor = updated_message(algorithm, cache, factors, edge)
-    new_message = operator(
-        new_tensor, linknames(BraView(factors), edge), linknames(KetView(factors), edge)
-    )
+    branames = linknames(branetwork(factors), edge)
+    ketnames = linknames(ketnetwork(factors), edge)
+    new_message = operator(new_tensor, branames, ketnames)
     if algorithm.normalize
         message_norm = tr(new_message)
         iszero(message_norm) || (new_message /= message_norm)

src/beliefpropagation/messagecache.jl

@@ -132,10 +132,7 @@ end
 
 function vertex_scalar(factors, messages, vertex; kwargs...)
     in_messages = incoming_edge_data(messages, [vertex])
-    # TODO: Remove `factor_tensors` once `contract_network` handles lazy tensors in
-    # contraction sequences properly.
-    tensors = [factor_tensors(factors, vertex); collect(in_messages)]
-    return contract_network(tensors; kwargs...)[]
+    return contract_network([[factors[vertex]]; collect(in_messages)]; kwargs...)[]
 end
 
 vertex_scalars(factors, messages) = vertex_scalars(factors, messages, keys(factors))
@@ -193,16 +190,16 @@ bethe_free_energy(factors, messages) = -bethe_free_entropy(factors, messages)
 function similar_message_environment(nn::NormNetwork)
     messages = mapmany(vertices(nn)) do vertex
         return map(in_incident_edges(nn, vertex)) do edge
-            braview = BraView(nn)
-            ketview = KetView(nn)
+            bra = branetwork(nn)
+            ket = ketnetwork(nn)
 
-            ketnames = linknames(ketview, edge)
-            branames = linknames(braview, edge)
-            braaxis = unnamed.(linkaxes(braview, edge))
+            ketnames = linknames(ket, edge)
+            branames = linknames(bra, edge)
+            braaxis = unnamed.(linkaxes(bra, edge))
 
             # Bra leg = operator output, ket leg = input, the bipartition in which the message
             # is positive semidefinite.
-            message = similar_operator(ketview[vertex], braaxis, branames, ketnames)
+            message = similar_operator(ket[vertex], braaxis, branames, ketnames)
 
             return edge => message
         end

src/bilinearforms/abstractbilinearformnetwork.jl

@@ -0,0 +1,194 @@
+using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
+using Dictionaries: Dictionaries, Dictionary, isinsertable, issettable
+using Graphs: Graphs, edges, vertices
+using ITensorBase: ITensorBase, conj, inds, name, rename
+using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
+
+"""
+    abstract type AbstractBilinearFormNetwork{T, V, I} <: AbstractITensorNetwork{T, V}
+
+Supertype of the lazy multi-layer networks built from a ket layer of type
+`ITensorNetwork{T, V, I}` and a ket→bra index name mapping.
+
+A subtype supplies its own graph structure, implements [`braname`](@ref), and returns an
+[`AbstractGramian`](@ref) from `getindex`; [`kettensor`](@ref), [`bratensor`](@ref) and,
+where the subtype has an operator layer, [`operatortensor`](@ref) read a vertex's layers from
+that Gramian. The layers as whole networks are returned by [`ketnetwork`](@ref),
+[`branetwork`](@ref) and [`operatornetwork`](@ref).
+"""
+abstract type AbstractBilinearFormNetwork{T, V, I} <: AbstractITensorNetwork{T, V} end
+
+"""
+    abstract type AbstractGramian
+
+The layers of an `AbstractBilinearFormNetwork` at one vertex. A subtype implements
+[`kettensor`](@ref), [`braname`](@ref), `layertensors` and `layerinds`; the bra tensor is built
+from the ket tensor and the name map each time it is requested. Its `inds`, `names` and `axes`
+are the indices of its layers that no other layer shares, those the layer product leaves open.
+"""
+abstract type AbstractGramian end
+
+# ====================================== Graphs.jl ======================================= #
+
+Graphs.edges(bn::AbstractBilinearFormNetwork) = edges(ketnetwork(bn))
+Graphs.vertices(bn::AbstractBilinearFormNetwork) = vertices(ketnetwork(bn))
+
+# ==================================== NamedGraphs.jl ==================================== #
+
+function NamedGraphs.encoded_vertex(bn::AbstractBilinearFormNetwork, vertex)
+    return encoded_vertex(ketnetwork(bn), vertex)
+end
+function NamedGraphs.decoded_vertex(bn::AbstractBilinearFormNetwork, code::Integer)
+    return decoded_vertex(ketnetwork(bn), code)
+end
+NamedGraphs.encoded_graph(bn::AbstractBilinearFormNetwork) = encoded_graph(ketnetwork(bn))
+
+# ==================================== DataGraphs.jl ===================================== #
+
+function DataGraphs.is_vertex_assigned(bn::AbstractBilinearFormNetwork, vertex)
+    return isassigned(ketnetwork(bn), vertex)
+end
+
+# =================================== Dictionaries.jl ==================================== #
+
+Dictionaries.issettable(::AbstractBilinearFormNetwork) = false
+Dictionaries.isinsertable(::AbstractBilinearFormNetwork) = false
+
+# ====================================== interface ======================================= #
+
+"""
+    braname(bn::AbstractBilinearFormNetwork, name)
+    braname(g::AbstractGramian, name)
+
+The bra-layer index name corresponding to the ket-layer index name `name`. The `AbstractGramian`
+form maps a name absent from its name map to itself, without checking it belongs to the network.
+"""
+function braname end
+function braname(bn::AbstractBilinearFormNetwork, name)
+    if !has_dimname(ketnetwork(bn), name)
+        error("index name $name not found underlying tensor network.")
+    end
+    # A name absent from the map has no separate bra copy and maps to itself: a site index of a
+    # norm network, or a site index a quadratic form's operator does not act on.
+    return get(branamemap(bn), name, name)
+end
+braname(g::AbstractGramian, name) = get(branamemap(g), name, name)
+
+"""
+    branamemap(bn::AbstractBilinearFormNetwork)
+    branamemap(g::AbstractGramian)
+
+The ket→bra name map, holding a bra name for each ket index name that has a separate bra copy.
+"""
+function branamemap end
+
+# A link name, or a name in `acted`, gets its bra name from `map`; every other name has none.
+function select_branames(ket::ITensorNetwork{T, V, I}, map, acted) where {T, V, I}
+    braname = Dictionary{I, I}()
+    for (name, vertices) in pairs(ket.dimname_vertices)
+        if length(vertices) == 2 || name in acted
+            insert!(braname, name, map[name])
+        end
+    end
+    return braname
+end
+
+"""
+    kettensor(g::AbstractGramian)
+
+The ket-layer tensor of the Gramian `g`.
+"""
+function kettensor end
+
+"""
+    operatortensor(g::AbstractGramian)
+
+The operator-layer tensor of the Gramian `g`, with its index names renamed so that its input
+legs meet the ket layer and its output legs meet the bra layer.
+"""
+function operatortensor end
+
+conj_bratensor(g::AbstractGramian) = rename(n -> braname(g, n), kettensor(g))
+
+"""
+    bratensor(g::AbstractGramian)
+
+The bra-layer tensor of the Gramian `g`.
+"""
+bratensor(g::AbstractGramian) = conj(conj_bratensor(g))
+
+# Read from `conj_bratensor`, which only renames, so the tensor data is not conjugated.
+brainds(g::AbstractGramian) = conj.(inds(conj_bratensor(g)))
+
+function ITensorBase.inds(g::AbstractGramian)
+    layer_inds = reduce(vcat, collect.(layerinds(g)))
+    layer_names = name.(layer_inds)
+    return [i for i in layer_inds if count(==(name(i)), layer_names) == 1]
+end
+ITensorBase.names(g::AbstractGramian) = name.(inds(g))
+Base.axes(g::AbstractGramian) = Tuple(inds(g))
+
+"""
+    ketnetwork(bn::AbstractBilinearFormNetwork)
+
+The ket-layer network of `bn`.
+"""
+function ketnetwork end
+
+"""
+    operatornetwork(bn::AbstractBilinearFormNetwork)
+
+The operator-layer network of `bn`, for a subtype that has an operator layer.
+"""
+function operatornetwork end
+
+"""
+    branetwork(bn::AbstractBilinearFormNetwork)
+
+The bra-layer network of `bn`. Unless a subtype stores its bra layer as a network, this is a
+`BraView`, whose tensors are built by [`bratensor`](@ref) when accessed.
+"""
+branetwork(bn::AbstractBilinearFormNetwork) = BraView(bn)
+
+"""
+    struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <: AbstractITensorNetwork{T, V}
+
+The bra layer of the bilinear-form network `parent(view)`, with each vertex tensor built by
+[`bratensor`](@ref) when accessed. Its graph structure and mutability are those of the parent.
+"""
+struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
+    AbstractITensorNetwork{T, V}
+    parent::P
+    function BraView(parent::AbstractBilinearFormNetwork{T, V, I}) where {T, V, I}
+        return new{T, V, I, typeof(parent)}(parent)
+    end
+end
+
+Base.parent(nnv::BraView) = nnv.parent
+
+# ==================================== DataGraphs.jl ===================================== #
+
+DataGraphs.get_vertex_data(nnv::BraView, vertex) = bratensor(parent(nnv)[vertex])
+function DataGraphs.is_vertex_assigned(nnv::BraView, vertex)
+    return is_vertex_assigned(parent(nnv), vertex)
+end
+
+# ====================================== Graphs.jl ======================================= #
+
+Graphs.edges(nnv::BraView) = edges(parent(nnv))
+Graphs.vertices(nnv::BraView) = vertices(parent(nnv))
+
+# ==================================== NamedGraphs.jl ==================================== #
+
+function NamedGraphs.encoded_vertex(nnv::BraView, vertex)
+    return encoded_vertex(parent(nnv), vertex)
+end
+function NamedGraphs.decoded_vertex(nnv::BraView, code::Integer)
+    return decoded_vertex(parent(nnv), code)
+end
+NamedGraphs.encoded_graph(nnv::BraView) = encoded_graph(parent(nnv))
+
+# =================================== Dictionaries.jl ==================================== #
+
+Dictionaries.issettable(nnv::BraView) = issettable(parent(nnv))
+Dictionaries.isinsertable(nnv::BraView) = isinsertable(parent(nnv))

src/bilinearforms/normnetwork.jl

@@ -0,0 +1,65 @@
+using Dictionaries: Dictionary
+using ITensorBase: similar_operator, uniquename
+using ITensorNetworksNext
+
+"""
+    struct NormNetwork{T, V, I} <: AbstractBilinearFormNetwork{T, V, I}
+
+Lazy wrapper representing the norm `⟨tn|tn⟩` of `tn::ITensorNetwork{T, V, I}`,
+together with a per-edge ket→bra name mapping that, for each index in the ket layer, defines
+the name of the corresponding index in the bra layer.
+"""
+struct NormNetwork{T, V, I} <: AbstractBilinearFormNetwork{T, V, I}
+    ket::ITensorNetwork{T, V, I}
+    braname::Dictionary{I, I}
+    function NormNetwork(
+            ket::ITensorNetwork{T, V, I},
+            map::Dictionary{I, I}
+        ) where {T, V, I}
+        return new{T, V, I}(ket, select_branames(ket, map, ()))
+    end
+end
+
+"""
+    struct NormGramian{T, I} <: AbstractGramian
+
+The layers of a `NormNetwork` at one vertex: the ket tensor and the network's ket→bra name map,
+from which the bra tensor is built when requested.
+"""
+struct NormGramian{T, I} <: AbstractGramian
+    ket::T
+    braname::Dictionary{I, I}
+end
+
+kettensor(g::NormGramian) = g.ket
+branamemap(g::NormGramian) = g.braname
+layertensors(g::NormGramian) = (; ket = kettensor(g), bra = bratensor(g))
+layerinds(g::NormGramian) = (inds(kettensor(g)), brainds(g))
+
+Base.eltype(::Type{<:NormNetwork{T, V, I}}) where {T, V, I} = NormGramian{T, I}
+
+function NormNetwork(tn::ITensorNetwork)
+    return NormNetwork(tn, map(uniquename, keys(tn.dimname_vertices)))
+end
+
+# ==================================== DataGraphs.jl ===================================== #
+
+function DataGraphs.get_vertex_data(nn::NormNetwork{T, V, I}, vertex) where {T, V, I}
+    return NormGramian{T, I}(nn.ket[vertex], nn.braname)
+end
+
+# ====================================== interface ======================================= #
+
+ketnetwork(nn::NormNetwork) = nn.ket
+branamemap(nn::NormNetwork) = nn.braname
+
+"""
+    normnetwork(tn::ITensorNetwork, [braname]) -> NormNetwork
+
+Build the double-layer norm network `⟨tn|tn⟩`, represented lazily as a `NomnNetwork` object.
+The optional second argument `braname` should implement `braname[ketdimname] = bradimname` for
+every link dimension name `ketdimname` in `tn`. If this is not specified, then a name is
+generated via the `ITensorBase.uniquename` function.
+"""
+normnetwork(tn::ITensorNetwork) = NormNetwork(tn)
+normnetwork(tn::ITensorNetwork, braname) = NormNetwork(tn, braname)