From 7150d5c9eebe768852a8f374b6f578ca43c800f1 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Wed, 9 Sep 2026 11:02:33 -0400
Subject: [PATCH 01/18] Add `QuadraticFormNetwork` and a shared bilinear-form
 interface
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`QuadraticFormNetwork` represents `⟨tn|op|tn⟩` lazily, with an operator layer
between the bra and ket layers. A ket index name is renamed in the bra layer
when it is a link name or appears in the operator's `inputnames`, so a ket index
the operator does not act on still contracts straight onto the bra.
`operatortensor` renames each operator output name to the bra name of the input
name it is paired with, leaving the input names to meet the ket. The constructor
throws unless `vertices(operator)` and `vertices(ket)` are set-equal.

`AbstractBilinearFormNetwork` is the new supertype of `NormNetwork` and
`QuadraticFormNetwork`. A subtype implements `braname`, `kettensor` and, where it
has an operator layer, `operatortensor`; `bratensor`, `conj_bratensor` and
`indmap` are defined once in terms of those, so the two concrete files now share
them instead of each carrying a copy. `AbstractBilinearFormNetworkView` does the
same for `KetView`, `BraView` and the new `OperatorView`: each implements
`Base.parent` and `get_vertex_data`, and takes its graph structure, assignment
and mutability from the parent network.

`test/test_quadraticformnetwork.jl` covers the name map, the three views, and
three contractions: an identity operator layer reproduces
`contract_network(NormNetwork(tn))`, a random on-site product operator matches a
dense `⟨ψ|O|ψ⟩` reference, and a `2I` layer on three vertices scales the norm by
8.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/ITensorNetworksNext.jl             |   3 +
 src/abstractbilinearformnetwork.jl     |  57 +++++++
 src/abstractbilinearformnetworkview.jl |  43 ++++++
 src/normnetwork.jl                     |  25 +---
 src/normnetworkview.jl                 |  55 ++++---
 src/quadraticformnetwork.jl            | 110 ++++++++++++++
 test/test_quadraticformnetwork.jl      | 199 +++++++++++++++++++++++++
 7 files changed, 444 insertions(+), 48 deletions(-)
 create mode 100644 src/abstractbilinearformnetwork.jl
 create mode 100644 src/abstractbilinearformnetworkview.jl
 create mode 100644 src/quadraticformnetwork.jl
 create mode 100644 test/test_quadraticformnetwork.jl

diff --git a/src/ITensorNetworksNext.jl b/src/ITensorNetworksNext.jl
index fe1b2540..6a24040f 100644
--- a/src/ITensorNetworksNext.jl
+++ b/src/ITensorNetworksNext.jl
@@ -13,7 +13,10 @@ include("select_algorithm.jl")
 include("AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl")
 include("abstracttensornetwork.jl")
 include("tensornetwork.jl")
+include("abstractbilinearformnetwork.jl")
+include("abstractbilinearformnetworkview.jl")
 include("normnetwork.jl")
+include("quadraticformnetwork.jl")
 include("normnetworkview.jl")
 include("ITensorNetworkGenerators/ITensorNetworkGenerators.jl")
 include("contract_network.jl")
diff --git a/src/abstractbilinearformnetwork.jl b/src/abstractbilinearformnetwork.jl
new file mode 100644
index 00000000..6d1e894e
--- /dev/null
+++ b/src/abstractbilinearformnetwork.jl
@@ -0,0 +1,57 @@
+using Dictionaries: Dictionaries
+using ITensorBase: conj, name, rename, setname
+
+"""
+    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 and implements [`braname`](@ref) together with one
+accessor per layer: [`kettensor`](@ref), [`bratensor`](@ref) and, where the subtype has an
+operator layer, [`operatortensor`](@ref). `bratensor` has a default built from `kettensor`
+and `braname`.
+"""
+abstract type AbstractBilinearFormNetwork{T, V, I} <: AbstractITensorNetwork{T, V} end
+
+# =================================== Dictionaries.jl ==================================== #
+
+Dictionaries.issettable(::AbstractBilinearFormNetwork) = false
+Dictionaries.isinsertable(::AbstractBilinearFormNetwork) = false
+
+# ====================================== interface ======================================= #
+
+"""
+    braname(bn::AbstractBilinearFormNetwork, name)
+
+The bra-layer index name corresponding to the ket-layer index name `name`.
+"""
+function braname end
+
+"""
+    kettensor(bn::AbstractBilinearFormNetwork, vertex)
+
+The ket-layer tensor at `vertex`.
+"""
+function kettensor end
+
+"""
+    operatortensor(bn::AbstractBilinearFormNetwork, vertex)
+
+The operator-layer tensor at `vertex`, 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
+
+function conj_bratensor(bn::AbstractBilinearFormNetwork, vertex)
+    return rename(n -> braname(bn, n), kettensor(bn, vertex))
+end
+
+"""
+    bratensor(bn::AbstractBilinearFormNetwork, vertex)
+
+The bra-layer tensor at `vertex`.
+"""
+bratensor(bn::AbstractBilinearFormNetwork, vertex) = conj(conj_bratensor(bn, vertex))
+
+indmap(bn::AbstractBilinearFormNetwork, ind) = setname(conj(ind), braname(bn, name(ind)))
diff --git a/src/abstractbilinearformnetworkview.jl b/src/abstractbilinearformnetworkview.jl
new file mode 100644
index 00000000..828cbce1
--- /dev/null
+++ b/src/abstractbilinearformnetworkview.jl
@@ -0,0 +1,43 @@
+using DataGraphs: DataGraphs, is_vertex_assigned
+using Dictionaries: Dictionaries, isinsertable, issettable
+using Graphs: Graphs, edges, vertices
+using NamedGraphs: NamedGraphs, ordered_vertices, position_graph, vertex_positions
+
+"""
+    abstract type AbstractBilinearFormNetworkView{T, V, I} <: AbstractITensorNetwork{T, V}
+
+Supertype of the single-layer views of an `AbstractBilinearFormNetwork{T, V, I}`.
+
+A subtype implements `Base.parent`, returning the network it views, and
+`DataGraphs.get_vertex_data`, returning that layer's tensor at a vertex. Its graph structure
+and mutability are those of the parent network.
+"""
+abstract type AbstractBilinearFormNetworkView{T, V, I} <: AbstractITensorNetwork{T, V} end
+
+# ====================================== Graphs.jl ======================================= #
+
+Graphs.edges(nnv::AbstractBilinearFormNetworkView) = edges(parent(nnv))
+Graphs.vertices(nnv::AbstractBilinearFormNetworkView) = vertices(parent(nnv))
+
+# ==================================== NamedGraphs.jl ==================================== #
+
+function NamedGraphs.vertex_positions(nnv::AbstractBilinearFormNetworkView)
+    return vertex_positions(parent(nnv))
+end
+function NamedGraphs.ordered_vertices(nnv::AbstractBilinearFormNetworkView)
+    return ordered_vertices(parent(nnv))
+end
+function NamedGraphs.position_graph(nnv::AbstractBilinearFormNetworkView)
+    return position_graph(parent(nnv))
+end
+
+# ==================================== DataGraphs.jl ===================================== #
+
+function DataGraphs.is_vertex_assigned(nnv::AbstractBilinearFormNetworkView, vertex)
+    return is_vertex_assigned(parent(nnv), vertex)
+end
+
+# =================================== Dictionaries.jl ==================================== #
+
+Dictionaries.issettable(nnv::AbstractBilinearFormNetworkView) = issettable(parent(nnv))
+Dictionaries.isinsertable(nnv::AbstractBilinearFormNetworkView) = isinsertable(parent(nnv))
diff --git a/src/normnetwork.jl b/src/normnetwork.jl
index 46725507..6428644f 100644
--- a/src/normnetwork.jl
+++ b/src/normnetwork.jl
@@ -1,15 +1,15 @@
 using Dictionaries: Dictionary
-using ITensorBase: LazyNamedTensor, lazy, rename, setname, similar_operator, uniquename
+using ITensorBase: LazyNamedTensor, lazy, similar_operator, uniquename
 using ITensorNetworksNext
 
 """
-    struct NormNetwork{T, V, I} <: AbstractITensorNetwork{T, V}
+    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} <: AbstractITensorNetwork{T, V}
+struct NormNetwork{T, V, I} <: AbstractBilinearFormNetwork{T, V, I}
     ket::ITensorNetwork{T, V, I}
     braname::Dictionary{I, I}
     function NormNetwork(
@@ -45,6 +45,10 @@ NamedGraphs.encoded_graph(nn::NormNetwork) = encoded_graph(nn.ket)
 
 # ==================================== DataGraphs.jl ===================================== #
 
+function DataGraphs.is_vertex_assigned(nn::NormNetwork, vertex)
+    return isassigned(nn.ket, vertex)
+end
+
 function DataGraphs.get_vertex_data(nn::NormNetwork, vertex)
     A = kettensor(nn, vertex)
     B = conj_bratensor(nn, vertex)
@@ -52,14 +56,6 @@ function DataGraphs.get_vertex_data(nn::NormNetwork, vertex)
     return lazy(A) * lazy(conj(B))
 end
 
-function DataGraphs.is_vertex_assigned(nn::NormNetwork, vertex)
-    return isassigned(nn.ket, vertex)
-end
-# =================================== Dictionaries.jl ==================================== #
-
-Dictionaries.issettable(::NormNetwork) = false
-Dictionaries.isinsertable(::NormNetwork) = false
-
 # ====================================== interface ======================================= #
 
 function braname(nn::NormNetwork, name)
@@ -71,14 +67,7 @@ function braname(nn::NormNetwork, name)
     return get(nn.braname, name, name)
 end
 
-indmap(nn::NormNetwork, ind) = setname(conj(ind), braname(nn, name(ind)))
-
 kettensor(nn::NormNetwork, vertex) = nn.ket[vertex]
-function conj_bratensor(nn::NormNetwork, vertex)
-    return rename(n -> braname(nn, n), kettensor(nn, vertex))
-end
-
-bratensor(nn::NormNetwork, vertex) = conj(conj_bratensor(nn, vertex))
 
 """
     normnetwork(tn::ITensorNetwork, [braname]) -> NormNetwork
diff --git a/src/normnetworkview.jl b/src/normnetworkview.jl
index 75c57926..3b197aeb 100644
--- a/src/normnetworkview.jl
+++ b/src/normnetworkview.jl
@@ -3,20 +3,11 @@ using Dictionaries: Dictionaries, isinsertable, issettable
 using Graphs: Graphs, edges, vertices
 using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
 
-struct KetView{T, V, I} <: AbstractITensorNetwork{T, V}
-    parent::NormNetwork{T, V, I}
-end
-
-struct BraView{T, V, I} <: AbstractITensorNetwork{T, V}
-    parent::NormNetwork{T, V, I}
-end
-
-# ====================================== Graphs.jl ======================================= #
-
-for View in (:KetView, :BraView)
-    @eval begin
-        Graphs.edges(nnv::$View) = edges(nnv.parent)
-        Graphs.vertices(nnv::$View) = vertices(nnv.parent)
+struct KetView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
+    AbstractBilinearFormNetworkView{T, V, I}
+    parent::P
+    function KetView(parent::AbstractBilinearFormNetwork{T, V, I}) where {T, V, I}
+        return new{T, V, I, typeof(parent)}(parent)
     end
 end
 
@@ -35,24 +26,28 @@ for View in (:KetView, :BraView)
     end
 end
 
-# ==================================== DataGraphs.jl ===================================== #
-
-DataGraphs.get_vertex_data(nn::KetView, vertex) = kettensor(nn.parent, vertex)
-DataGraphs.get_vertex_data(nn::BraView, vertex) = bratensor(nn.parent, vertex)
-
-for View in (:KetView, :BraView)
-    @eval begin
-        function DataGraphs.is_vertex_assigned(nnv::$View, vertex)
-            return isassigned(nnv.parent.ket, vertex)
-        end
+struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
+    AbstractBilinearFormNetworkView{T, V, I}
+    parent::P
+    function BraView(parent::AbstractBilinearFormNetwork{T, V, I}) where {T, V, I}
+        return new{T, V, I, typeof(parent)}(parent)
     end
 end
 
-# =================================== Dictionaries.jl ==================================== #
-
-for View in (:KetView, :BraView)
-    @eval begin
-        Dictionaries.issettable(nnv::$View) = issettable(nnv.parent)
-        Dictionaries.isinsertable(nnv::$View) = isinsertable(nnv.parent)
+struct OperatorView{T, V, I, P <: QuadraticFormNetwork{T, V, I}} <:
+    AbstractBilinearFormNetworkView{T, V, I}
+    parent::P
+    function OperatorView(parent::QuadraticFormNetwork{T, V, I}) where {T, V, I}
+        return new{T, V, I, typeof(parent)}(parent)
     end
 end
+
+Base.parent(nnv::KetView) = nnv.parent
+Base.parent(nnv::BraView) = nnv.parent
+Base.parent(nnv::OperatorView) = nnv.parent
+
+# ==================================== DataGraphs.jl ===================================== #
+
+DataGraphs.get_vertex_data(nnv::KetView, vertex) = kettensor(parent(nnv), vertex)
+DataGraphs.get_vertex_data(nnv::BraView, vertex) = bratensor(parent(nnv), vertex)
+DataGraphs.get_vertex_data(nnv::OperatorView, vertex) = operatortensor(parent(nnv), vertex)
diff --git a/src/quadraticformnetwork.jl b/src/quadraticformnetwork.jl
new file mode 100644
index 00000000..1b803476
--- /dev/null
+++ b/src/quadraticformnetwork.jl
@@ -0,0 +1,110 @@
+using Dictionaries: Dictionary
+using ITensorBase:
+    LazyNamedTensor, inputnames, lazy, outputnames, rename, state, uniquename
+using ITensorNetworksNext
+
+"""
+    struct QuadraticFormNetwork{T, V, I, O} <: AbstractBilinearFormNetwork{T, V, I}
+
+Lazy wrapper representing the quadratic form `⟨tn|op|tn⟩` of `tn::ITensorNetwork{T, V, I}`
+sandwiched around the operator layer `op::ITensorNetwork{O, V, I}`, together with a per-index
+ket→bra name mapping that, for each index in the ket layer, defines the name of the
+corresponding index in the bra layer.
+
+The operator layer must be defined on every vertex of `tn`. At each vertex the operator's
+input names are ket index names, and each paired output name is renamed to the bra layer, so
+the operator's remaining index names must be distinct from every index name of `tn`.
+"""
+struct QuadraticFormNetwork{T, V, I, O} <: AbstractBilinearFormNetwork{T, V, I}
+    ket::ITensorNetwork{T, V, I}
+    operator::ITensorNetwork{O, V, I}
+    braname::Dictionary{I, I}
+    function QuadraticFormNetwork(
+            ket::ITensorNetwork{T, V, I},
+            operator::ITensorNetwork{O, V, I},
+            map::Dictionary{I, I}
+        ) where {T, V, I, O}
+        if !issetequal(vertices(operator), vertices(ket))
+            error("the operator layer must be defined on every vertex of the ket layer.")
+        end
+        acted = Set{I}()
+        for vertex in vertices(operator)
+            union!(acted, inputnames(operator[vertex]))
+        end
+        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 new{T, V, I, O}(ket, operator, braname)
+    end
+end
+
+Base.eltype(::Type{<:QuadraticFormNetwork{T, V, I}}) where {T, V, I} = LazyNamedTensor{I, T}
+
+function QuadraticFormNetwork(ket::ITensorNetwork, operator::ITensorNetwork)
+    return QuadraticFormNetwork(ket, operator, map(uniquename, keys(ket.dimname_vertices)))
+end
+
+# ====================================== Graphs.jl ======================================= #
+
+Graphs.edges(qf::QuadraticFormNetwork) = edges(qf.ket)
+Graphs.vertices(qf::QuadraticFormNetwork) = vertices(qf.ket)
+
+# ==================================== NamedGraphs.jl ==================================== #
+
+NamedGraphs.vertex_positions(qf::QuadraticFormNetwork) = vertex_positions(qf.ket)
+NamedGraphs.ordered_vertices(qf::QuadraticFormNetwork) = ordered_vertices(qf.ket)
+NamedGraphs.position_graph(qf::QuadraticFormNetwork) = position_graph(qf.ket)
+
+# ==================================== DataGraphs.jl ===================================== #
+
+function DataGraphs.get_vertex_data(qf::QuadraticFormNetwork, vertex)
+    A = kettensor(qf, vertex)
+    O = operatortensor(qf, vertex)
+    B = conj_bratensor(qf, vertex)
+    # TODO: implement and use a lazy `conj` via `LazyNamedDimsArrays` here?
+    return lazy(A) * lazy(state(O)) * lazy(conj(B))
+end
+
+function DataGraphs.is_vertex_assigned(qf::QuadraticFormNetwork, vertex)
+    return isassigned(qf.ket, vertex) && isassigned(qf.operator, vertex)
+end
+
+# ====================================== interface ======================================= #
+
+function braname(qf::QuadraticFormNetwork, name)
+    if !has_dimname(qf.ket, name)
+        error("index name $name not found underlying tensor network.")
+    end
+    # The indices not stored in `qf.braname` are the dangling ket indices the operator does
+    # not act on, which get mapped to themselves.
+    return get(qf.braname, name, name)
+end
+
+kettensor(qf::QuadraticFormNetwork, vertex) = qf.ket[vertex]
+
+# Each output name is renamed to the bra name of the input name it is paired with, so the
+# operator's output legs meet the bra layer and its input legs meet the ket layer.
+function operatortensor(qf::QuadraticFormNetwork, vertex)
+    tensor = qf.operator[vertex]
+    names = map(inputnames(tensor)) do n
+        return braname(qf, n)
+    end
+    return rename(tensor, (outputnames(tensor) .=> names)...)
+end
+
+"""
+    quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetwork, [braname])
+
+Build the triple-layer network `⟨tn|op|tn⟩`, represented lazily as a
+`QuadraticFormNetwork` object. The optional third argument `braname` should implement
+`braname[ketdimname] = bradimname` for every link dimension name `ketdimname` in `tn` and
+every dimension name of `tn` the operator acts on. If this is not specified, then a name is
+generated via the `ITensorBase.uniquename` function.
+"""
+quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetwork) = QuadraticFormNetwork(tn, op)
+function quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetwork, braname)
+    return QuadraticFormNetwork(tn, op, braname)
+end
diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl
new file mode 100644
index 00000000..f303b259
--- /dev/null
+++ b/test/test_quadraticformnetwork.jl
@@ -0,0 +1,199 @@
+using DataGraphs: is_vertex_assigned
+using Dictionaries: isinsertable, issettable
+using Graphs: edges, vertices
+using ITensorBase: ITensor, Index, IndexName, conj, inds, inputnames, name, operator,
+    outputnames, rename, setname, uniquename
+using ITensorNetworksNext: BraView, ITensorNetwork, KetView, NormNetwork, OperatorView,
+    QuadraticFormNetwork, braname, bratensor, conj_bratensor, contract_network, indmap,
+    kettensor, operatortensor, quadraticformnetwork, tensornetwork
+using LinearAlgebra: I, norm
+using NamedGraphs.GraphsExtensions: incident_edges
+using NamedGraphs.NamedGraphGenerators: named_grid, named_path_graph
+using NamedGraphs: NamedEdge
+using Test: @test, @test_throws, @testset
+
+# Build a random `ITensorNetwork` state on the graph `g` with site dimension `d` and
+# bond dimension `χ`.
+function random_state(::Type{T}, g; d = 2, χ = 2) where {T}
+    l = Dict(e => Index(χ) for e in edges(g))
+    l = merge(l, Dict(reverse(e) => l[e] for e in edges(g)))
+    s = Dict(v => Index(d) for v in vertices(g))
+    tn = tensornetwork(vertices(g)) do v
+        is = map(e -> l[e], incident_edges(g, v))
+        return randn(T, (s[v], is...))
+    end
+    return tn, l, s
+end
+
+# Build a bondless operator layer on the vertices of `g`, with `f(v)` supplying the matrix
+# acting on the site index `s[v]`.
+function product_operator(f, g, s; d = 2)
+    op = tensornetwork(vertices(g)) do v
+        out = Index(d)
+        return operator(ITensor(f(v), (out, s[v])), (name(out),), (name(s[v]),))
+    end
+    return op
+end
+
+identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s; d)
+
+@testset "`QuadraticFormNetwork`" begin
+    @testset "Basics" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        op = identity_operator(g, s)
+        qf = QuadraticFormNetwork(tn, op)
+
+        # `quadraticformnetwork` is the public constructor and agrees with the type.
+        @test quadraticformnetwork(tn, op) isa QuadraticFormNetwork
+        @test qf isa QuadraticFormNetwork
+
+        # The quadratic form shares the graph structure of the ket layer.
+        @test issetequal(vertices(qf), vertices(tn))
+        @test issetequal(edges(qf), edges(tn))
+
+        # `eltype` is the type of the (lazy triple-layer) vertex data.
+        @test eltype(qf) === typeof(qf[1])
+
+        # Vertex data is assigned wherever both layers are.
+        @test is_vertex_assigned(qf, 1)
+
+        # The quadratic form is neither settable nor insertable (it is a lazy view).
+        @test !issettable(qf)
+        @test !isinsertable(qf)
+
+        # The operator layer must cover every vertex of the ket layer.
+        g4 = named_path_graph(4)
+        _, _, s4 = random_state(Float64, g4)
+        @test_throws ErrorException QuadraticFormNetwork(tn, identity_operator(g4, s4))
+    end
+
+    @testset "layer tensors and the name map" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        op = identity_operator(g, s)
+        qf = QuadraticFormNetwork(tn, op)
+
+        # `kettensor` returns the underlying tensor untouched.
+        @test kettensor(qf, 2) === tn[2]
+
+        # Unlike the norm network, the site indices *are* renamed in the bra layer: the
+        # operator sits between the two layers, so they no longer contract directly.
+        sname = name(s[2])
+        @test braname(qf, sname) != sname
+        @test sname in name.(inds(kettensor(qf, 2)))
+        @test !(sname in name.(inds(conj_bratensor(qf, 2))))
+        @test braname(qf, sname) in name.(inds(conj_bratensor(qf, 2)))
+
+        # Link indices are shared by two tensors, so they are renamed in the bra layer to
+        # keep the two layers' bonds distinct.
+        lname = name(l[NamedEdge(1 => 2)])
+        @test braname(qf, lname) != lname
+        @test lname in name.(inds(kettensor(qf, 2)))
+        @test !(lname in name.(inds(conj_bratensor(qf, 2))))
+        @test braname(qf, lname) in name.(inds(conj_bratensor(qf, 2)))
+
+        # The operator's input name meets the ket and its output name is renamed to meet
+        # the bra.
+        o = operatortensor(qf, 2)
+        @test inputnames(o) == [sname]
+        @test outputnames(o) == [braname(qf, sname)]
+
+        # `bratensor` is the elementwise conjugate of `conj_bratensor` and carries the same
+        # indices.
+        @test inds(bratensor(qf, 2)) == inds(conj_bratensor(qf, 2))
+
+        # `indmap` conjugates an index and renames it according to the name map.
+        ind = only(i for i in inds(kettensor(qf, 2)) if name(i) == lname)
+        @test name(indmap(qf, ind)) == braname(qf, name(ind))
+        @test indmap(qf, ind) == setname(conj(ind), braname(qf, name(ind)))
+
+        # Querying the name map with an index name absent from the ket layer errors.
+        @test_throws ErrorException braname(qf, name(Index(2)))
+    end
+
+    @testset "custom name map" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        op = identity_operator(g, s)
+
+        # A user-supplied map dictates the bra-layer name for each renamed index.
+        custom = map(uniquename, keys(tn.dimname_vertices))
+        qf = quadraticformnetwork(tn, op, custom)
+
+        lname = name(l[NamedEdge(1 => 2)])
+        @test braname(qf, lname) == custom[lname]
+        @test braname(qf, lname) in name.(inds(conj_bratensor(qf, 2)))
+
+        sname = name(s[2])
+        @test braname(qf, sname) == custom[sname]
+    end
+
+    @testset "`KetView` / `OperatorView` / `BraView`" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        op = identity_operator(g, s)
+        qf = QuadraticFormNetwork(tn, op)
+
+        kv = KetView(qf)
+        ov = OperatorView(qf)
+        bv = BraView(qf)
+
+        # Views share the graph structure of the ket layer.
+        for view in (kv, ov, bv)
+            @test issetequal(vertices(view), vertices(tn))
+            @test issetequal(edges(view), edges(tn))
+            @test !issettable(view)
+            @test !isinsertable(view)
+            @test is_vertex_assigned(view, 1)
+        end
+
+        for v in vertices(tn)
+            @test kv[v] === kettensor(qf, v)
+            @test inds(ov[v]) == inds(operatortensor(qf, v))
+            @test inds(bv[v]) == inds(bratensor(qf, v))
+        end
+    end
+
+    @testset "contraction / physics" begin
+        @testset "identity operator layer reproduces the norm network" begin
+            g = named_grid((2, 2))
+            tn, l, s = random_state(Float64, g)
+            op = identity_operator(g, s)
+
+            @test contract_network(QuadraticFormNetwork(tn, op))[] ≈
+                contract_network(NormNetwork(tn))[]
+        end
+
+        @testset "$T" for T in (Float64, ComplexF64)
+            g = named_grid((2, 2))
+            tn, l, s = random_state(T, g)
+
+            # A product of on-site matrices, contracted densely for the reference value.
+            mats = Dict(v => randn(T, 2, 2) for v in vertices(g))
+            op = product_operator(v -> mats[v], g, s)
+            qf = QuadraticFormNetwork(tn, op)
+
+            # ⟨ψ|O|ψ⟩ built by applying each on-site matrix to the dense ket and
+            # overlapping with the dense bra.
+            psi = prod(tn)
+            ket = psi
+            for v in vertices(g)
+                out = Index(2)
+                gate = ITensor(mats[v], (out, s[v]))
+                ket = rename(gate * ket, name(out) => name(s[v]))
+            end
+            @test contract_network(qf)[] ≈ (conj(psi) * ket)[]
+        end
+
+        @testset "scaled identity" begin
+            g = named_path_graph(3)
+            tn, l, s = random_state(Float64, g)
+            op = product_operator(v -> 2.0 * Matrix(1.0I, 2, 2), g, s)
+
+            # A factor of 2 at each of the three vertices scales the norm by 2³.
+            @test contract_network(QuadraticFormNetwork(tn, op))[] ≈
+                8 * contract_network(NormNetwork(tn))[]
+        end
+    end
+end

From 73b45b7aff01d5d8fbd44ce944d918856d42a3c2 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Wed, 9 Sep 2026 11:02:54 -0400
Subject: [PATCH 02/18] Move the bilinear-form networks and views to
 `src/bilinearforms/`

The five files are byte-identical at their new paths; only the `include` paths in
`src/ITensorNetworksNext.jl` change. `src/normnetworkview.jl` becomes
`bilinearformnetworkview.jl`, since it now holds the views of both `NormNetwork`
and `QuadraticFormNetwork`.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/ITensorNetworksNext.jl                             | 10 +++++-----
 src/{ => bilinearforms}/abstractbilinearformnetwork.jl |  0
 .../abstractbilinearformnetworkview.jl                 |  0
 .../bilinearformnetworkview.jl}                        |  0
 src/{ => bilinearforms}/normnetwork.jl                 |  0
 src/{ => bilinearforms}/quadraticformnetwork.jl        |  0
 6 files changed, 5 insertions(+), 5 deletions(-)
 rename src/{ => bilinearforms}/abstractbilinearformnetwork.jl (100%)
 rename src/{ => bilinearforms}/abstractbilinearformnetworkview.jl (100%)
 rename src/{normnetworkview.jl => bilinearforms/bilinearformnetworkview.jl} (100%)
 rename src/{ => bilinearforms}/normnetwork.jl (100%)
 rename src/{ => bilinearforms}/quadraticformnetwork.jl (100%)

diff --git a/src/ITensorNetworksNext.jl b/src/ITensorNetworksNext.jl
index 6a24040f..6b851182 100644
--- a/src/ITensorNetworksNext.jl
+++ b/src/ITensorNetworksNext.jl
@@ -13,11 +13,11 @@ include("select_algorithm.jl")
 include("AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl")
 include("abstracttensornetwork.jl")
 include("tensornetwork.jl")
-include("abstractbilinearformnetwork.jl")
-include("abstractbilinearformnetworkview.jl")
-include("normnetwork.jl")
-include("quadraticformnetwork.jl")
-include("normnetworkview.jl")
+include("bilinearforms/abstractbilinearformnetwork.jl")
+include("bilinearforms/abstractbilinearformnetworkview.jl")
+include("bilinearforms/normnetwork.jl")
+include("bilinearforms/quadraticformnetwork.jl")
+include("bilinearforms/bilinearformnetworkview.jl")
 include("ITensorNetworkGenerators/ITensorNetworkGenerators.jl")
 include("contract_network.jl")
 
diff --git a/src/abstractbilinearformnetwork.jl b/src/bilinearforms/abstractbilinearformnetwork.jl
similarity index 100%
rename from src/abstractbilinearformnetwork.jl
rename to src/bilinearforms/abstractbilinearformnetwork.jl
diff --git a/src/abstractbilinearformnetworkview.jl b/src/bilinearforms/abstractbilinearformnetworkview.jl
similarity index 100%
rename from src/abstractbilinearformnetworkview.jl
rename to src/bilinearforms/abstractbilinearformnetworkview.jl
diff --git a/src/normnetworkview.jl b/src/bilinearforms/bilinearformnetworkview.jl
similarity index 100%
rename from src/normnetworkview.jl
rename to src/bilinearforms/bilinearformnetworkview.jl
diff --git a/src/normnetwork.jl b/src/bilinearforms/normnetwork.jl
similarity index 100%
rename from src/normnetwork.jl
rename to src/bilinearforms/normnetwork.jl
diff --git a/src/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl
similarity index 100%
rename from src/quadraticformnetwork.jl
rename to src/bilinearforms/quadraticformnetwork.jl

From 1e87e8b70000fbe36f52fd972a720db1de463d13 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Wed, 9 Sep 2026 13:12:18 -0400
Subject: [PATCH 03/18] Reject `NormNetwork` in `OperatorView` with an
 `ArgumentError`

`OperatorView` now accepts any `AbstractBilinearFormNetwork` rather than only a
`QuadraticFormNetwork`, and a second constructor throws an `ArgumentError` for a
`NormNetwork`, which reports the reason rather than surfacing a `MethodError`.

A `NormNetwork` maps every site index name to itself in its bra layer, so an
identity operator between the layers would need its output name to equal its
input name, and `ITensor` rejects that: "Dimension names must be distinct".

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/bilinearforms/bilinearformnetworkview.jl | 10 ++++++++--
 test/test_quadraticformnetwork.jl            |  3 +++
 2 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/src/bilinearforms/bilinearformnetworkview.jl b/src/bilinearforms/bilinearformnetworkview.jl
index 3b197aeb..4c940778 100644
--- a/src/bilinearforms/bilinearformnetworkview.jl
+++ b/src/bilinearforms/bilinearformnetworkview.jl
@@ -34,14 +34,20 @@ struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
     end
 end
 
-struct OperatorView{T, V, I, P <: QuadraticFormNetwork{T, V, I}} <:
+struct OperatorView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
     AbstractBilinearFormNetworkView{T, V, I}
     parent::P
-    function OperatorView(parent::QuadraticFormNetwork{T, V, I}) where {T, V, I}
+    function OperatorView(parent::AbstractBilinearFormNetwork{T, V, I}) where {T, V, I}
         return new{T, V, I, typeof(parent)}(parent)
     end
 end
 
+# A `NormNetwork`'s bra and ket layers share their site index names, so there is no distinct
+# bra site index for an identity operator to carry as its output name.
+function OperatorView(::NormNetwork)
+    throw(ArgumentError("a `NormNetwork` has no operator layer to view."))
+end
+
 Base.parent(nnv::KetView) = nnv.parent
 Base.parent(nnv::BraView) = nnv.parent
 Base.parent(nnv::OperatorView) = nnv.parent
diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl
index f303b259..cd90bad1 100644
--- a/test/test_quadraticformnetwork.jl
+++ b/test/test_quadraticformnetwork.jl
@@ -148,6 +148,9 @@ identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s;
             @test is_vertex_assigned(view, 1)
         end
 
+        # A `NormNetwork` has no operator layer, so `OperatorView` rejects one.
+        @test_throws ArgumentError OperatorView(NormNetwork(tn))
+
         for v in vertices(tn)
             @test kv[v] === kettensor(qf, v)
             @test inds(ov[v]) == inds(operatortensor(qf, v))

From d0d48f63c1416fcf9e2c4d26873e030e3ac4f2f6 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Wed, 9 Sep 2026 13:57:58 -0400
Subject: [PATCH 04/18] Define the NamedGraphs vertex encoding once on the
 abstract view type

`encoded_vertex`, `decoded_vertex` and `encoded_graph` move from a
`for View in (:KetView, :BraView)` loop in `bilinearformnetworkview.jl` to
`AbstractBilinearFormNetworkView` in `abstractbilinearformnetworkview.jl`,
replacing the `vertex_positions`/`ordered_vertices`/`position_graph` methods that
NamedGraphs 0.13 removed. `OperatorView` is now covered as well, as it was by the
methods being replaced. The loop also `@eval`ed against `BraView` above that
type's own definition, so including the file threw `UndefVarError: BraView`.

`QuadraticFormNetwork` gets the same rename; it still defined the three removed
names, so the package did not load. `test_quadraticformnetwork.jl` takes
`incident_edges`, `named_grid` and `named_path_graph` from `NamedGraphs` itself,
since the `GraphsExtensions` and `NamedGraphGenerators` submodules are gone.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 .../abstractbilinearformnetworkview.jl           | 14 ++++++--------
 src/bilinearforms/bilinearformnetworkview.jl     | 16 ----------------
 src/bilinearforms/quadraticformnetwork.jl        | 10 +++++++---
 test/test_quadraticformnetwork.jl                |  4 +---
 4 files changed, 14 insertions(+), 30 deletions(-)

diff --git a/src/bilinearforms/abstractbilinearformnetworkview.jl b/src/bilinearforms/abstractbilinearformnetworkview.jl
index 828cbce1..262bc4e7 100644
--- a/src/bilinearforms/abstractbilinearformnetworkview.jl
+++ b/src/bilinearforms/abstractbilinearformnetworkview.jl
@@ -1,7 +1,7 @@
 using DataGraphs: DataGraphs, is_vertex_assigned
 using Dictionaries: Dictionaries, isinsertable, issettable
 using Graphs: Graphs, edges, vertices
-using NamedGraphs: NamedGraphs, ordered_vertices, position_graph, vertex_positions
+using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
 
 """
     abstract type AbstractBilinearFormNetworkView{T, V, I} <: AbstractITensorNetwork{T, V}
@@ -21,15 +21,13 @@ Graphs.vertices(nnv::AbstractBilinearFormNetworkView) = vertices(parent(nnv))
 
 # ==================================== NamedGraphs.jl ==================================== #
 
-function NamedGraphs.vertex_positions(nnv::AbstractBilinearFormNetworkView)
-    return vertex_positions(parent(nnv))
+function NamedGraphs.encoded_vertex(nnv::AbstractBilinearFormNetworkView, vertex)
+    return encoded_vertex(parent(nnv), vertex)
 end
-function NamedGraphs.ordered_vertices(nnv::AbstractBilinearFormNetworkView)
-    return ordered_vertices(parent(nnv))
-end
-function NamedGraphs.position_graph(nnv::AbstractBilinearFormNetworkView)
-    return position_graph(parent(nnv))
+function NamedGraphs.decoded_vertex(nnv::AbstractBilinearFormNetworkView, code::Integer)
+    return decoded_vertex(parent(nnv), code)
 end
+NamedGraphs.encoded_graph(nnv::AbstractBilinearFormNetworkView) = encoded_graph(parent(nnv))
 
 # ==================================== DataGraphs.jl ===================================== #
 
diff --git a/src/bilinearforms/bilinearformnetworkview.jl b/src/bilinearforms/bilinearformnetworkview.jl
index 4c940778..dbeb81bf 100644
--- a/src/bilinearforms/bilinearformnetworkview.jl
+++ b/src/bilinearforms/bilinearformnetworkview.jl
@@ -1,7 +1,6 @@
 using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
 using Dictionaries: Dictionaries, isinsertable, issettable
 using Graphs: Graphs, edges, vertices
-using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
 
 struct KetView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
     AbstractBilinearFormNetworkView{T, V, I}
@@ -11,21 +10,6 @@ struct KetView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
     end
 end
 
-# ==================================== NamedGraphs.jl ==================================== #
-
-for View in (:KetView, :BraView)
-    @eval begin
-        function NamedGraphs.encoded_vertex(nnv::$View, vertex)
-            return encoded_vertex(nnv.parent, vertex)
-        end
-        function NamedGraphs.decoded_vertex(nnv::$View, code::Integer)
-            return decoded_vertex(nnv.parent, code)
-        end
-
-        NamedGraphs.encoded_graph(nnv::$View) = encoded_graph(nnv.parent)
-    end
-end
-
 struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
     AbstractBilinearFormNetworkView{T, V, I}
     parent::P
diff --git a/src/bilinearforms/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl
index 1b803476..39c2dd0e 100644
--- a/src/bilinearforms/quadraticformnetwork.jl
+++ b/src/bilinearforms/quadraticformnetwork.jl
@@ -54,9 +54,13 @@ Graphs.vertices(qf::QuadraticFormNetwork) = vertices(qf.ket)
 
 # ==================================== NamedGraphs.jl ==================================== #
 
-NamedGraphs.vertex_positions(qf::QuadraticFormNetwork) = vertex_positions(qf.ket)
-NamedGraphs.ordered_vertices(qf::QuadraticFormNetwork) = ordered_vertices(qf.ket)
-NamedGraphs.position_graph(qf::QuadraticFormNetwork) = position_graph(qf.ket)
+function NamedGraphs.encoded_vertex(qf::QuadraticFormNetwork, vertex)
+    return encoded_vertex(qf.ket, vertex)
+end
+function NamedGraphs.decoded_vertex(qf::QuadraticFormNetwork, code::Integer)
+    return decoded_vertex(qf.ket, code)
+end
+NamedGraphs.encoded_graph(qf::QuadraticFormNetwork) = encoded_graph(qf.ket)
 
 # ==================================== DataGraphs.jl ===================================== #
 
diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl
index cd90bad1..7c036bd2 100644
--- a/test/test_quadraticformnetwork.jl
+++ b/test/test_quadraticformnetwork.jl
@@ -7,9 +7,7 @@ using ITensorNetworksNext: BraView, ITensorNetwork, KetView, NormNetwork, Operat
     QuadraticFormNetwork, braname, bratensor, conj_bratensor, contract_network, indmap,
     kettensor, operatortensor, quadraticformnetwork, tensornetwork
 using LinearAlgebra: I, norm
-using NamedGraphs.GraphsExtensions: incident_edges
-using NamedGraphs.NamedGraphGenerators: named_grid, named_path_graph
-using NamedGraphs: NamedEdge
+using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph
 using Test: @test, @test_throws, @testset
 
 # Build a random `ITensorNetwork` state on the graph `g` with site dimension `d` and

From e383b2a584aa810a76b762c8173bcfb81a42bbed Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 14 Sep 2026 11:56:24 -0400
Subject: [PATCH 05/18] Add `ITensorNetworkOperator` and use it as the operator
 layer

`ITensorNetworkOperator` is the network equivalent of `ITensorBase.ITensorOperator`:
a tensor network together with a positional pairing of its dangling index names,
`outputnames[i]` with `inputnames[i]`. The constructor, reached through
`ITensorBase.operator(tn, out, in)`, throws an `ArgumentError` when the two name
lists differ in length and when a named index is not dangling, since
`dimnamevertices` must return exactly one vertex for an operator leg.

A pair may straddle two vertices, as it does for a swap or a translation, so
`getindex` wraps each vertex tensor as an `ITensorOperator` carrying only the pairs
whose two halves both sit on that vertex; the remaining legs are dangling on the
wrapper, which `NamedTensorOperator` permits for names in neither list. The wrapper
is built even when no pair is local, so `eltype` is `NamedTensorOperator{I, T}` for
every vertex.

`QuadraticFormNetwork` takes one as its operator layer. Its `acted` set is now
`Set{I}(inputnames(operator))` rather than a loop over per-vertex input names, and
`operatortensor` reads the network-level pairing: reading `qf.operator[vertex]`
would skip a crossing pair's output leg, which is dangling on that wrapper, and the
bra layer would not meet it. `operatortensor` therefore returns a plain tensor and
`get_vertex_data` drops its `state` call.

`supportof` gains a method for the new type. Its body moves to
`supportof_dimnames(tn, opnames)` in `abstracttensornetwork.jl` so both methods
share it; a `Union` signature is not possible there because
`ITensorNetworkOperator` is not yet defined at that point in the include order.

`replacedimnames` is not implemented for the new type: `ITensorNetwork` has no
`replacedimnames` method to delegate to.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 src/ITensorNetworksNext.jl                |   1 +
 src/abstracttensornetwork.jl              |   8 +-
 src/bilinearforms/quadraticformnetwork.jl |  55 ++++++-----
 src/itensornetworkoperator.jl             | 114 +++++++++++++++++++++
 test/test_itensornetworkoperator.jl       | 115 ++++++++++++++++++++++
 test/test_quadraticformnetwork.jl         |  16 +--
 6 files changed, 276 insertions(+), 33 deletions(-)
 create mode 100644 src/itensornetworkoperator.jl
 create mode 100644 test/test_itensornetworkoperator.jl

diff --git a/src/ITensorNetworksNext.jl b/src/ITensorNetworksNext.jl
index 6b851182..87c20a5a 100644
--- a/src/ITensorNetworksNext.jl
+++ b/src/ITensorNetworksNext.jl
@@ -13,6 +13,7 @@ include("select_algorithm.jl")
 include("AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl")
 include("abstracttensornetwork.jl")
 include("tensornetwork.jl")
+include("itensornetworkoperator.jl")
 include("bilinearforms/abstractbilinearformnetwork.jl")
 include("bilinearforms/abstractbilinearformnetworkview.jl")
 include("bilinearforms/normnetwork.jl")
diff --git a/src/abstracttensornetwork.jl b/src/abstracttensornetwork.jl
index 341d7a6d..a8d959fa 100644
--- a/src/abstracttensornetwork.jl
+++ b/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
diff --git a/src/bilinearforms/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl
index 39c2dd0e..74808803 100644
--- a/src/bilinearforms/quadraticformnetwork.jl
+++ b/src/bilinearforms/quadraticformnetwork.jl
@@ -1,49 +1,47 @@
 using Dictionaries: Dictionary
 using ITensorBase:
-    LazyNamedTensor, inputnames, lazy, outputnames, rename, state, uniquename
+    LazyNamedTensor, inputnames, lazy, names, outputnames, rename, state, uniquename
 using ITensorNetworksNext
 
 """
     struct QuadraticFormNetwork{T, V, I, O} <: AbstractBilinearFormNetwork{T, V, I}
 
 Lazy wrapper representing the quadratic form `⟨tn|op|tn⟩` of `tn::ITensorNetwork{T, V, I}`
-sandwiched around the operator layer `op::ITensorNetwork{O, V, I}`, together with a per-index
+sandwiched around the operator layer `op::ITensorNetworkOperator`, together with a per-index
 ket→bra name mapping that, for each index in the ket layer, defines the name of the
 corresponding index in the bra layer.
 
-The operator layer must be defined on every vertex of `tn`. At each vertex the operator's
-input names are ket index names, and each paired output name is renamed to the bra layer, so
-the operator's remaining index names must be distinct from every index name of `tn`.
+The operator layer must be defined on every vertex of `tn`. Its input names are ket index
+names, and each paired output name is renamed to the bra layer, so the operator's remaining
+index names must be distinct from every index name of `tn`.
 """
-struct QuadraticFormNetwork{T, V, I, O} <: AbstractBilinearFormNetwork{T, V, I}
+struct QuadraticFormNetwork{T, V, I, O <: ITensorNetworkOperator} <:
+    AbstractBilinearFormNetwork{T, V, I}
     ket::ITensorNetwork{T, V, I}
-    operator::ITensorNetwork{O, V, I}
+    operator::O
     braname::Dictionary{I, I}
     function QuadraticFormNetwork(
             ket::ITensorNetwork{T, V, I},
-            operator::ITensorNetwork{O, V, I},
+            operator::ITensorNetworkOperator,
             map::Dictionary{I, I}
-        ) where {T, V, I, O}
+        ) where {T, V, I}
         if !issetequal(vertices(operator), vertices(ket))
             error("the operator layer must be defined on every vertex of the ket layer.")
         end
-        acted = Set{I}()
-        for vertex in vertices(operator)
-            union!(acted, inputnames(operator[vertex]))
-        end
+        acted = Set{I}(inputnames(operator))
         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 new{T, V, I, O}(ket, operator, braname)
+        return new{T, V, I, typeof(operator)}(ket, operator, braname)
     end
 end
 
 Base.eltype(::Type{<:QuadraticFormNetwork{T, V, I}}) where {T, V, I} = LazyNamedTensor{I, T}
 
-function QuadraticFormNetwork(ket::ITensorNetwork, operator::ITensorNetwork)
+function QuadraticFormNetwork(ket::ITensorNetwork, operator::ITensorNetworkOperator)
     return QuadraticFormNetwork(ket, operator, map(uniquename, keys(ket.dimname_vertices)))
 end
 
@@ -69,7 +67,7 @@ function DataGraphs.get_vertex_data(qf::QuadraticFormNetwork, vertex)
     O = operatortensor(qf, vertex)
     B = conj_bratensor(qf, vertex)
     # TODO: implement and use a lazy `conj` via `LazyNamedDimsArrays` here?
-    return lazy(A) * lazy(state(O)) * lazy(conj(B))
+    return lazy(A) * lazy(O) * lazy(conj(B))
 end
 
 function DataGraphs.is_vertex_assigned(qf::QuadraticFormNetwork, vertex)
@@ -90,17 +88,22 @@ end
 kettensor(qf::QuadraticFormNetwork, vertex) = qf.ket[vertex]
 
 # Each output name is renamed to the bra name of the input name it is paired with, so the
-# operator's output legs meet the bra layer and its input legs meet the ket layer.
+# operator's output legs meet the bra layer and its input legs meet the ket layer. The pairing
+# is read from the operator network rather than from `qf.operator[vertex]`, whose wrapper drops
+# a pair whose input sits on another vertex.
 function operatortensor(qf::QuadraticFormNetwork, vertex)
-    tensor = qf.operator[vertex]
-    names = map(inputnames(tensor)) do n
-        return braname(qf, n)
-    end
-    return rename(tensor, (outputnames(tensor) .=> names)...)
+    tensor = state(qf.operator)[vertex]
+    tensor_names = names(tensor)
+    replacements = [
+        output => braname(qf, input) for
+            (output, input) in zip(outputnames(qf.operator), inputnames(qf.operator))
+            if output in tensor_names
+    ]
+    return rename(tensor, replacements...)
 end
 
 """
-    quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetwork, [braname])
+    quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetworkOperator, [braname])
 
 Build the triple-layer network `⟨tn|op|tn⟩`, represented lazily as a
 `QuadraticFormNetwork` object. The optional third argument `braname` should implement
@@ -108,7 +111,9 @@ Build the triple-layer network `⟨tn|op|tn⟩`, represented lazily as a
 every dimension name of `tn` the operator acts on. If this is not specified, then a name is
 generated via the `ITensorBase.uniquename` function.
 """
-quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetwork) = QuadraticFormNetwork(tn, op)
-function quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetwork, braname)
+function quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetworkOperator)
+    return QuadraticFormNetwork(tn, op)
+end
+function quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetworkOperator, braname)
     return QuadraticFormNetwork(tn, op, braname)
 end
diff --git a/src/itensornetworkoperator.jl b/src/itensornetworkoperator.jl
new file mode 100644
index 00000000..c64d1208
--- /dev/null
+++ b/src/itensornetworkoperator.jl
@@ -0,0 +1,114 @@
+using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
+using Dictionaries: Dictionaries
+using Graphs: Graphs, AbstractGraph, edges, vertices
+using ITensorBase: ITensorBase, NamedTensorOperator, inputnames, names, nametype, operator,
+    outputnames, state
+using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
+
+"""
+    struct ITensorNetworkOperator{T, V, I, P} <: AbstractITensorNetwork{T, V}
+
+The network equivalent of `ITensorBase.ITensorOperator`: a tensor network of type
+`P <: AbstractITensorNetwork{T, V}` together with a pairing of its dangling index names,
+where `outputnames[i]` is paired with `inputnames[i]`. Applying the operator contracts over
+the input names and leaves the output names.
+
+A pair may straddle two vertices, as it does for a swap or a translation. Indexing returns
+the vertex tensor wrapped as an `ITensorOperator` carrying only the pairs whose two halves
+both sit on that vertex; the remaining legs are dangling on the wrapper.
+"""
+struct ITensorNetworkOperator{T, V, I, P <: AbstractITensorNetwork{T, V}} <:
+    AbstractITensorNetwork{T, V}
+    parent::P
+    outputnames::Vector{I}
+    inputnames::Vector{I}
+    function ITensorNetworkOperator(
+            parent::AbstractITensorNetwork{T, V}, outputnames, inputnames
+        ) where {T, V}
+        I = nametype(T)
+        outputnames = collect(I, outputnames)
+        inputnames = collect(I, inputnames)
+        if length(outputnames) != length(inputnames)
+            throw(
+                ArgumentError(
+                    "Operator `outputnames` and `inputnames` must have equal length " *
+                        "(positional pairing), got $(length(outputnames)) and " *
+                        "$(length(inputnames))."
+                )
+            )
+        end
+        for opname in Iterators.flatten((outputnames, inputnames))
+            nvertices = length(dimnamevertices(parent, opname))
+            if nvertices != 1
+                throw(
+                    ArgumentError(
+                        "operator dim name $opname is associated with $nvertices vertices " *
+                            "in the tensor network; an operator leg must be a dangling index."
+                    )
+                )
+            end
+        end
+        return new{T, V, I, typeof(parent)}(parent, outputnames, inputnames)
+    end
+end
+
+function ITensorBase.operator(tn::AbstractITensorNetwork, outputnames, inputnames)
+    return ITensorNetworkOperator(tn, outputnames, inputnames)
+end
+
+function Base.eltype(::Type{<:ITensorNetworkOperator{T, V, I}}) where {T, V, I}
+    return NamedTensorOperator{I, T}
+end
+
+# ====================================== Graphs.jl ======================================= #
+
+Graphs.edges(op::ITensorNetworkOperator) = edges(state(op))
+Graphs.vertices(op::ITensorNetworkOperator) = vertices(state(op))
+
+# ==================================== NamedGraphs.jl ==================================== #
+
+function NamedGraphs.encoded_vertex(op::ITensorNetworkOperator, vertex)
+    return encoded_vertex(state(op), vertex)
+end
+function NamedGraphs.decoded_vertex(op::ITensorNetworkOperator, code::Integer)
+    return decoded_vertex(state(op), code)
+end
+NamedGraphs.encoded_graph(op::ITensorNetworkOperator) = encoded_graph(state(op))
+
+# ==================================== DataGraphs.jl ===================================== #
+
+# A pair whose two halves sit on different vertices has no bijection to give this vertex, so
+# only the pairs local to `vertex` become the wrapper's pairing and the rest stay dangling.
+function DataGraphs.get_vertex_data(op::ITensorNetworkOperator, vertex)
+    tensor = state(op)[vertex]
+    tensor_names = names(tensor)
+    outputs = similar(outputnames(op), 0)
+    inputs = similar(inputnames(op), 0)
+    for (output, input) in zip(outputnames(op), inputnames(op))
+        if output in tensor_names && input in tensor_names
+            push!(outputs, output)
+            push!(inputs, input)
+        end
+    end
+    return operator(tensor, outputs, inputs)
+end
+
+function DataGraphs.is_vertex_assigned(op::ITensorNetworkOperator, vertex)
+    return is_vertex_assigned(state(op), vertex)
+end
+
+# =================================== Dictionaries.jl ==================================== #
+
+Dictionaries.issettable(::ITensorNetworkOperator) = false
+Dictionaries.isinsertable(::ITensorNetworkOperator) = false
+
+# ====================================== interface ======================================= #
+
+ITensorBase.state(op::ITensorNetworkOperator) = op.parent
+Base.parent(op::ITensorNetworkOperator) = state(op)
+ITensorBase.outputnames(op::ITensorNetworkOperator) = op.outputnames
+ITensorBase.inputnames(op::ITensorNetworkOperator) = op.inputnames
+
+function operator_support(tn::AbstractGraph, op::ITensorNetworkOperator)
+    return operator_support_names(tn, inputnames(op))
+end
diff --git a/test/test_itensornetworkoperator.jl b/test/test_itensornetworkoperator.jl
new file mode 100644
index 00000000..68b0b6c4
--- /dev/null
+++ b/test/test_itensornetworkoperator.jl
@@ -0,0 +1,115 @@
+using DataGraphs: is_vertex_assigned
+using Dictionaries: isinsertable, issettable
+using Graphs: edges, vertices
+using ITensorBase: ITensor, Index, inputnames, name, names, operator, outputnames, state
+using ITensorNetworksNext:
+    ITensorNetwork, ITensorNetworkOperator, operator_support, tensornetwork
+using NamedGraphs: incident_edges, named_path_graph
+using Test: @test, @test_throws, @testset
+
+# A bondless operator network on `g`, with one input/output pair per vertex.
+function local_operator(g; d = 2)
+    inp = Dict(v => Index(d) for v in vertices(g))
+    out = Dict(v => Index(d) for v in vertices(g))
+    tn = tensornetwork(vertices(g)) do v
+        return randn((out[v], inp[v]))
+    end
+    vs = collect(vertices(g))
+    op = operator(tn, [name(out[v]) for v in vs], [name(inp[v]) for v in vs])
+    return op, inp, out
+end
+
+@testset "`ITensorNetworkOperator`" begin
+    @testset "Basics" begin
+        g = named_path_graph(3)
+        op, inp, out = local_operator(g)
+
+        @test op isa ITensorNetworkOperator
+
+        # The operator shares the graph structure of the network it wraps.
+        @test issetequal(vertices(op), vertices(state(op)))
+        @test issetequal(edges(op), edges(state(op)))
+        @test parent(op) === state(op)
+
+        # The pairing is positional: `outputnames[i]` goes with `inputnames[i]`.
+        vs = collect(vertices(g))
+        @test outputnames(op) == [name(out[v]) for v in vs]
+        @test inputnames(op) == [name(inp[v]) for v in vs]
+
+        @test is_vertex_assigned(op, first(vs))
+
+        # It is a lazy wrapper, so it is neither settable nor insertable.
+        @test !issettable(op)
+        @test !isinsertable(op)
+    end
+
+    @testset "`getindex` wraps the vertex tensor" begin
+        g = named_path_graph(3)
+        op, inp, out = local_operator(g)
+
+        # Both halves of each pair sit on the same vertex here, so the wrapper carries them.
+        for v in vertices(g)
+            @test outputnames(op[v]) == [name(out[v])]
+            @test inputnames(op[v]) == [name(inp[v])]
+            @test state(op[v]) === state(op)[v]
+        end
+
+        # `eltype` is the type of the wrapped vertex data.
+        @test eltype(op) === typeof(op[first(vertices(g))])
+    end
+
+    @testset "a pair may straddle two vertices" begin
+        # A swap: the output at vertex 1 is paired with the input at vertex 2, and vice versa.
+        a, b = Index(2), Index(2)
+        a′, b′ = Index(2), Index(2)
+        tn = ITensorNetwork(Dict(1 => randn((a′, a)), 2 => randn((b′, b))))
+        op = operator(tn, [name(a′), name(b′)], [name(b), name(a)])
+
+        @test outputnames(op) == [name(a′), name(b′)]
+        @test inputnames(op) == [name(b), name(a)]
+
+        # Neither pair is local to a vertex, so each vertex wrapper has an empty pairing and
+        # both of its legs stay dangling.
+        for v in (1, 2)
+            @test isempty(outputnames(op[v]))
+            @test isempty(inputnames(op[v]))
+            @test issetequal(names(op[v]), names(state(op)[v]))
+        end
+    end
+
+    @testset "constructor validation" begin
+        g = named_path_graph(2)
+        op, inp, out = local_operator(g)
+        tn = state(op)
+        vs = collect(vertices(g))
+
+        # The pairing is a bijection, so the two name lists must have equal length.
+        @test_throws ArgumentError operator(tn, [name(out[vs[1]])], [])
+        @test_throws ArgumentError operator(
+            tn, [name(out[v]) for v in vs], [name(inp[vs[1]])]
+        )
+
+        # An operator leg must be dangling: a link name touches two vertices.
+        l = Index(2)
+        linked = tensornetwork(vertices(g)) do v
+            return randn((inp[v], l))
+        end
+        @test_throws ArgumentError operator(linked, [name(l)], [name(inp[vs[1]])])
+    end
+
+    @testset "`operator_support`" begin
+        g = named_path_graph(3)
+        s = Dict(v => Index(2) for v in vertices(g))
+        state_tn = tensornetwork(vertices(g)) do v
+            is = map(e -> Index(2), incident_edges(g, v))
+            return randn((s[v], is...))
+        end
+
+        # An operator acting on sites 1 and 3 is supported on exactly those vertices.
+        out1, out3 = Index(2), Index(2)
+        optn = ITensorNetwork(Dict(1 => randn((out1, s[1])), 3 => randn((out3, s[3]))))
+        op = operator(optn, [name(out1), name(out3)], [name(s[1]), name(s[3])])
+
+        @test issetequal(operator_support(state_tn, op), Set([1, 3]))
+    end
+end
diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl
index 7c036bd2..1866e985 100644
--- a/test/test_quadraticformnetwork.jl
+++ b/test/test_quadraticformnetwork.jl
@@ -1,7 +1,7 @@
 using DataGraphs: is_vertex_assigned
 using Dictionaries: isinsertable, issettable
 using Graphs: edges, vertices
-using ITensorBase: ITensor, Index, IndexName, conj, inds, inputnames, name, operator,
+using ITensorBase: ITensor, Index, IndexName, conj, inds, inputnames, name, names, operator,
     outputnames, rename, setname, uniquename
 using ITensorNetworksNext: BraView, ITensorNetwork, KetView, NormNetwork, OperatorView,
     QuadraticFormNetwork, braname, bratensor, conj_bratensor, contract_network, indmap,
@@ -26,11 +26,12 @@ end
 # Build a bondless operator layer on the vertices of `g`, with `f(v)` supplying the matrix
 # acting on the site index `s[v]`.
 function product_operator(f, g, s; d = 2)
-    op = tensornetwork(vertices(g)) do v
-        out = Index(d)
-        return operator(ITensor(f(v), (out, s[v])), (name(out),), (name(s[v]),))
+    out = Dict(v => Index(d) for v in vertices(g))
+    tn = tensornetwork(vertices(g)) do v
+        return ITensor(f(v), (out[v], s[v]))
     end
-    return op
+    vs = collect(vertices(g))
+    return operator(tn, [name(out[v]) for v in vs], [name(s[v]) for v in vs])
 end
 
 identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s; d)
@@ -94,8 +95,9 @@ identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s;
         # The operator's input name meets the ket and its output name is renamed to meet
         # the bra.
         o = operatortensor(qf, 2)
-        @test inputnames(o) == [sname]
-        @test outputnames(o) == [braname(qf, sname)]
+        @test sname in names(o)
+        @test braname(qf, sname) in names(o)
+        @test inputnames(op) == [name(s[v]) for v in vertices(g)]
 
         # `bratensor` is the elementwise conjugate of `conj_bratensor` and carries the same
         # indices.

From 023d7b89f1e323b77bc8808143c606025e6e91c3 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 11:54:10 -0400
Subject: [PATCH 06/18] Replace `KetView` and `OperatorView` with `ketnetwork`
 and `operatornetwork`

`ketnetwork` and `operatornetwork` return the stored ket network and `ITensorNetworkOperator`.
`branetwork` defaults to a `BraView`, since a bilinear-form network does not store its bra layer.
The belief propagation code for `NormNetwork` now reads its link names through these accessors.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/beliefpropagation/beliefpropagation.jl    |  6 +--
 src/beliefpropagation/messagecache.jl         | 12 +++---
 .../abstractbilinearformnetwork.jl            | 25 +++++++++++-
 src/bilinearforms/bilinearformnetworkview.jl  | 26 -------------
 src/bilinearforms/normnetwork.jl              |  1 +
 src/bilinearforms/quadraticformnetwork.jl     |  2 +
 test/test_normnetwork.jl                      | 25 +++++-------
 test/test_quadraticformnetwork.jl             | 38 ++++++++-----------
 8 files changed, 61 insertions(+), 74 deletions(-)

diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl
index 3741fb3a..7dba8aa2 100644
--- a/src/beliefpropagation/beliefpropagation.jl
+++ b/src/beliefpropagation/beliefpropagation.jl
@@ -275,9 +275,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)
diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index 650a5858..96e7f417 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/src/beliefpropagation/messagecache.jl
@@ -193,16 +193,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
diff --git a/src/bilinearforms/abstractbilinearformnetwork.jl b/src/bilinearforms/abstractbilinearformnetwork.jl
index 6d1e894e..6d0dd6c2 100644
--- a/src/bilinearforms/abstractbilinearformnetwork.jl
+++ b/src/bilinearforms/abstractbilinearformnetwork.jl
@@ -10,7 +10,8 @@ Supertype of the lazy multi-layer networks built from a ket layer of type
 A subtype supplies its own graph structure and implements [`braname`](@ref) together with one
 accessor per layer: [`kettensor`](@ref), [`bratensor`](@ref) and, where the subtype has an
 operator layer, [`operatortensor`](@ref). `bratensor` has a default built from `kettensor`
-and `braname`.
+and `braname`. 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
 
@@ -54,4 +55,26 @@ The bra-layer tensor at `vertex`.
 """
 bratensor(bn::AbstractBilinearFormNetwork, vertex) = conj(conj_bratensor(bn, vertex))
 
+"""
+    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)
+
 indmap(bn::AbstractBilinearFormNetwork, ind) = setname(conj(ind), braname(bn, name(ind)))
diff --git a/src/bilinearforms/bilinearformnetworkview.jl b/src/bilinearforms/bilinearformnetworkview.jl
index dbeb81bf..d9622fee 100644
--- a/src/bilinearforms/bilinearformnetworkview.jl
+++ b/src/bilinearforms/bilinearformnetworkview.jl
@@ -2,14 +2,6 @@ using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
 using Dictionaries: Dictionaries, isinsertable, issettable
 using Graphs: Graphs, edges, vertices
 
-struct KetView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
-    AbstractBilinearFormNetworkView{T, V, I}
-    parent::P
-    function KetView(parent::AbstractBilinearFormNetwork{T, V, I}) where {T, V, I}
-        return new{T, V, I, typeof(parent)}(parent)
-    end
-end
-
 struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
     AbstractBilinearFormNetworkView{T, V, I}
     parent::P
@@ -18,26 +10,8 @@ struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
     end
 end
 
-struct OperatorView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <:
-    AbstractBilinearFormNetworkView{T, V, I}
-    parent::P
-    function OperatorView(parent::AbstractBilinearFormNetwork{T, V, I}) where {T, V, I}
-        return new{T, V, I, typeof(parent)}(parent)
-    end
-end
-
-# A `NormNetwork`'s bra and ket layers share their site index names, so there is no distinct
-# bra site index for an identity operator to carry as its output name.
-function OperatorView(::NormNetwork)
-    throw(ArgumentError("a `NormNetwork` has no operator layer to view."))
-end
-
-Base.parent(nnv::KetView) = nnv.parent
 Base.parent(nnv::BraView) = nnv.parent
-Base.parent(nnv::OperatorView) = nnv.parent
 
 # ==================================== DataGraphs.jl ===================================== #
 
-DataGraphs.get_vertex_data(nnv::KetView, vertex) = kettensor(parent(nnv), vertex)
 DataGraphs.get_vertex_data(nnv::BraView, vertex) = bratensor(parent(nnv), vertex)
-DataGraphs.get_vertex_data(nnv::OperatorView, vertex) = operatortensor(parent(nnv), vertex)
diff --git a/src/bilinearforms/normnetwork.jl b/src/bilinearforms/normnetwork.jl
index 6428644f..03add5a0 100644
--- a/src/bilinearforms/normnetwork.jl
+++ b/src/bilinearforms/normnetwork.jl
@@ -68,6 +68,7 @@ function braname(nn::NormNetwork, name)
 end
 
 kettensor(nn::NormNetwork, vertex) = nn.ket[vertex]
+ketnetwork(nn::NormNetwork) = nn.ket
 
 """
     normnetwork(tn::ITensorNetwork, [braname]) -> NormNetwork
diff --git a/src/bilinearforms/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl
index 74808803..a01b5887 100644
--- a/src/bilinearforms/quadraticformnetwork.jl
+++ b/src/bilinearforms/quadraticformnetwork.jl
@@ -86,6 +86,8 @@ function braname(qf::QuadraticFormNetwork, name)
 end
 
 kettensor(qf::QuadraticFormNetwork, vertex) = qf.ket[vertex]
+ketnetwork(qf::QuadraticFormNetwork) = qf.ket
+operatornetwork(qf::QuadraticFormNetwork) = qf.operator
 
 # Each output name is renamed to the bra name of the input name it is paired with, so the
 # operator's output legs meet the bra layer and its input legs meet the ket layer. The pairing
diff --git a/test/test_normnetwork.jl b/test/test_normnetwork.jl
index 5007be69..6f2b9ff3 100644
--- a/test/test_normnetwork.jl
+++ b/test/test_normnetwork.jl
@@ -3,8 +3,8 @@ using Dictionaries: isinsertable, issettable
 using Graphs: edges, vertices
 using ITensorBase:
     ITensor, Index, IndexName, LazyITensor, conj, inds, name, setname, uniquename
-using ITensorNetworksNext: BraView, ITensorNetwork, KetView, NormNetwork, braname,
-    bratensor, conj_bratensor, contract_network, indmap, kettensor, normnetwork,
+using ITensorNetworksNext: BraView, ITensorNetwork, NormNetwork, braname, branetwork,
+    bratensor, conj_bratensor, contract_network, indmap, ketnetwork, kettensor, normnetwork,
     tensornetwork
 using LinearAlgebra: norm
 using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph
@@ -96,32 +96,25 @@ end
         @test braname(nn, lname) in name.(inds(conj_bratensor(nn, 2)))
     end
 
-    @testset "`KetView` / `BraView`" begin
+    @testset "`ketnetwork` / `branetwork`" begin
         g = named_path_graph(3)
         tn, l, s = random_state(Float64, g)
         nn = NormNetwork(tn)
 
-        kv = KetView(nn)
-        bv = BraView(nn)
+        # The ket layer is the network the norm network was built from.
+        @test ketnetwork(nn) === tn
 
-        # Views share the graph structure of the underlying network.
-        @test issetequal(vertices(kv), vertices(tn))
+        # The bra layer is not stored, so it is a view sharing the ket layer's graph structure.
+        bv = branetwork(nn)
+        @test bv isa BraView
         @test issetequal(vertices(bv), vertices(tn))
-        @test issetequal(edges(kv), edges(tn))
         @test issetequal(edges(bv), edges(tn))
-
-        # The ket view exposes the bare ket tensors; the bra view exposes the bra tensors.
         for v in vertices(tn)
-            @test kv[v] === kettensor(nn, v)
             @test inds(bv[v]) == inds(bratensor(nn, v))
         end
-
-        @test is_vertex_assigned(kv, 1)
         @test is_vertex_assigned(bv, 1)
 
-        # Views inherit the (non-)mutability of their parent norm network.
-        @test !issettable(kv)
-        @test !isinsertable(kv)
+        # The view inherits the (non-)mutability of its parent norm network.
         @test !issettable(bv)
         @test !isinsertable(bv)
     end
diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl
index 1866e985..fd18c2e6 100644
--- a/test/test_quadraticformnetwork.jl
+++ b/test/test_quadraticformnetwork.jl
@@ -3,9 +3,9 @@ using Dictionaries: isinsertable, issettable
 using Graphs: edges, vertices
 using ITensorBase: ITensor, Index, IndexName, conj, inds, inputnames, name, names, operator,
     outputnames, rename, setname, uniquename
-using ITensorNetworksNext: BraView, ITensorNetwork, KetView, NormNetwork, OperatorView,
-    QuadraticFormNetwork, braname, bratensor, conj_bratensor, contract_network, indmap,
-    kettensor, operatortensor, quadraticformnetwork, tensornetwork
+using ITensorNetworksNext: BraView, ITensorNetwork, NormNetwork, QuadraticFormNetwork,
+    braname, branetwork, bratensor, conj_bratensor, contract_network, indmap, ketnetwork,
+    kettensor, operatornetwork, operatortensor, quadraticformnetwork, tensornetwork
 using LinearAlgebra: I, norm
 using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph
 using Test: @test, @test_throws, @testset
@@ -129,31 +129,25 @@ identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s;
         @test braname(qf, sname) == custom[sname]
     end
 
-    @testset "`KetView` / `OperatorView` / `BraView`" begin
+    @testset "`ketnetwork` / `operatornetwork` / `branetwork`" begin
         g = named_path_graph(3)
         tn, l, s = random_state(Float64, g)
         op = identity_operator(g, s)
         qf = QuadraticFormNetwork(tn, op)
 
-        kv = KetView(qf)
-        ov = OperatorView(qf)
-        bv = BraView(qf)
-
-        # Views share the graph structure of the ket layer.
-        for view in (kv, ov, bv)
-            @test issetequal(vertices(view), vertices(tn))
-            @test issetequal(edges(view), edges(tn))
-            @test !issettable(view)
-            @test !isinsertable(view)
-            @test is_vertex_assigned(view, 1)
-        end
-
-        # A `NormNetwork` has no operator layer, so `OperatorView` rejects one.
-        @test_throws ArgumentError OperatorView(NormNetwork(tn))
-
+        # The ket and operator layers are the networks the quadratic form was built from.
+        @test ketnetwork(qf) === tn
+        @test operatornetwork(qf) === op
+
+        # The bra layer is not stored, so it is a view sharing the ket layer's graph structure.
+        bv = branetwork(qf)
+        @test bv isa BraView
+        @test issetequal(vertices(bv), vertices(tn))
+        @test issetequal(edges(bv), edges(tn))
+        @test !issettable(bv)
+        @test !isinsertable(bv)
+        @test is_vertex_assigned(bv, 1)
         for v in vertices(tn)
-            @test kv[v] === kettensor(qf, v)
-            @test inds(ov[v]) == inds(operatortensor(qf, v))
             @test inds(bv[v]) == inds(bratensor(qf, v))
         end
     end

From 3de553adf8181f6fdd1ec6bc4b101ba96860542a Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 11:57:41 -0400
Subject: [PATCH 07/18] Remove `AbstractBilinearFormNetworkView`; define its
 forwarding methods on `BraView`

`BraView` was the only remaining subtype, so it now subtypes `AbstractITensorNetwork`
directly and carries the graph, encoding, assignment and mutability forwarding itself.
The include order places `bilinearformnetworkview.jl` first, since the forwarding methods
now name `BraView`.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/ITensorNetworksNext.jl                    |  2 +-
 .../abstractbilinearformnetworkview.jl        | 27 ++++++-------------
 src/bilinearforms/bilinearformnetworkview.jl  |  8 +++++-
 3 files changed, 16 insertions(+), 21 deletions(-)

diff --git a/src/ITensorNetworksNext.jl b/src/ITensorNetworksNext.jl
index 87c20a5a..a262eabf 100644
--- a/src/ITensorNetworksNext.jl
+++ b/src/ITensorNetworksNext.jl
@@ -15,10 +15,10 @@ include("abstracttensornetwork.jl")
 include("tensornetwork.jl")
 include("itensornetworkoperator.jl")
 include("bilinearforms/abstractbilinearformnetwork.jl")
-include("bilinearforms/abstractbilinearformnetworkview.jl")
 include("bilinearforms/normnetwork.jl")
 include("bilinearforms/quadraticformnetwork.jl")
 include("bilinearforms/bilinearformnetworkview.jl")
+include("bilinearforms/abstractbilinearformnetworkview.jl")
 include("ITensorNetworkGenerators/ITensorNetworkGenerators.jl")
 include("contract_network.jl")
 
diff --git a/src/bilinearforms/abstractbilinearformnetworkview.jl b/src/bilinearforms/abstractbilinearformnetworkview.jl
index 262bc4e7..02bc8607 100644
--- a/src/bilinearforms/abstractbilinearformnetworkview.jl
+++ b/src/bilinearforms/abstractbilinearformnetworkview.jl
@@ -3,39 +3,28 @@ using Dictionaries: Dictionaries, isinsertable, issettable
 using Graphs: Graphs, edges, vertices
 using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
 
-"""
-    abstract type AbstractBilinearFormNetworkView{T, V, I} <: AbstractITensorNetwork{T, V}
-
-Supertype of the single-layer views of an `AbstractBilinearFormNetwork{T, V, I}`.
-
-A subtype implements `Base.parent`, returning the network it views, and
-`DataGraphs.get_vertex_data`, returning that layer's tensor at a vertex. Its graph structure
-and mutability are those of the parent network.
-"""
-abstract type AbstractBilinearFormNetworkView{T, V, I} <: AbstractITensorNetwork{T, V} end
-
 # ====================================== Graphs.jl ======================================= #
 
-Graphs.edges(nnv::AbstractBilinearFormNetworkView) = edges(parent(nnv))
-Graphs.vertices(nnv::AbstractBilinearFormNetworkView) = vertices(parent(nnv))
+Graphs.edges(nnv::BraView) = edges(parent(nnv))
+Graphs.vertices(nnv::BraView) = vertices(parent(nnv))
 
 # ==================================== NamedGraphs.jl ==================================== #
 
-function NamedGraphs.encoded_vertex(nnv::AbstractBilinearFormNetworkView, vertex)
+function NamedGraphs.encoded_vertex(nnv::BraView, vertex)
     return encoded_vertex(parent(nnv), vertex)
 end
-function NamedGraphs.decoded_vertex(nnv::AbstractBilinearFormNetworkView, code::Integer)
+function NamedGraphs.decoded_vertex(nnv::BraView, code::Integer)
     return decoded_vertex(parent(nnv), code)
 end
-NamedGraphs.encoded_graph(nnv::AbstractBilinearFormNetworkView) = encoded_graph(parent(nnv))
+NamedGraphs.encoded_graph(nnv::BraView) = encoded_graph(parent(nnv))
 
 # ==================================== DataGraphs.jl ===================================== #
 
-function DataGraphs.is_vertex_assigned(nnv::AbstractBilinearFormNetworkView, vertex)
+function DataGraphs.is_vertex_assigned(nnv::BraView, vertex)
     return is_vertex_assigned(parent(nnv), vertex)
 end
 
 # =================================== Dictionaries.jl ==================================== #
 
-Dictionaries.issettable(nnv::AbstractBilinearFormNetworkView) = issettable(parent(nnv))
-Dictionaries.isinsertable(nnv::AbstractBilinearFormNetworkView) = isinsertable(parent(nnv))
+Dictionaries.issettable(nnv::BraView) = issettable(parent(nnv))
+Dictionaries.isinsertable(nnv::BraView) = isinsertable(parent(nnv))
diff --git a/src/bilinearforms/bilinearformnetworkview.jl b/src/bilinearforms/bilinearformnetworkview.jl
index d9622fee..277529a5 100644
--- a/src/bilinearforms/bilinearformnetworkview.jl
+++ b/src/bilinearforms/bilinearformnetworkview.jl
@@ -2,8 +2,14 @@ using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
 using Dictionaries: Dictionaries, isinsertable, issettable
 using Graphs: Graphs, edges, vertices
 
+"""
+    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}} <:
-    AbstractBilinearFormNetworkView{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)

From 38198a012b824f9afb1265879313be2de42840b8 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 11:59:17 -0400
Subject: [PATCH 08/18] Move `BraView` into `abstractbilinearformnetwork.jl`

Pure move: the contents of `bilinearformnetworkview.jl` and
`abstractbilinearformnetworkview.jl` are appended verbatim, and their `using` lines are
merged into the existing imports. `src/bilinearforms/` now holds the interface, `NormNetwork`
and `QuadraticFormNetwork`.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/ITensorNetworksNext.jl                    |  2 -
 .../abstractbilinearformnetwork.jl            | 51 ++++++++++++++++++-
 .../abstractbilinearformnetworkview.jl        | 30 -----------
 src/bilinearforms/bilinearformnetworkview.jl  | 23 ---------
 4 files changed, 50 insertions(+), 56 deletions(-)
 delete mode 100644 src/bilinearforms/abstractbilinearformnetworkview.jl
 delete mode 100644 src/bilinearforms/bilinearformnetworkview.jl

diff --git a/src/ITensorNetworksNext.jl b/src/ITensorNetworksNext.jl
index a262eabf..ac62144b 100644
--- a/src/ITensorNetworksNext.jl
+++ b/src/ITensorNetworksNext.jl
@@ -17,8 +17,6 @@ include("itensornetworkoperator.jl")
 include("bilinearforms/abstractbilinearformnetwork.jl")
 include("bilinearforms/normnetwork.jl")
 include("bilinearforms/quadraticformnetwork.jl")
-include("bilinearforms/bilinearformnetworkview.jl")
-include("bilinearforms/abstractbilinearformnetworkview.jl")
 include("ITensorNetworkGenerators/ITensorNetworkGenerators.jl")
 include("contract_network.jl")
 
diff --git a/src/bilinearforms/abstractbilinearformnetwork.jl b/src/bilinearforms/abstractbilinearformnetwork.jl
index 6d0dd6c2..10f54661 100644
--- a/src/bilinearforms/abstractbilinearformnetwork.jl
+++ b/src/bilinearforms/abstractbilinearformnetwork.jl
@@ -1,5 +1,8 @@
-using Dictionaries: Dictionaries
+using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
+using Dictionaries: Dictionaries, isinsertable, issettable
+using Graphs: Graphs, edges, vertices
 using ITensorBase: conj, name, rename, setname
+using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
 
 """
     abstract type AbstractBilinearFormNetwork{T, V, I} <: AbstractITensorNetwork{T, V}
@@ -78,3 +81,49 @@ The bra-layer network of `bn`. Unless a subtype stores its bra layer as a networ
 branetwork(bn::AbstractBilinearFormNetwork) = BraView(bn)
 
 indmap(bn::AbstractBilinearFormNetwork, ind) = setname(conj(ind), braname(bn, name(ind)))
+
+"""
+    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)
+
+# ====================================== 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))
+
+# ==================================== DataGraphs.jl ===================================== #
+
+function DataGraphs.is_vertex_assigned(nnv::BraView, vertex)
+    return is_vertex_assigned(parent(nnv), vertex)
+end
+
+# =================================== Dictionaries.jl ==================================== #
+
+Dictionaries.issettable(nnv::BraView) = issettable(parent(nnv))
+Dictionaries.isinsertable(nnv::BraView) = isinsertable(parent(nnv))
diff --git a/src/bilinearforms/abstractbilinearformnetworkview.jl b/src/bilinearforms/abstractbilinearformnetworkview.jl
deleted file mode 100644
index 02bc8607..00000000
--- a/src/bilinearforms/abstractbilinearformnetworkview.jl
+++ /dev/null
@@ -1,30 +0,0 @@
-using DataGraphs: DataGraphs, is_vertex_assigned
-using Dictionaries: Dictionaries, isinsertable, issettable
-using Graphs: Graphs, edges, vertices
-using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
-
-# ====================================== 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))
-
-# ==================================== DataGraphs.jl ===================================== #
-
-function DataGraphs.is_vertex_assigned(nnv::BraView, vertex)
-    return is_vertex_assigned(parent(nnv), vertex)
-end
-
-# =================================== Dictionaries.jl ==================================== #
-
-Dictionaries.issettable(nnv::BraView) = issettable(parent(nnv))
-Dictionaries.isinsertable(nnv::BraView) = isinsertable(parent(nnv))
diff --git a/src/bilinearforms/bilinearformnetworkview.jl b/src/bilinearforms/bilinearformnetworkview.jl
deleted file mode 100644
index 277529a5..00000000
--- a/src/bilinearforms/bilinearformnetworkview.jl
+++ /dev/null
@@ -1,23 +0,0 @@
-using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
-using Dictionaries: Dictionaries, isinsertable, issettable
-using Graphs: Graphs, edges, vertices
-
-"""
-    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)

From f0a4580796a316ab73ac5609c875efb5bee9e25f Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 12:56:44 -0400
Subject: [PATCH 09/18] Combine the two `DataGraphs.jl` sections of `BraView`

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/bilinearforms/abstractbilinearformnetwork.jl | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/src/bilinearforms/abstractbilinearformnetwork.jl b/src/bilinearforms/abstractbilinearformnetwork.jl
index 10f54661..6a8d8c9b 100644
--- a/src/bilinearforms/abstractbilinearformnetwork.jl
+++ b/src/bilinearforms/abstractbilinearformnetwork.jl
@@ -101,6 +101,9 @@ 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 ======================================= #
 
@@ -117,12 +120,6 @@ function NamedGraphs.decoded_vertex(nnv::BraView, code::Integer)
 end
 NamedGraphs.encoded_graph(nnv::BraView) = encoded_graph(parent(nnv))
 
-# ==================================== DataGraphs.jl ===================================== #
-
-function DataGraphs.is_vertex_assigned(nnv::BraView, vertex)
-    return is_vertex_assigned(parent(nnv), vertex)
-end
-
 # =================================== Dictionaries.jl ==================================== #
 
 Dictionaries.issettable(nnv::BraView) = issettable(parent(nnv))

From 1c0774d2088098e2df2c8e7b1cbe9b09a09760ee Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 13:35:07 -0400
Subject: [PATCH 10/18] Return an `AbstractGramian` from bilinear-form networks
 and split it in `contract_network`

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .../abstractbilinearformnetwork.jl            | 30 ++++++---
 src/bilinearforms/normnetwork.jl              | 27 +++++---
 src/bilinearforms/quadraticformnetwork.jl     | 63 +++++++++++++------
 src/contract_network.jl                       | 19 ++++++
 test/test_beliefpropagation.jl                |  4 +-
 test/test_normnetwork.jl                      | 23 ++++++-
 test/test_quadraticformnetwork.jl             | 52 +++++++++++++--
 7 files changed, 174 insertions(+), 44 deletions(-)

diff --git a/src/bilinearforms/abstractbilinearformnetwork.jl b/src/bilinearforms/abstractbilinearformnetwork.jl
index 6a8d8c9b..074b54f9 100644
--- a/src/bilinearforms/abstractbilinearformnetwork.jl
+++ b/src/bilinearforms/abstractbilinearformnetwork.jl
@@ -10,14 +10,23 @@ using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
 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 and implements [`braname`](@ref) together with one
-accessor per layer: [`kettensor`](@ref), [`bratensor`](@ref) and, where the subtype has an
-operator layer, [`operatortensor`](@ref). `bratensor` has a default built from `kettensor`
-and `braname`. The layers as whole networks are returned by [`ketnetwork`](@ref),
+A subtype supplies its own graph structure, implements [`braname`](@ref), and returns an
+[`AbstractGramian`](@ref) from `getindex`. The per-vertex accessors [`kettensor`](@ref),
+[`bratensor`](@ref) and, where the subtype has an operator layer, [`operatortensor`](@ref) read
+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) and `layertensors`; the bra tensor is built from the ket
+tensor and the name map each time it is requested.
+"""
+abstract type AbstractGramian end
+
 # =================================== Dictionaries.jl ==================================== #
 
 Dictionaries.issettable(::AbstractBilinearFormNetwork) = false
@@ -34,29 +43,34 @@ function braname end
 
 """
     kettensor(bn::AbstractBilinearFormNetwork, vertex)
+    kettensor(g::AbstractGramian)
 
 The ket-layer tensor at `vertex`.
 """
 function kettensor end
+kettensor(bn::AbstractBilinearFormNetwork, vertex) = kettensor(bn[vertex])
 
 """
     operatortensor(bn::AbstractBilinearFormNetwork, vertex)
+    operatortensor(g::AbstractGramian)
 
 The operator-layer tensor at `vertex`, 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
+operatortensor(bn::AbstractBilinearFormNetwork, vertex) = operatortensor(bn[vertex])
 
-function conj_bratensor(bn::AbstractBilinearFormNetwork, vertex)
-    return rename(n -> braname(bn, n), kettensor(bn, vertex))
-end
+conj_bratensor(g::AbstractGramian) = rename(n -> braname(g, n), kettensor(g))
+conj_bratensor(bn::AbstractBilinearFormNetwork, vertex) = conj_bratensor(bn[vertex])
 
 """
     bratensor(bn::AbstractBilinearFormNetwork, vertex)
+    bratensor(g::AbstractGramian)
 
 The bra-layer tensor at `vertex`.
 """
-bratensor(bn::AbstractBilinearFormNetwork, vertex) = conj(conj_bratensor(bn, vertex))
+bratensor(g::AbstractGramian) = conj(conj_bratensor(g))
+bratensor(bn::AbstractBilinearFormNetwork, vertex) = bratensor(bn[vertex])
 
 """
     ketnetwork(bn::AbstractBilinearFormNetwork)
diff --git a/src/bilinearforms/normnetwork.jl b/src/bilinearforms/normnetwork.jl
index 03add5a0..2a34f09f 100644
--- a/src/bilinearforms/normnetwork.jl
+++ b/src/bilinearforms/normnetwork.jl
@@ -1,5 +1,5 @@
 using Dictionaries: Dictionary
-using ITensorBase: LazyNamedTensor, lazy, similar_operator, uniquename
+using ITensorBase: similar_operator, uniquename
 using ITensorNetworksNext
 
 """
@@ -26,7 +26,22 @@ struct NormNetwork{T, V, I} <: AbstractBilinearFormNetwork{T, V, I}
     end
 end
 
-Base.eltype(::Type{<:NormNetwork{T, V, I}}) where {T, V, I} = LazyNamedTensor{I, T}
+"""
+    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
+braname(g::NormGramian, name) = get(g.braname, name, name)
+layertensors(g::NormGramian) = (; ket = kettensor(g), bra = bratensor(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)))
@@ -49,11 +64,8 @@ function DataGraphs.is_vertex_assigned(nn::NormNetwork, vertex)
     return isassigned(nn.ket, vertex)
 end
 
-function DataGraphs.get_vertex_data(nn::NormNetwork, vertex)
-    A = kettensor(nn, vertex)
-    B = conj_bratensor(nn, vertex)
-    # TODO: implement and use a lazy `conj` via `LazyNamedDimsArrays` here?
-    return lazy(A) * lazy(conj(B))
+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 ======================================= #
@@ -67,7 +79,6 @@ function braname(nn::NormNetwork, name)
     return get(nn.braname, name, name)
 end
 
-kettensor(nn::NormNetwork, vertex) = nn.ket[vertex]
 ketnetwork(nn::NormNetwork) = nn.ket
 
 """
diff --git a/src/bilinearforms/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl
index a01b5887..75fdcab2 100644
--- a/src/bilinearforms/quadraticformnetwork.jl
+++ b/src/bilinearforms/quadraticformnetwork.jl
@@ -1,6 +1,5 @@
 using Dictionaries: Dictionary
-using ITensorBase:
-    LazyNamedTensor, inputnames, lazy, names, outputnames, rename, state, uniquename
+using ITensorBase: inputnames, names, outputnames, rename, state, uniquename
 using ITensorNetworksNext
 
 """
@@ -28,6 +27,9 @@ struct QuadraticFormNetwork{T, V, I, O <: ITensorNetworkOperator} <:
         if !issetequal(vertices(operator), vertices(ket))
             error("the operator layer must be defined on every vertex of the ket layer.")
         end
+        if !issubset(inputnames(operator), keys(ket.dimname_vertices))
+            error("every operator input name must be an index name of the ket layer.")
+        end
         acted = Set{I}(inputnames(operator))
         braname = Dictionary{I, I}()
         for (name, vertices) in pairs(ket.dimname_vertices)
@@ -39,7 +41,32 @@ struct QuadraticFormNetwork{T, V, I, O <: ITensorNetworkOperator} <:
     end
 end
 
-Base.eltype(::Type{<:QuadraticFormNetwork{T, V, I}}) where {T, V, I} = LazyNamedTensor{I, T}
+"""
+    struct QuadraticFormGramian{T, O, I} <: AbstractGramian
+
+The layers of a `QuadraticFormNetwork` at one vertex: the ket tensor, the operator network's tensor,
+the ket→bra name map and the operator's name pairing, from which the bra tensor and the renamed
+operator tensor are built when requested.
+"""
+struct QuadraticFormGramian{T, O, I} <: AbstractGramian
+    ket::T
+    operator::O
+    braname::Dictionary{I, I}
+    outputnames::Vector{I}
+    inputnames::Vector{I}
+end
+
+kettensor(g::QuadraticFormGramian) = g.ket
+braname(g::QuadraticFormGramian, name) = get(g.braname, name, name)
+function layertensors(g::QuadraticFormGramian)
+    return (; ket = kettensor(g), operator = operatortensor(g), bra = bratensor(g))
+end
+
+function Base.eltype(
+        ::Type{<:QuadraticFormNetwork{T, V, I, O}}
+    ) where {T, V, I, TO, O <: ITensorNetworkOperator{TO}}
+    return QuadraticFormGramian{T, TO, I}
+end
 
 function QuadraticFormNetwork(ket::ITensorNetwork, operator::ITensorNetworkOperator)
     return QuadraticFormNetwork(ket, operator, map(uniquename, keys(ket.dimname_vertices)))
@@ -62,12 +89,13 @@ NamedGraphs.encoded_graph(qf::QuadraticFormNetwork) = encoded_graph(qf.ket)
 
 # ==================================== DataGraphs.jl ===================================== #
 
-function DataGraphs.get_vertex_data(qf::QuadraticFormNetwork, vertex)
-    A = kettensor(qf, vertex)
-    O = operatortensor(qf, vertex)
-    B = conj_bratensor(qf, vertex)
-    # TODO: implement and use a lazy `conj` via `LazyNamedDimsArrays` here?
-    return lazy(A) * lazy(O) * lazy(conj(B))
+function DataGraphs.get_vertex_data(
+        qf::QuadraticFormNetwork{T, V, I, O}, vertex
+    ) where {T, V, I, TO, O <: ITensorNetworkOperator{TO}}
+    return QuadraticFormGramian{T, TO, I}(
+        qf.ket[vertex], state(qf.operator)[vertex], qf.braname,
+        outputnames(qf.operator), inputnames(qf.operator)
+    )
 end
 
 function DataGraphs.is_vertex_assigned(qf::QuadraticFormNetwork, vertex)
@@ -85,23 +113,20 @@ function braname(qf::QuadraticFormNetwork, name)
     return get(qf.braname, name, name)
 end
 
-kettensor(qf::QuadraticFormNetwork, vertex) = qf.ket[vertex]
 ketnetwork(qf::QuadraticFormNetwork) = qf.ket
 operatornetwork(qf::QuadraticFormNetwork) = qf.operator
 
 # Each output name is renamed to the bra name of the input name it is paired with, so the
 # operator's output legs meet the bra layer and its input legs meet the ket layer. The pairing
-# is read from the operator network rather than from `qf.operator[vertex]`, whose wrapper drops
-# a pair whose input sits on another vertex.
-function operatortensor(qf::QuadraticFormNetwork, vertex)
-    tensor = state(qf.operator)[vertex]
-    tensor_names = names(tensor)
+# is read from the operator network rather than from the operator network's per-vertex wrapper,
+# which drops a pair whose input sits on another vertex.
+function operatortensor(g::QuadraticFormGramian)
+    tensor_names = names(g.operator)
     replacements = [
-        output => braname(qf, input) for
-            (output, input) in zip(outputnames(qf.operator), inputnames(qf.operator))
-            if output in tensor_names
+        output => braname(g, input) for
+            (output, input) in zip(g.outputnames, g.inputnames) if output in tensor_names
     ]
-    return rename(tensor, replacements...)
+    return rename(g.operator, replacements...)
 end
 
 """
diff --git a/src/contract_network.jl b/src/contract_network.jl
index f13b36c8..73402a46 100644
--- a/src/contract_network.jl
+++ b/src/contract_network.jl
@@ -1,5 +1,6 @@
 using Base.Broadcast: materialize
 using Base: @kwdef
+using Dictionaries: Dictionary
 using ITensorBase: EvaluationOrderAlgorithm, Greedy, Mul, lazy, optimize_evaluation_order,
     substitute, symnamedtensor
 
@@ -35,7 +36,25 @@ end
 # operator, so mixing operators and plain tensors is the common case) widens the symbolic `Mul`
 # container to a `UnionAll` it cannot construct. `promote_type`/`convert` keep an all-plain network
 # at the plain type (the promotion is a no-op), so its fast path is unchanged.
+# A Gramian enters the contraction as its separate layer tensors, so the contraction order can
+# absorb other operands into one layer before the layers are joined.
+function split_gramians(tn)
+    any(t -> t isa AbstractGramian, tn) || return tn
+    operands = Dictionary{Any, Any}()
+    for (key, t) in pairs(tn)
+        if t isa AbstractGramian
+            for (layer, tensor) in pairs(layertensors(t))
+                insert!(operands, (key, layer), tensor)
+            end
+        else
+            insert!(operands, key, t)
+        end
+    end
+    return operands
+end
+
 function contract_network(alg::Exact, tn)
+    tn = split_gramians(tn)
     order = get_order(alg, tn)
     T = mapreduce(typeof, promote_type, tn)
     syms_to_ts = Dict(
diff --git a/test/test_beliefpropagation.jl b/test/test_beliefpropagation.jl
index f838f5d3..0b413602 100644
--- a/test/test_beliefpropagation.jl
+++ b/test/test_beliefpropagation.jl
@@ -371,7 +371,7 @@ end
 
             # A doubled vertex splits into its two layers, and the split is faithful.
             @test length(factor_tensors(nn, v)) == 2
-            @test prod(factor_tensors(nn, v)) ≈ materialize(nn[v])
+            @test prod(factor_tensors(nn, v)) ≈ contract_network([nn[v]])
             # A single-layer network's factor is a single operand.
             @test factor_tensors(network, v) == [network[v]]
 
@@ -388,7 +388,7 @@ end
             message = updated_message(algorithm, cache, nn, edge)
             # `v` has degree 4, so 3 incoming messages plus the ket and bra layers.
             @test only(counts) == 5
-            @test message ≈ contract_network([messages; [nn[v]]])
+            @test message ≈ contract_network([messages; [contract_network([nn[v]])]])
         end
     end
 end
diff --git a/test/test_normnetwork.jl b/test/test_normnetwork.jl
index 6f2b9ff3..b0b6256c 100644
--- a/test/test_normnetwork.jl
+++ b/test/test_normnetwork.jl
@@ -3,9 +3,9 @@ using Dictionaries: isinsertable, issettable
 using Graphs: edges, vertices
 using ITensorBase:
     ITensor, Index, IndexName, LazyITensor, conj, inds, name, setname, uniquename
-using ITensorNetworksNext: BraView, ITensorNetwork, NormNetwork, braname, branetwork,
-    bratensor, conj_bratensor, contract_network, indmap, ketnetwork, kettensor, normnetwork,
-    tensornetwork
+using ITensorNetworksNext: ITensorNetworksNext, BraView, ITensorNetwork, NormGramian,
+    NormNetwork, braname, branetwork, bratensor, conj_bratensor, contract_network, indmap,
+    ketnetwork, kettensor, normnetwork, tensornetwork
 using LinearAlgebra: norm
 using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph
 using Test: @test, @test_throws, @testset
@@ -119,6 +119,23 @@ end
         @test !isinsertable(bv)
     end
 
+    @testset "`NormGramian`" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        nn = NormNetwork(tn)
+        gram = nn[2]
+
+        @test gram isa NormGramian
+        @test eltype(nn) === typeof(gram)
+        # The Gramian holds the network's ket tensor and name map, not copies.
+        @test kettensor(gram) === tn[2]
+        @test gram.braname === nn.braname
+        @test inds(bratensor(gram)) == inds(conj_bratensor(gram))
+        @test keys(ITensorNetworksNext.layertensors(gram)) == (:ket, :bra)
+        # Contracting a Gramian contracts its layers.
+        @test contract_network([gram]) ≈ kettensor(gram) * bratensor(gram)
+    end
+
     @testset "contraction / physics" begin
         @testset "single normalized tensor contracts to 1" begin
             s = Index(3)
diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl
index fd18c2e6..b5203647 100644
--- a/test/test_quadraticformnetwork.jl
+++ b/test/test_quadraticformnetwork.jl
@@ -2,10 +2,11 @@ using DataGraphs: is_vertex_assigned
 using Dictionaries: isinsertable, issettable
 using Graphs: edges, vertices
 using ITensorBase: ITensor, Index, IndexName, conj, inds, inputnames, name, names, operator,
-    outputnames, rename, setname, uniquename
-using ITensorNetworksNext: BraView, ITensorNetwork, NormNetwork, QuadraticFormNetwork,
-    braname, branetwork, bratensor, conj_bratensor, contract_network, indmap, ketnetwork,
-    kettensor, operatornetwork, operatortensor, quadraticformnetwork, tensornetwork
+    outputnames, rename, setname, state, uniquename
+using ITensorNetworksNext: ITensorNetworksNext, BraView, ITensorNetwork, NormNetwork,
+    QuadraticFormGramian, QuadraticFormNetwork, braname, branetwork, bratensor,
+    conj_bratensor, contract_network, indmap, ketnetwork, kettensor, operatornetwork,
+    operatortensor, quadraticformnetwork, tensornetwork
 using LinearAlgebra: I, norm
 using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph
 using Test: @test, @test_throws, @testset
@@ -67,6 +68,49 @@ identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s;
         @test_throws ErrorException QuadraticFormNetwork(tn, identity_operator(g4, s4))
     end
 
+    @testset "`QuadraticFormGramian`" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        op = identity_operator(g, s)
+        qf = QuadraticFormNetwork(tn, op)
+        gram = qf[2]
+
+        @test gram isa QuadraticFormGramian
+        @test eltype(qf) === typeof(gram)
+        @test kettensor(gram) === tn[2]
+        @test gram.operator === state(op)[2]
+        @test keys(ITensorNetworksNext.layertensors(gram)) == (:ket, :operator, :bra)
+        @test contract_network([gram]) ≈
+            kettensor(gram) * operatortensor(gram) * bratensor(gram)
+    end
+
+    @testset "crossing operator pair" begin
+        g = named_path_graph(2)
+        tn, l, s = random_state(Float64, g)
+        # Each output is paired with the input on the other vertex.
+        out1, out2 = Index(2), Index(2)
+        optn = ITensorNetwork(Dict(1 => randn((out1, s[1])), 2 => randn((out2, s[2]))))
+        op = operator(optn, [name(out1), name(out2)], [name(s[2]), name(s[1])])
+        qf = QuadraticFormNetwork(tn, op)
+
+        @test braname(qf, name(s[2])) in names(operatortensor(qf[1]))
+        @test braname(qf, name(s[1])) in names(operatortensor(qf[2]))
+
+        psi = prod(tn)
+        bra = rename(conj(psi), name(s[2]) => name(out1), name(s[1]) => name(out2))
+        @test contract_network(qf)[] ≈ (bra * (optn[1] * optn[2] * psi))[]
+    end
+
+    @testset "operator input outside the ket network" begin
+        g = named_path_graph(2)
+        tn, l, s = random_state(Float64, g)
+        stray = Index(2)
+        out1, out2 = Index(2), Index(2)
+        optn = ITensorNetwork(Dict(1 => randn((out1, stray)), 2 => randn((out2, s[2]))))
+        op = operator(optn, [name(out1), name(out2)], [name(stray), name(s[2])])
+        @test_throws ErrorException QuadraticFormNetwork(tn, op)
+    end
+
     @testset "layer tensors and the name map" begin
         g = named_path_graph(3)
         tn, l, s = random_state(Float64, g)

From 5d224574544c82dc1681b7980b095925694aab97 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 13:41:09 -0400
Subject: [PATCH 11/18] Key every operand of a split collection as a tuple

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/contract_network.jl        | 29 ++++++++++++-----------------
 test/test_beliefpropagation.jl |  2 +-
 2 files changed, 13 insertions(+), 18 deletions(-)

diff --git a/src/contract_network.jl b/src/contract_network.jl
index 73402a46..2f2d807d 100644
--- a/src/contract_network.jl
+++ b/src/contract_network.jl
@@ -31,28 +31,23 @@ function get_order(alg::Exact, tn)
         Dict(symnamedtensor(i) => symnamedtensor(i, Tuple(axes(t))) for (i, t) in pairs(tn))
     return substitute(order, subs)
 end
-# Promote the operands to their common type before lowering to the lazy expression, so every lazy
-# operand shares one concrete type. Otherwise a network of mixed types (a plain tensor is a trivial
-# operator, so mixing operators and plain tensors is the common case) widens the symbolic `Mul`
-# container to a `UnionAll` it cannot construct. `promote_type`/`convert` keep an all-plain network
-# at the plain type (the promotion is a no-op), so its fast path is unchanged.
 # A Gramian enters the contraction as its separate layer tensors, so the contraction order can
-# absorb other operands into one layer before the layers are joined.
+# absorb other operands into one layer before the layers are joined. Every operand is keyed by a
+# `(key, layer)` tuple so all keys share one concrete type even when only some entries split.
 function split_gramians(tn)
     any(t -> t isa AbstractGramian, tn) || return tn
-    operands = Dictionary{Any, Any}()
-    for (key, t) in pairs(tn)
-        if t isa AbstractGramian
-            for (layer, tensor) in pairs(layertensors(t))
-                insert!(operands, (key, layer), tensor)
-            end
-        else
-            insert!(operands, key, t)
-        end
-    end
-    return operands
+    pairs_split = [
+        (key, layer) => tensor for (key, t) in pairs(tn) for
+            (layer, tensor) in (t isa AbstractGramian ? pairs(layertensors(t)) : [:tensor => t])
+    ]
+    return Dictionary(first.(pairs_split), last.(pairs_split))
 end
 
+# Promote the operands to their common type before lowering to the lazy expression, so every lazy
+# operand shares one concrete type. Otherwise a network of mixed types (a plain tensor is a trivial
+# operator, so mixing operators and plain tensors is the common case) widens the symbolic `Mul`
+# container to a `UnionAll` it cannot construct. `promote_type`/`convert` keep an all-plain network
+# at the plain type (the promotion is a no-op), so its fast path is unchanged.
 function contract_network(alg::Exact, tn)
     tn = split_gramians(tn)
     order = get_order(alg, tn)
diff --git a/test/test_beliefpropagation.jl b/test/test_beliefpropagation.jl
index 0b413602..04760408 100644
--- a/test/test_beliefpropagation.jl
+++ b/test/test_beliefpropagation.jl
@@ -388,7 +388,7 @@ end
             message = updated_message(algorithm, cache, nn, edge)
             # `v` has degree 4, so 3 incoming messages plus the ket and bra layers.
             @test only(counts) == 5
-            @test message ≈ contract_network([messages; [contract_network([nn[v]])]])
+            @test message ≈ contract_network([messages; [nn[v]]])
         end
     end
 end

From 329e1c2aa2113d8d18859b77ecce094d9f7dc946 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 13:46:12 -0400
Subject: [PATCH 12/18] Shorten the `split_gramians` comment

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/contract_network.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/contract_network.jl b/src/contract_network.jl
index 2f2d807d..bc13854c 100644
--- a/src/contract_network.jl
+++ b/src/contract_network.jl
@@ -31,9 +31,8 @@ function get_order(alg::Exact, tn)
         Dict(symnamedtensor(i) => symnamedtensor(i, Tuple(axes(t))) for (i, t) in pairs(tn))
     return substitute(order, subs)
 end
-# A Gramian enters the contraction as its separate layer tensors, so the contraction order can
-# absorb other operands into one layer before the layers are joined. Every operand is keyed by a
-# `(key, layer)` tuple so all keys share one concrete type even when only some entries split.
+# A Gramian enters as its separate layer tensors, so the order can place other operands between
+# layers; every operand gets a `(key, layer)` key so all keys share one concrete type.
 function split_gramians(tn)
     any(t -> t isa AbstractGramian, tn) || return tn
     pairs_split = [

From 15a61630c5dcb4295a283879ce1566511c58194b Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 13:56:13 -0400
Subject: [PATCH 13/18] Delete `factor_tensors`; `contract_network` splits
 Gramians itself

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/beliefpropagation/beliefpropagation.jl | 15 ++-------------
 src/beliefpropagation/messagecache.jl      |  5 +----
 test/test_beliefpropagation.jl             | 18 ++++++++----------
 3 files changed, 11 insertions(+), 27 deletions(-)

diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl
index 7dba8aa2..75f118f4 100644
--- a/src/beliefpropagation/beliefpropagation.jl
+++ b/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
 
diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index 96e7f417..be902628 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/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))
diff --git a/test/test_beliefpropagation.jl b/test/test_beliefpropagation.jl
index 04760408..0d034e9d 100644
--- a/test/test_beliefpropagation.jl
+++ b/test/test_beliefpropagation.jl
@@ -1,5 +1,4 @@
 import AlgorithmsInterface as AI
-using Base.Broadcast: materialize
 using DataGraphs: DataGraphs, DataGraph, edge_data, edge_data_type
 using Dictionaries: Dictionary, dictionary, set!
 using GradedArrays: U1, gradedrange, isdual
@@ -10,9 +9,9 @@ using ITensorBase:
 using ITensorNetworksNext: ITensorNetworksNext, Exact, ITensorNetwork, MessageCache,
     NormNetwork, SimpleMessageUpdate, StopWhenConverged, beliefpropagation,
     bethe_free_energy, bethe_free_entropy, bratensor, contract_network, contraction_order,
-    edge_scalar, edge_scalars, factor_tensors, incoming_messages, insertlink!, kettensor,
-    linkaxes, linkinds, message_environment, messagecache, region_scalar, subgraph,
-    tensornetwork, updated_message, vertex_scalar, vertex_scalars
+    edge_scalar, edge_scalars, incoming_messages, insertlink!, kettensor, linkaxes,
+    linkinds, message_environment, messagecache, region_scalar, subgraph, tensornetwork,
+    updated_message, vertex_scalar, vertex_scalars
 using LinearAlgebra: LinearAlgebra, norm, tr
 using NamedGraphs: NamedEdge, all_edges, incident_edges, named_comb_tree, named_grid,
     named_path_graph, vertextype
@@ -369,11 +368,10 @@ end
             nn = NormNetwork(network)
             v = (2, 2)
 
-            # A doubled vertex splits into its two layers, and the split is faithful.
-            @test length(factor_tensors(nn, v)) == 2
-            @test prod(factor_tensors(nn, v)) ≈ contract_network([nn[v]])
-            # A single-layer network's factor is a single operand.
-            @test factor_tensors(network, v) == [network[v]]
+            # A doubled vertex is a Gramian whose layers contract to the multiplied-out vertex.
+            gram = nn[v]
+            doubled = kettensor(gram) * bratensor(gram)
+            @test contract_network([gram]) ≈ doubled
 
             # The message update passes the layers to `contract_network` as separate operands, so
             # the contraction order can interleave the incoming messages between them, and the
@@ -388,7 +386,7 @@ end
             message = updated_message(algorithm, cache, nn, edge)
             # `v` has degree 4, so 3 incoming messages plus the ket and bra layers.
             @test only(counts) == 5
-            @test message ≈ contract_network([messages; [nn[v]]])
+            @test message ≈ contract_network([messages; [doubled]])
         end
     end
 end

From 965df6b42dc02c8d0ccd11b96bf1b4e286ac687a Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 14:13:44 -0400
Subject: [PATCH 14/18] Split Gramians in `contraction_order`; test the
 pass-through; tidy docs

`contraction_order(tn; alg)` now splits Gramians before dispatching, matching
`contract_network(alg::Exact, tn)`; previously `contraction_order(nn)` on a
`NormNetwork` threw a `MethodError` calling `size` on a Gramian. Adds tests for
that fix and for `split_gramians`' pass-through on inputs with no Gramian.
Also extends the `braname` docstring to cover its `AbstractGramian` method and
shortens the `operatortensor(g::QuadraticFormGramian)` comment.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .../abstractbilinearformnetwork.jl              |  4 +++-
 src/bilinearforms/quadraticformnetwork.jl       |  6 ++----
 src/contract_network.jl                         |  2 +-
 test/test_contract_network.jl                   | 11 +++++++++--
 test/test_normnetwork.jl                        | 17 ++++++++++++++---
 5 files changed, 29 insertions(+), 11 deletions(-)

diff --git a/src/bilinearforms/abstractbilinearformnetwork.jl b/src/bilinearforms/abstractbilinearformnetwork.jl
index 074b54f9..244b13ad 100644
--- a/src/bilinearforms/abstractbilinearformnetwork.jl
+++ b/src/bilinearforms/abstractbilinearformnetwork.jl
@@ -36,8 +36,10 @@ Dictionaries.isinsertable(::AbstractBilinearFormNetwork) = false
 
 """
     braname(bn::AbstractBilinearFormNetwork, name)
+    braname(g::AbstractGramian, name)
 
-The bra-layer index name corresponding to the ket-layer index name `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
 
diff --git a/src/bilinearforms/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl
index 75fdcab2..91994b9a 100644
--- a/src/bilinearforms/quadraticformnetwork.jl
+++ b/src/bilinearforms/quadraticformnetwork.jl
@@ -116,10 +116,8 @@ end
 ketnetwork(qf::QuadraticFormNetwork) = qf.ket
 operatornetwork(qf::QuadraticFormNetwork) = qf.operator
 
-# Each output name is renamed to the bra name of the input name it is paired with, so the
-# operator's output legs meet the bra layer and its input legs meet the ket layer. The pairing
-# is read from the operator network rather than from the operator network's per-vertex wrapper,
-# which drops a pair whose input sits on another vertex.
+# Each output name is renamed to the bra name of the input name it is paired with, using the
+# whole operator's pairing: the per-vertex wrapper drops a pair whose input sits on another vertex.
 function operatortensor(g::QuadraticFormGramian)
     tensor_names = names(g.operator)
     replacements = [
diff --git a/src/contract_network.jl b/src/contract_network.jl
index bc13854c..c4897f3f 100644
--- a/src/contract_network.jl
+++ b/src/contract_network.jl
@@ -61,7 +61,7 @@ end
 # `contraction_order`
 function contraction_order end
 function contraction_order(tn; alg = Greedy())
-    return contraction_order(alg, tn)
+    return contraction_order(alg, split_gramians(tn))
 end
 # Convert the tensor network to a flat symbolic multiplication expression.
 struct Flat end
diff --git a/test/test_contract_network.jl b/test/test_contract_network.jl
index 667b7d54..c5dd719f 100644
--- a/test/test_contract_network.jl
+++ b/test/test_contract_network.jl
@@ -1,8 +1,8 @@
 using Graphs: edges, vertices
 using ITensorBase:
     Greedy, Index, NamedTensorOperator, inputnames, operator, outputnames, state
-using ITensorNetworksNext: Exact, ITensorNetwork, LeftAssociative, contract_network,
-    linkinds, siteinds, tensornetwork
+using ITensorNetworksNext: ITensorNetworksNext, Exact, ITensorNetwork, LeftAssociative,
+    contract_network, linkinds, siteinds, tensornetwork
 using NamedGraphs: incident_edges, named_grid
 using OMEinsumContractionOrders: ExhaustiveSearch, GreedyMethod, TreeSA
 using Test: @test, @testset
@@ -17,6 +17,10 @@ using Test: @test, @testset
         C = [5.0, 1.0][j]
         D = [-2.0, 3.0, 4.0, 5.0, 1.0][k]
 
+        # A `Vector` of tensors holding no Gramian is returned unchanged.
+        ABCD = [A, B, C, D]
+        @test ITensorNetworksNext.split_gramians(ABCD) === ABCD
+
         ABCD_1 = contract_network([A, B, C, D]; alg = orderalg(LeftAssociative()))
         ABCD_2 = contract_network([A, B, C, D]; alg = orderalg(Greedy()))
         ABCD_3 = contract_network([A, B, C, D]; alg = orderalg(ExhaustiveSearch()))
@@ -37,6 +41,9 @@ using Test: @test, @testset
             return randn(Tuple(is))
         end
 
+        # A plain `ITensorNetwork` holds no Gramian, so `split_gramians` returns it unchanged.
+        @test ITensorNetworksNext.split_gramians(tn) === tn
+
         z1 = contract_network(tn; alg = orderalg(LeftAssociative()))[]
         z2 = contract_network(tn; alg = orderalg(Greedy()))[]
         z3 = contract_network(tn; alg = orderalg(ExhaustiveSearch()))[]
diff --git a/test/test_normnetwork.jl b/test/test_normnetwork.jl
index b0b6256c..441e2ed0 100644
--- a/test/test_normnetwork.jl
+++ b/test/test_normnetwork.jl
@@ -3,9 +3,9 @@ using Dictionaries: isinsertable, issettable
 using Graphs: edges, vertices
 using ITensorBase:
     ITensor, Index, IndexName, LazyITensor, conj, inds, name, setname, uniquename
-using ITensorNetworksNext: ITensorNetworksNext, BraView, ITensorNetwork, NormGramian,
-    NormNetwork, braname, branetwork, bratensor, conj_bratensor, contract_network, indmap,
-    ketnetwork, kettensor, normnetwork, tensornetwork
+using ITensorNetworksNext: ITensorNetworksNext, BraView, Exact, ITensorNetwork, NormGramian,
+    NormNetwork, braname, branetwork, bratensor, conj_bratensor, contract_network,
+    contraction_order, indmap, ketnetwork, kettensor, normnetwork, tensornetwork
 using LinearAlgebra: norm
 using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph
 using Test: @test, @test_throws, @testset
@@ -136,6 +136,17 @@ end
         @test contract_network([gram]) ≈ kettensor(gram) * bratensor(gram)
     end
 
+    @testset "`contraction_order` on a `NormNetwork`" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        nn = NormNetwork(tn)
+
+        # `contraction_order` splits the Gramians before computing an order, so it does not
+        # throw trying to call `size` on a `NormGramian`.
+        order = contraction_order(nn)
+        @test contract_network(nn; alg = Exact(; order))[] ≈ contract_network(nn)[]
+    end
+
     @testset "contraction / physics" begin
         @testset "single normalized tensor contracts to 1" begin
             s = Index(3)

From 3222934ca5068cae229ff91a8861ffe0700d51ae Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 14:47:13 -0400
Subject: [PATCH 15/18] Add `inds`, `names` and `axes` to Gramians

A Gramian's indices are the layer indices no other layer at that vertex shares, so the
network-level `linkinds`, `linknames`, `linkaxes`, `siteinds`, `sitenames`, `siteaxes` and
`dimnamevertices` work on bilinear-form networks again. The bra indices are read from
`conj_bratensor`, so no tensor data is conjugated.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .../abstractbilinearformnetwork.jl            | 18 ++++++++++---
 src/bilinearforms/normnetwork.jl              |  1 +
 src/bilinearforms/quadraticformnetwork.jl     |  3 +++
 test/test_normnetwork.jl                      | 27 +++++++++++++++++--
 test/test_quadraticformnetwork.jl             | 13 +++++++++
 5 files changed, 57 insertions(+), 5 deletions(-)

diff --git a/src/bilinearforms/abstractbilinearformnetwork.jl b/src/bilinearforms/abstractbilinearformnetwork.jl
index 244b13ad..395b809f 100644
--- a/src/bilinearforms/abstractbilinearformnetwork.jl
+++ b/src/bilinearforms/abstractbilinearformnetwork.jl
@@ -1,7 +1,7 @@
 using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
 using Dictionaries: Dictionaries, isinsertable, issettable
 using Graphs: Graphs, edges, vertices
-using ITensorBase: conj, name, rename, setname
+using ITensorBase: ITensorBase, conj, inds, name, rename, setname
 using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
 
 """
@@ -22,8 +22,9 @@ abstract type AbstractBilinearFormNetwork{T, V, I} <: AbstractITensorNetwork{T,
     abstract type AbstractGramian
 
 The layers of an `AbstractBilinearFormNetwork` at one vertex. A subtype implements
-[`kettensor`](@ref), [`braname`](@ref) and `layertensors`; the bra tensor is built from the ket
-tensor and the name map each time it is requested.
+[`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
 
@@ -74,6 +75,17 @@ The bra-layer tensor at `vertex`.
 bratensor(g::AbstractGramian) = conj(conj_bratensor(g))
 bratensor(bn::AbstractBilinearFormNetwork, vertex) = bratensor(bn[vertex])
 
+# 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)
 
diff --git a/src/bilinearforms/normnetwork.jl b/src/bilinearforms/normnetwork.jl
index 2a34f09f..868931c0 100644
--- a/src/bilinearforms/normnetwork.jl
+++ b/src/bilinearforms/normnetwork.jl
@@ -40,6 +40,7 @@ end
 kettensor(g::NormGramian) = g.ket
 braname(g::NormGramian, name) = get(g.braname, name, name)
 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}
 
diff --git a/src/bilinearforms/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl
index 91994b9a..a43af479 100644
--- a/src/bilinearforms/quadraticformnetwork.jl
+++ b/src/bilinearforms/quadraticformnetwork.jl
@@ -61,6 +61,9 @@ braname(g::QuadraticFormGramian, name) = get(g.braname, name, name)
 function layertensors(g::QuadraticFormGramian)
     return (; ket = kettensor(g), operator = operatortensor(g), bra = bratensor(g))
 end
+function layerinds(g::QuadraticFormGramian)
+    return (inds(kettensor(g)), inds(operatortensor(g)), brainds(g))
+end
 
 function Base.eltype(
         ::Type{<:QuadraticFormNetwork{T, V, I, O}}
diff --git a/test/test_normnetwork.jl b/test/test_normnetwork.jl
index 441e2ed0..65c110e0 100644
--- a/test/test_normnetwork.jl
+++ b/test/test_normnetwork.jl
@@ -2,10 +2,11 @@ using DataGraphs: is_vertex_assigned
 using Dictionaries: isinsertable, issettable
 using Graphs: edges, vertices
 using ITensorBase:
-    ITensor, Index, IndexName, LazyITensor, conj, inds, name, setname, uniquename
+    ITensor, Index, IndexName, LazyITensor, conj, inds, name, names, setname, uniquename
 using ITensorNetworksNext: ITensorNetworksNext, BraView, Exact, ITensorNetwork, NormGramian,
     NormNetwork, braname, branetwork, bratensor, conj_bratensor, contract_network,
-    contraction_order, indmap, ketnetwork, kettensor, normnetwork, tensornetwork
+    contraction_order, dimnamevertices, indmap, ketnetwork, kettensor, linkaxes, linkinds,
+    linknames, normnetwork, siteaxes, siteinds, sitenames, tensornetwork
 using LinearAlgebra: norm
 using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph
 using Test: @test, @test_throws, @testset
@@ -134,6 +135,28 @@ end
         @test keys(ITensorNetworksNext.layertensors(gram)) == (:ket, :bra)
         # Contracting a Gramian contracts its layers.
         @test contract_network([gram]) ≈ kettensor(gram) * bratensor(gram)
+        # A Gramian's indices are those its ket and bra layers do not share.
+        @test issetequal(inds(gram), inds(kettensor(gram) * bratensor(gram)))
+        @test names(gram) == name.(inds(gram))
+        @test axes(gram) == Tuple(inds(gram))
+    end
+
+    @testset "index queries on a `NormNetwork`" begin
+        g = named_path_graph(3)
+        tn, l, s = random_state(Float64, g)
+        nn = NormNetwork(tn)
+        e = NamedEdge(1 => 2)
+        lname = name(l[e])
+
+        # A link of the norm network is the ket link together with its bra-layer copy.
+        @test issetequal(linknames(nn, e), [lname, braname(nn, lname)])
+        @test issetequal(name.(linkinds(nn, e)), linknames(nn, e))
+        @test issetequal(name.(linkaxes(nn, e)), linknames(nn, e))
+        @test issetequal(dimnamevertices(nn, lname), [1, 2])
+        # The site index contracts between the layers, so no vertex has a site index.
+        @test isempty(siteinds(nn, 2))
+        @test isempty(sitenames(nn, 2))
+        @test isempty(siteaxes(nn, 2))
     end
 
     @testset "`contraction_order` on a `NormNetwork`" begin
diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl
index b5203647..726de774 100644
--- a/test/test_quadraticformnetwork.jl
+++ b/test/test_quadraticformnetwork.jl
@@ -82,6 +82,12 @@ identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s;
         @test keys(ITensorNetworksNext.layertensors(gram)) == (:ket, :operator, :bra)
         @test contract_network([gram]) ≈
             kettensor(gram) * operatortensor(gram) * bratensor(gram)
+        # A Gramian's indices are those no two of its layers share.
+        @test issetequal(
+            inds(gram), inds(kettensor(gram) * operatortensor(gram) * bratensor(gram))
+        )
+        @test names(gram) == name.(inds(gram))
+        @test axes(gram) == Tuple(inds(gram))
     end
 
     @testset "crossing operator pair" begin
@@ -95,6 +101,13 @@ identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s;
 
         @test braname(qf, name(s[2])) in names(operatortensor(qf[1]))
         @test braname(qf, name(s[1])) in names(operatortensor(qf[2]))
+        # Both bra site names dangle at vertex 1: one from the operator, one from the bra.
+        @test braname(qf, name(s[2])) in names(qf[1])
+        @test braname(qf, name(s[1])) in names(qf[1])
+        gram = qf[1]
+        @test issetequal(
+            inds(gram), inds(kettensor(gram) * operatortensor(gram) * bratensor(gram))
+        )
 
         psi = prod(tn)
         bra = rename(conj(psi), name(s[2]) => name(out1), name(s[1]) => name(out2))

From c197677e61e230c47ce2d121363db6c051fd40e3 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 15:03:05 -0400
Subject: [PATCH 16/18] Store the vertex operator in `QuadraticFormGramian`;
 reject crossing operator pairs

`ITensorNetworkOperator` now throws when a pair's output and input sit on different vertices,
so each vertex operator carries its full pairing. `QuadraticFormGramian` stores that operator
in place of the raw tensor and the network-wide name vectors.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/bilinearforms/quadraticformnetwork.jl | 33 ++++++++++-------------
 src/itensornetworkoperator.jl             | 20 +++++++++-----
 test/test_itensornetworkoperator.jl       | 19 +++----------
 test/test_quadraticformnetwork.jl         | 29 +++-----------------
 4 files changed, 36 insertions(+), 65 deletions(-)

diff --git a/src/bilinearforms/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl
index a43af479..b0ff57f8 100644
--- a/src/bilinearforms/quadraticformnetwork.jl
+++ b/src/bilinearforms/quadraticformnetwork.jl
@@ -1,5 +1,5 @@
 using Dictionaries: Dictionary
-using ITensorBase: inputnames, names, outputnames, rename, state, uniquename
+using ITensorBase: inputnames, outputnames, rename, state, uniquename
 using ITensorNetworksNext
 
 """
@@ -44,16 +44,14 @@ end
 """
     struct QuadraticFormGramian{T, O, I} <: AbstractGramian
 
-The layers of a `QuadraticFormNetwork` at one vertex: the ket tensor, the operator network's tensor,
-the ket→bra name map and the operator's name pairing, from which the bra tensor and the renamed
-operator tensor are built when requested.
+The layers of a `QuadraticFormNetwork` at one vertex: the ket tensor, the operator at that vertex
+and the ket→bra name map, from which the bra tensor and the renamed operator tensor are built
+when requested.
 """
 struct QuadraticFormGramian{T, O, I} <: AbstractGramian
     ket::T
     operator::O
     braname::Dictionary{I, I}
-    outputnames::Vector{I}
-    inputnames::Vector{I}
 end
 
 kettensor(g::QuadraticFormGramian) = g.ket
@@ -65,10 +63,8 @@ function layerinds(g::QuadraticFormGramian)
     return (inds(kettensor(g)), inds(operatortensor(g)), brainds(g))
 end
 
-function Base.eltype(
-        ::Type{<:QuadraticFormNetwork{T, V, I, O}}
-    ) where {T, V, I, TO, O <: ITensorNetworkOperator{TO}}
-    return QuadraticFormGramian{T, TO, I}
+function Base.eltype(::Type{<:QuadraticFormNetwork{T, V, I, O}}) where {T, V, I, O}
+    return QuadraticFormGramian{T, eltype(O), I}
 end
 
 function QuadraticFormNetwork(ket::ITensorNetwork, operator::ITensorNetworkOperator)
@@ -94,10 +90,9 @@ NamedGraphs.encoded_graph(qf::QuadraticFormNetwork) = encoded_graph(qf.ket)
 
 function DataGraphs.get_vertex_data(
         qf::QuadraticFormNetwork{T, V, I, O}, vertex
-    ) where {T, V, I, TO, O <: ITensorNetworkOperator{TO}}
-    return QuadraticFormGramian{T, TO, I}(
-        qf.ket[vertex], state(qf.operator)[vertex], qf.braname,
-        outputnames(qf.operator), inputnames(qf.operator)
+    ) where {T, V, I, O}
+    return QuadraticFormGramian{T, eltype(O), I}(
+        qf.ket[vertex], qf.operator[vertex], qf.braname
     )
 end
 
@@ -119,15 +114,15 @@ end
 ketnetwork(qf::QuadraticFormNetwork) = qf.ket
 operatornetwork(qf::QuadraticFormNetwork) = qf.operator
 
-# Each output name is renamed to the bra name of the input name it is paired with, using the
-# whole operator's pairing: the per-vertex wrapper drops a pair whose input sits on another vertex.
+# Each output name is renamed to the bra name of the input it is paired with, so the output legs
+# meet the bra layer and the input legs meet the ket layer.
 function operatortensor(g::QuadraticFormGramian)
-    tensor_names = names(g.operator)
+    op = g.operator
     replacements = [
         output => braname(g, input) for
-            (output, input) in zip(g.outputnames, g.inputnames) if output in tensor_names
+            (output, input) in zip(outputnames(op), inputnames(op))
     ]
-    return rename(g.operator, replacements...)
+    return rename(state(op), replacements...)
 end
 
 """
diff --git a/src/itensornetworkoperator.jl b/src/itensornetworkoperator.jl
index c64d1208..48fe077a 100644
--- a/src/itensornetworkoperator.jl
+++ b/src/itensornetworkoperator.jl
@@ -13,9 +13,8 @@ The network equivalent of `ITensorBase.ITensorOperator`: a tensor network of typ
 where `outputnames[i]` is paired with `inputnames[i]`. Applying the operator contracts over
 the input names and leaves the output names.
 
-A pair may straddle two vertices, as it does for a swap or a translation. Indexing returns
-the vertex tensor wrapped as an `ITensorOperator` carrying only the pairs whose two halves
-both sit on that vertex; the remaining legs are dangling on the wrapper.
+The output and input of each pair must sit on the same vertex. Indexing returns the vertex
+tensor wrapped as an `ITensorOperator` carrying the pairs at that vertex.
 """
 struct ITensorNetworkOperator{T, V, I, P <: AbstractITensorNetwork{T, V}} <:
     AbstractITensorNetwork{T, V}
@@ -48,6 +47,16 @@ struct ITensorNetworkOperator{T, V, I, P <: AbstractITensorNetwork{T, V}} <:
                 )
             end
         end
+        for (output, input) in zip(outputnames, inputnames)
+            if only(dimnamevertices(parent, output)) != only(dimnamevertices(parent, input))
+                throw(
+                    ArgumentError(
+                        "operator output $output and its paired input $input must sit on " *
+                            "the same vertex."
+                    )
+                )
+            end
+        end
         return new{T, V, I, typeof(parent)}(parent, outputnames, inputnames)
     end
 end
@@ -77,15 +86,14 @@ NamedGraphs.encoded_graph(op::ITensorNetworkOperator) = encoded_graph(state(op))
 
 # ==================================== DataGraphs.jl ===================================== #
 
-# A pair whose two halves sit on different vertices has no bijection to give this vertex, so
-# only the pairs local to `vertex` become the wrapper's pairing and the rest stay dangling.
+# Both names of a pair sit on one vertex, so the pairs whose output is on `vertex` are its pairing.
 function DataGraphs.get_vertex_data(op::ITensorNetworkOperator, vertex)
     tensor = state(op)[vertex]
     tensor_names = names(tensor)
     outputs = similar(outputnames(op), 0)
     inputs = similar(inputnames(op), 0)
     for (output, input) in zip(outputnames(op), inputnames(op))
-        if output in tensor_names && input in tensor_names
+        if output in tensor_names
             push!(outputs, output)
             push!(inputs, input)
         end
diff --git a/test/test_itensornetworkoperator.jl b/test/test_itensornetworkoperator.jl
index 68b0b6c4..0fb4363c 100644
--- a/test/test_itensornetworkoperator.jl
+++ b/test/test_itensornetworkoperator.jl
@@ -1,7 +1,7 @@
 using DataGraphs: is_vertex_assigned
 using Dictionaries: isinsertable, issettable
 using Graphs: edges, vertices
-using ITensorBase: ITensor, Index, inputnames, name, names, operator, outputnames, state
+using ITensorBase: ITensor, Index, inputnames, name, operator, outputnames, state
 using ITensorNetworksNext:
     ITensorNetwork, ITensorNetworkOperator, operator_support, tensornetwork
 using NamedGraphs: incident_edges, named_path_graph
@@ -58,23 +58,12 @@ end
         @test eltype(op) === typeof(op[first(vertices(g))])
     end
 
-    @testset "a pair may straddle two vertices" begin
-        # A swap: the output at vertex 1 is paired with the input at vertex 2, and vice versa.
+    @testset "a pair must sit on one vertex" begin
+        # A swap pairs the output at vertex 1 with the input at vertex 2, and vice versa.
         a, b = Index(2), Index(2)
         a′, b′ = Index(2), Index(2)
         tn = ITensorNetwork(Dict(1 => randn((a′, a)), 2 => randn((b′, b))))
-        op = operator(tn, [name(a′), name(b′)], [name(b), name(a)])
-
-        @test outputnames(op) == [name(a′), name(b′)]
-        @test inputnames(op) == [name(b), name(a)]
-
-        # Neither pair is local to a vertex, so each vertex wrapper has an empty pairing and
-        # both of its legs stay dangling.
-        for v in (1, 2)
-            @test isempty(outputnames(op[v]))
-            @test isempty(inputnames(op[v]))
-            @test issetequal(names(op[v]), names(state(op)[v]))
-        end
+        @test_throws ArgumentError operator(tn, [name(a′), name(b′)], [name(b), name(a)])
     end
 
     @testset "constructor validation" begin
diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl
index 726de774..c306e664 100644
--- a/test/test_quadraticformnetwork.jl
+++ b/test/test_quadraticformnetwork.jl
@@ -78,7 +78,10 @@ identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s;
         @test gram isa QuadraticFormGramian
         @test eltype(qf) === typeof(gram)
         @test kettensor(gram) === tn[2]
-        @test gram.operator === state(op)[2]
+        # The Gramian holds the operator at its vertex, with that vertex's pairing.
+        @test state(gram.operator) === state(op)[2]
+        @test outputnames(gram.operator) == outputnames(op[2])
+        @test inputnames(gram.operator) == [name(s[2])]
         @test keys(ITensorNetworksNext.layertensors(gram)) == (:ket, :operator, :bra)
         @test contract_network([gram]) ≈
             kettensor(gram) * operatortensor(gram) * bratensor(gram)
@@ -90,30 +93,6 @@ identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s;
         @test axes(gram) == Tuple(inds(gram))
     end
 
-    @testset "crossing operator pair" begin
-        g = named_path_graph(2)
-        tn, l, s = random_state(Float64, g)
-        # Each output is paired with the input on the other vertex.
-        out1, out2 = Index(2), Index(2)
-        optn = ITensorNetwork(Dict(1 => randn((out1, s[1])), 2 => randn((out2, s[2]))))
-        op = operator(optn, [name(out1), name(out2)], [name(s[2]), name(s[1])])
-        qf = QuadraticFormNetwork(tn, op)
-
-        @test braname(qf, name(s[2])) in names(operatortensor(qf[1]))
-        @test braname(qf, name(s[1])) in names(operatortensor(qf[2]))
-        # Both bra site names dangle at vertex 1: one from the operator, one from the bra.
-        @test braname(qf, name(s[2])) in names(qf[1])
-        @test braname(qf, name(s[1])) in names(qf[1])
-        gram = qf[1]
-        @test issetequal(
-            inds(gram), inds(kettensor(gram) * operatortensor(gram) * bratensor(gram))
-        )
-
-        psi = prod(tn)
-        bra = rename(conj(psi), name(s[2]) => name(out1), name(s[1]) => name(out2))
-        @test contract_network(qf)[] ≈ (bra * (optn[1] * optn[2] * psi))[]
-    end
-
     @testset "operator input outside the ket network" begin
         g = named_path_graph(2)
         tn, l, s = random_state(Float64, g)

From 2903a180414f29cf598809e089d4391a72da1524 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 28 Sep 2026 15:18:24 -0400
Subject: [PATCH 17/18] Define shared bilinear-form methods on the abstract
 type; remove `indmap` and the network-level accessors

Graph forwarding, the ket part of `is_vertex_assigned` and `braname` are defined once on
`AbstractBilinearFormNetwork` through `ketnetwork` and a new `branamemap` accessor, and both
constructors build the name map with `select_branames`. `kettensor(bn, v)`, `bratensor(bn, v)`,
`conj_bratensor(bn, v)` and `operatortensor(bn, v)` are removed in favour of calling the
accessor on the Gramian `bn[v]`, and the unused `indmap` is removed.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .../abstractbilinearformnetwork.jl            | 78 ++++++++++++++-----
 src/bilinearforms/normnetwork.jl              | 35 +--------
 src/bilinearforms/quadraticformnetwork.jl     | 37 +--------
 test/test_beliefpropagation.jl                |  2 +-
 test/test_normnetwork.jl                      | 28 +++----
 test/test_quadraticformnetwork.jl             | 31 ++++----
 6 files changed, 91 insertions(+), 120 deletions(-)

diff --git a/src/bilinearforms/abstractbilinearformnetwork.jl b/src/bilinearforms/abstractbilinearformnetwork.jl
index 395b809f..f59823a9 100644
--- a/src/bilinearforms/abstractbilinearformnetwork.jl
+++ b/src/bilinearforms/abstractbilinearformnetwork.jl
@@ -1,7 +1,7 @@
 using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
-using Dictionaries: Dictionaries, isinsertable, issettable
+using Dictionaries: Dictionaries, Dictionary, isinsertable, issettable
 using Graphs: Graphs, edges, vertices
-using ITensorBase: ITensorBase, conj, inds, name, rename, setname
+using ITensorBase: ITensorBase, conj, inds, name, rename
 using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
 
 """
@@ -11,9 +11,9 @@ 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`. The per-vertex accessors [`kettensor`](@ref),
-[`bratensor`](@ref) and, where the subtype has an operator layer, [`operatortensor`](@ref) read
-from that Gramian. The layers as whole networks are returned by [`ketnetwork`](@ref),
+[`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
@@ -28,6 +28,27 @@ are the indices of its layers that no other layer shares, those the layer produc
 """
 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
@@ -43,37 +64,58 @@ The bra-layer index name corresponding to the ket-layer index name `name`. The `
 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(bn::AbstractBilinearFormNetwork, vertex)
     kettensor(g::AbstractGramian)
 
-The ket-layer tensor at `vertex`.
+The ket-layer tensor of the Gramian `g`.
 """
 function kettensor end
-kettensor(bn::AbstractBilinearFormNetwork, vertex) = kettensor(bn[vertex])
 
 """
-    operatortensor(bn::AbstractBilinearFormNetwork, vertex)
     operatortensor(g::AbstractGramian)
 
-The operator-layer tensor at `vertex`, with its index names renamed so that its input legs
-meet the ket layer and its output legs meet the bra layer.
+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
-operatortensor(bn::AbstractBilinearFormNetwork, vertex) = operatortensor(bn[vertex])
 
 conj_bratensor(g::AbstractGramian) = rename(n -> braname(g, n), kettensor(g))
-conj_bratensor(bn::AbstractBilinearFormNetwork, vertex) = conj_bratensor(bn[vertex])
 
 """
-    bratensor(bn::AbstractBilinearFormNetwork, vertex)
     bratensor(g::AbstractGramian)
 
-The bra-layer tensor at `vertex`.
+The bra-layer tensor of the Gramian `g`.
 """
 bratensor(g::AbstractGramian) = conj(conj_bratensor(g))
-bratensor(bn::AbstractBilinearFormNetwork, vertex) = bratensor(bn[vertex])
 
 # Read from `conj_bratensor`, which only renames, so the tensor data is not conjugated.
 brainds(g::AbstractGramian) = conj.(inds(conj_bratensor(g)))
@@ -108,8 +150,6 @@ The bra-layer network of `bn`. Unless a subtype stores its bra layer as a networ
 """
 branetwork(bn::AbstractBilinearFormNetwork) = BraView(bn)
 
-indmap(bn::AbstractBilinearFormNetwork, ind) = setname(conj(ind), braname(bn, name(ind)))
-
 """
     struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <: AbstractITensorNetwork{T, V}
 
@@ -128,7 +168,7 @@ Base.parent(nnv::BraView) = nnv.parent
 
 # ==================================== DataGraphs.jl ===================================== #
 
-DataGraphs.get_vertex_data(nnv::BraView, vertex) = bratensor(parent(nnv), vertex)
+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
diff --git a/src/bilinearforms/normnetwork.jl b/src/bilinearforms/normnetwork.jl
index 868931c0..39da8153 100644
--- a/src/bilinearforms/normnetwork.jl
+++ b/src/bilinearforms/normnetwork.jl
@@ -16,13 +16,7 @@ struct NormNetwork{T, V, I} <: AbstractBilinearFormNetwork{T, V, I}
             ket::ITensorNetwork{T, V, I},
             map::Dictionary{I, I}
         ) where {T, V, I}
-        braname = Dictionary{I, I}()
-        for (name, vertices) in pairs(ket.dimname_vertices)
-            if length(vertices) == 2
-                insert!(braname, name, map[name])
-            end
-        end
-        return new{T, V, I}(ket, braname)
+        return new{T, V, I}(ket, select_branames(ket, map, ()))
     end
 end
 
@@ -38,7 +32,7 @@ struct NormGramian{T, I} <: AbstractGramian
 end
 
 kettensor(g::NormGramian) = g.ket
-braname(g::NormGramian, name) = get(g.braname, name, name)
+branamemap(g::NormGramian) = g.braname
 layertensors(g::NormGramian) = (; ket = kettensor(g), bra = bratensor(g))
 layerinds(g::NormGramian) = (inds(kettensor(g)), brainds(g))
 
@@ -48,39 +42,16 @@ function NormNetwork(tn::ITensorNetwork)
     return NormNetwork(tn, map(uniquename, keys(tn.dimname_vertices)))
 end
 
-# ====================================== Graphs.jl ======================================= #
-
-Graphs.edges(nn::NormNetwork) = edges(nn.ket)
-Graphs.vertices(nn::NormNetwork) = vertices(nn.ket)
-
-# ==================================== NamedGraphs.jl ==================================== #
-
-NamedGraphs.encoded_vertex(nn::NormNetwork, vertex) = encoded_vertex(nn.ket, vertex)
-NamedGraphs.decoded_vertex(nn::NormNetwork, code::Integer) = decoded_vertex(nn.ket, code)
-NamedGraphs.encoded_graph(nn::NormNetwork) = encoded_graph(nn.ket)
-
 # ==================================== DataGraphs.jl ===================================== #
 
-function DataGraphs.is_vertex_assigned(nn::NormNetwork, vertex)
-    return isassigned(nn.ket, vertex)
-end
-
 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 ======================================= #
 
-function braname(nn::NormNetwork, name)
-    if !has_dimname(nn.ket, name)
-        error("index name $name not found underlying tensor network.")
-    end
-    # The indices not stored in `nn.braname` are precisely the site indices, which
-    # get mapped to themselves.
-    return get(nn.braname, name, name)
-end
-
 ketnetwork(nn::NormNetwork) = nn.ket
+branamemap(nn::NormNetwork) = nn.braname
 
 """
     normnetwork(tn::ITensorNetwork, [braname]) -> NormNetwork
diff --git a/src/bilinearforms/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl
index b0ff57f8..ed50ae81 100644
--- a/src/bilinearforms/quadraticformnetwork.jl
+++ b/src/bilinearforms/quadraticformnetwork.jl
@@ -30,13 +30,7 @@ struct QuadraticFormNetwork{T, V, I, O <: ITensorNetworkOperator} <:
         if !issubset(inputnames(operator), keys(ket.dimname_vertices))
             error("every operator input name must be an index name of the ket layer.")
         end
-        acted = Set{I}(inputnames(operator))
-        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
+        braname = select_branames(ket, map, Set{I}(inputnames(operator)))
         return new{T, V, I, typeof(operator)}(ket, operator, braname)
     end
 end
@@ -55,7 +49,7 @@ struct QuadraticFormGramian{T, O, I} <: AbstractGramian
 end
 
 kettensor(g::QuadraticFormGramian) = g.ket
-braname(g::QuadraticFormGramian, name) = get(g.braname, name, name)
+branamemap(g::QuadraticFormGramian) = g.braname
 function layertensors(g::QuadraticFormGramian)
     return (; ket = kettensor(g), operator = operatortensor(g), bra = bratensor(g))
 end
@@ -71,21 +65,6 @@ function QuadraticFormNetwork(ket::ITensorNetwork, operator::ITensorNetworkOpera
     return QuadraticFormNetwork(ket, operator, map(uniquename, keys(ket.dimname_vertices)))
 end
 
-# ====================================== Graphs.jl ======================================= #
-
-Graphs.edges(qf::QuadraticFormNetwork) = edges(qf.ket)
-Graphs.vertices(qf::QuadraticFormNetwork) = vertices(qf.ket)
-
-# ==================================== NamedGraphs.jl ==================================== #
-
-function NamedGraphs.encoded_vertex(qf::QuadraticFormNetwork, vertex)
-    return encoded_vertex(qf.ket, vertex)
-end
-function NamedGraphs.decoded_vertex(qf::QuadraticFormNetwork, code::Integer)
-    return decoded_vertex(qf.ket, code)
-end
-NamedGraphs.encoded_graph(qf::QuadraticFormNetwork) = encoded_graph(qf.ket)
-
 # ==================================== DataGraphs.jl ===================================== #
 
 function DataGraphs.get_vertex_data(
@@ -97,21 +76,13 @@ function DataGraphs.get_vertex_data(
 end
 
 function DataGraphs.is_vertex_assigned(qf::QuadraticFormNetwork, vertex)
-    return isassigned(qf.ket, vertex) && isassigned(qf.operator, vertex)
+    return isassigned(ketnetwork(qf), vertex) && isassigned(operatornetwork(qf), vertex)
 end
 
 # ====================================== interface ======================================= #
 
-function braname(qf::QuadraticFormNetwork, name)
-    if !has_dimname(qf.ket, name)
-        error("index name $name not found underlying tensor network.")
-    end
-    # The indices not stored in `qf.braname` are the dangling ket indices the operator does
-    # not act on, which get mapped to themselves.
-    return get(qf.braname, name, name)
-end
-
 ketnetwork(qf::QuadraticFormNetwork) = qf.ket
+branamemap(qf::QuadraticFormNetwork) = qf.braname
 operatornetwork(qf::QuadraticFormNetwork) = qf.operator
 
 # Each output name is renamed to the bra name of the input it is paired with, so the output legs
diff --git a/test/test_beliefpropagation.jl b/test/test_beliefpropagation.jl
index 0d034e9d..d86bb00e 100644
--- a/test/test_beliefpropagation.jl
+++ b/test/test_beliefpropagation.jl
@@ -343,7 +343,7 @@ end
             ones = message_environment(one, nn)
             for (edge, rest) in ((1 => 2, 2:4), (4 => 3, 1:3))
                 layers =
-                    [[kettensor(nn, v) for v in rest]; [bratensor(nn, v) for v in rest]]
+                    [[kettensor(nn[v]) for v in rest]; [bratensor(nn[v]) for v in rest]]
                 z_rest = contract_network([state(ones[edge]); layers])[]
                 @test z_rest ≈ norm(prod([network[v] for v in rest]))^2 rtol =
                     eps(real(T))^(1 / 3)
diff --git a/test/test_normnetwork.jl b/test/test_normnetwork.jl
index 65c110e0..ec1d64f8 100644
--- a/test/test_normnetwork.jl
+++ b/test/test_normnetwork.jl
@@ -1,11 +1,10 @@
 using DataGraphs: is_vertex_assigned
 using Dictionaries: isinsertable, issettable
 using Graphs: edges, vertices
-using ITensorBase:
-    ITensor, Index, IndexName, LazyITensor, conj, inds, name, names, setname, uniquename
+using ITensorBase: ITensor, Index, IndexName, LazyITensor, inds, name, names, uniquename
 using ITensorNetworksNext: ITensorNetworksNext, BraView, Exact, ITensorNetwork, NormGramian,
     NormNetwork, braname, branetwork, bratensor, conj_bratensor, contract_network,
-    contraction_order, dimnamevertices, indmap, ketnetwork, kettensor, linkaxes, linkinds,
+    contraction_order, dimnamevertices, ketnetwork, kettensor, linkaxes, linkinds,
     linknames, normnetwork, siteaxes, siteinds, sitenames, tensornetwork
 using LinearAlgebra: norm
 using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph
@@ -55,30 +54,25 @@ end
         nn = NormNetwork(tn)
 
         # `kettensor` returns the underlying tensor untouched.
-        @test kettensor(nn, 2) === tn[2]
+        @test kettensor(nn[2]) === tn[2]
 
         # Site indices appear in a single tensor, so they are *not* renamed: the ket and
         # bra layers share them (they get contracted, forming the physical overlap).
         sname = name(s[2])
         @test braname(nn, sname) == sname
-        @test sname in name.(inds(kettensor(nn, 2)))
-        @test sname in name.(inds(conj_bratensor(nn, 2)))
+        @test sname in name.(inds(kettensor(nn[2])))
+        @test sname in name.(inds(conj_bratensor(nn[2])))
 
         # Link indices are shared by two tensors, so they *are* renamed in the bra layer
         # to keep the two layers' bonds distinct.
         lname = name(l[NamedEdge(1 => 2)])
         @test braname(nn, lname) != lname
-        @test lname in name.(inds(kettensor(nn, 2)))
-        @test !(lname in name.(inds(conj_bratensor(nn, 2))))
-        @test braname(nn, lname) in name.(inds(conj_bratensor(nn, 2)))
+        @test lname in name.(inds(kettensor(nn[2])))
+        @test !(lname in name.(inds(conj_bratensor(nn[2]))))
+        @test braname(nn, lname) in name.(inds(conj_bratensor(nn[2])))
 
         # `bra` is the elementwise conjugate of `conj_bratensor` and carries the same indices.
-        @test inds(bratensor(nn, 2)) == inds(conj_bratensor(nn, 2))
-
-        # `indmap` conjugates an index and renames it according to the name map.
-        ind = only(i for i in inds(kettensor(nn, 2)) if name(i) == lname)
-        @test name(indmap(nn, ind)) == braname(nn, name(ind))
-        @test indmap(nn, ind) == setname(conj(ind), braname(nn, name(ind)))
+        @test inds(bratensor(nn[2])) == inds(conj_bratensor(nn[2]))
 
         # Querying the name map with an index name absent from the network errors.
         @test_throws ErrorException braname(nn, name(Index(2)))
@@ -94,7 +88,7 @@ end
 
         lname = name(l[NamedEdge(1 => 2)])
         @test braname(nn, lname) == custom[lname]
-        @test braname(nn, lname) in name.(inds(conj_bratensor(nn, 2)))
+        @test braname(nn, lname) in name.(inds(conj_bratensor(nn[2])))
     end
 
     @testset "`ketnetwork` / `branetwork`" begin
@@ -111,7 +105,7 @@ end
         @test issetequal(vertices(bv), vertices(tn))
         @test issetequal(edges(bv), edges(tn))
         for v in vertices(tn)
-            @test inds(bv[v]) == inds(bratensor(nn, v))
+            @test inds(bv[v]) == inds(bratensor(nn[v]))
         end
         @test is_vertex_assigned(bv, 1)
 
diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl
index c306e664..4dd3122e 100644
--- a/test/test_quadraticformnetwork.jl
+++ b/test/test_quadraticformnetwork.jl
@@ -2,10 +2,10 @@ using DataGraphs: is_vertex_assigned
 using Dictionaries: isinsertable, issettable
 using Graphs: edges, vertices
 using ITensorBase: ITensor, Index, IndexName, conj, inds, inputnames, name, names, operator,
-    outputnames, rename, setname, state, uniquename
+    outputnames, rename, state, uniquename
 using ITensorNetworksNext: ITensorNetworksNext, BraView, ITensorNetwork, NormNetwork,
     QuadraticFormGramian, QuadraticFormNetwork, braname, branetwork, bratensor,
-    conj_bratensor, contract_network, indmap, ketnetwork, kettensor, operatornetwork,
+    conj_bratensor, contract_network, ketnetwork, kettensor, operatornetwork,
     operatortensor, quadraticformnetwork, tensornetwork
 using LinearAlgebra: I, norm
 using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph
@@ -110,39 +110,34 @@ identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s;
         qf = QuadraticFormNetwork(tn, op)
 
         # `kettensor` returns the underlying tensor untouched.
-        @test kettensor(qf, 2) === tn[2]
+        @test kettensor(qf[2]) === tn[2]
 
         # Unlike the norm network, the site indices *are* renamed in the bra layer: the
         # operator sits between the two layers, so they no longer contract directly.
         sname = name(s[2])
         @test braname(qf, sname) != sname
-        @test sname in name.(inds(kettensor(qf, 2)))
-        @test !(sname in name.(inds(conj_bratensor(qf, 2))))
-        @test braname(qf, sname) in name.(inds(conj_bratensor(qf, 2)))
+        @test sname in name.(inds(kettensor(qf[2])))
+        @test !(sname in name.(inds(conj_bratensor(qf[2]))))
+        @test braname(qf, sname) in name.(inds(conj_bratensor(qf[2])))
 
         # Link indices are shared by two tensors, so they are renamed in the bra layer to
         # keep the two layers' bonds distinct.
         lname = name(l[NamedEdge(1 => 2)])
         @test braname(qf, lname) != lname
-        @test lname in name.(inds(kettensor(qf, 2)))
-        @test !(lname in name.(inds(conj_bratensor(qf, 2))))
-        @test braname(qf, lname) in name.(inds(conj_bratensor(qf, 2)))
+        @test lname in name.(inds(kettensor(qf[2])))
+        @test !(lname in name.(inds(conj_bratensor(qf[2]))))
+        @test braname(qf, lname) in name.(inds(conj_bratensor(qf[2])))
 
         # The operator's input name meets the ket and its output name is renamed to meet
         # the bra.
-        o = operatortensor(qf, 2)
+        o = operatortensor(qf[2])
         @test sname in names(o)
         @test braname(qf, sname) in names(o)
         @test inputnames(op) == [name(s[v]) for v in vertices(g)]
 
         # `bratensor` is the elementwise conjugate of `conj_bratensor` and carries the same
         # indices.
-        @test inds(bratensor(qf, 2)) == inds(conj_bratensor(qf, 2))
-
-        # `indmap` conjugates an index and renames it according to the name map.
-        ind = only(i for i in inds(kettensor(qf, 2)) if name(i) == lname)
-        @test name(indmap(qf, ind)) == braname(qf, name(ind))
-        @test indmap(qf, ind) == setname(conj(ind), braname(qf, name(ind)))
+        @test inds(bratensor(qf[2])) == inds(conj_bratensor(qf[2]))
 
         # Querying the name map with an index name absent from the ket layer errors.
         @test_throws ErrorException braname(qf, name(Index(2)))
@@ -159,7 +154,7 @@ identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s;
 
         lname = name(l[NamedEdge(1 => 2)])
         @test braname(qf, lname) == custom[lname]
-        @test braname(qf, lname) in name.(inds(conj_bratensor(qf, 2)))
+        @test braname(qf, lname) in name.(inds(conj_bratensor(qf[2])))
 
         sname = name(s[2])
         @test braname(qf, sname) == custom[sname]
@@ -184,7 +179,7 @@ identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s;
         @test !isinsertable(bv)
         @test is_vertex_assigned(bv, 1)
         for v in vertices(tn)
-            @test inds(bv[v]) == inds(bratensor(qf, v))
+            @test inds(bv[v]) == inds(bratensor(qf[v]))
         end
     end
 

From f9bdc22ebf1a645e2b4f915c11d80440a71e4908 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Tue, 29 Sep 2026 10:54:24 -0400
Subject: [PATCH 18/18] Store the `ITensorNetworkOperator` pairing as a
 `Dictionary`

The output and input name vectors are replaced by one map from each output name to its input
name, so `getindex` looks up a vertex's pairs by name instead of scanning every pair. A repeated
output name now throws an `ArgumentError`.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/itensornetworkoperator.jl       | 29 ++++++++++++-----------------
 test/test_itensornetworkoperator.jl |  4 ++++
 2 files changed, 16 insertions(+), 17 deletions(-)

diff --git a/src/itensornetworkoperator.jl b/src/itensornetworkoperator.jl
index 48fe077a..c8b68686 100644
--- a/src/itensornetworkoperator.jl
+++ b/src/itensornetworkoperator.jl
@@ -1,5 +1,5 @@
 using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned
-using Dictionaries: Dictionaries
+using Dictionaries: Dictionaries, Dictionary
 using Graphs: Graphs, AbstractGraph, edges, vertices
 using ITensorBase: ITensorBase, NamedTensorOperator, inputnames, names, nametype, operator,
     outputnames, state
@@ -10,7 +10,7 @@ using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex
 
 The network equivalent of `ITensorBase.ITensorOperator`: a tensor network of type
 `P <: AbstractITensorNetwork{T, V}` together with a pairing of its dangling index names,
-where `outputnames[i]` is paired with `inputnames[i]`. Applying the operator contracts over
+stored as a map from each output name to its input name. Applying the operator contracts over
 the input names and leaves the output names.
 
 The output and input of each pair must sit on the same vertex. Indexing returns the vertex
@@ -19,8 +19,7 @@ tensor wrapped as an `ITensorOperator` carrying the pairs at that vertex.
 struct ITensorNetworkOperator{T, V, I, P <: AbstractITensorNetwork{T, V}} <:
     AbstractITensorNetwork{T, V}
     parent::P
-    outputnames::Vector{I}
-    inputnames::Vector{I}
+    pairing::Dictionary{I, I}
     function ITensorNetworkOperator(
             parent::AbstractITensorNetwork{T, V}, outputnames, inputnames
         ) where {T, V}
@@ -36,6 +35,9 @@ struct ITensorNetworkOperator{T, V, I, P <: AbstractITensorNetwork{T, V}} <:
                 )
             )
         end
+        if !allunique(outputnames)
+            throw(ArgumentError("each operator output name must appear only once."))
+        end
         for opname in Iterators.flatten((outputnames, inputnames))
             nvertices = length(dimnamevertices(parent, opname))
             if nvertices != 1
@@ -57,7 +59,8 @@ struct ITensorNetworkOperator{T, V, I, P <: AbstractITensorNetwork{T, V}} <:
                 )
             end
         end
-        return new{T, V, I, typeof(parent)}(parent, outputnames, inputnames)
+        pairing = Dictionary(outputnames, inputnames)
+        return new{T, V, I, typeof(parent)}(parent, pairing)
     end
 end
 
@@ -89,16 +92,8 @@ NamedGraphs.encoded_graph(op::ITensorNetworkOperator) = encoded_graph(state(op))
 # Both names of a pair sit on one vertex, so the pairs whose output is on `vertex` are its pairing.
 function DataGraphs.get_vertex_data(op::ITensorNetworkOperator, vertex)
     tensor = state(op)[vertex]
-    tensor_names = names(tensor)
-    outputs = similar(outputnames(op), 0)
-    inputs = similar(inputnames(op), 0)
-    for (output, input) in zip(outputnames(op), inputnames(op))
-        if output in tensor_names
-            push!(outputs, output)
-            push!(inputs, input)
-        end
-    end
-    return operator(tensor, outputs, inputs)
+    outputs = filter(name -> haskey(op.pairing, name), names(tensor))
+    return operator(tensor, outputs, [op.pairing[output] for output in outputs])
 end
 
 function DataGraphs.is_vertex_assigned(op::ITensorNetworkOperator, vertex)
@@ -114,8 +109,8 @@ Dictionaries.isinsertable(::ITensorNetworkOperator) = false
 
 ITensorBase.state(op::ITensorNetworkOperator) = op.parent
 Base.parent(op::ITensorNetworkOperator) = state(op)
-ITensorBase.outputnames(op::ITensorNetworkOperator) = op.outputnames
-ITensorBase.inputnames(op::ITensorNetworkOperator) = op.inputnames
+ITensorBase.outputnames(op::ITensorNetworkOperator) = collect(keys(op.pairing))
+ITensorBase.inputnames(op::ITensorNetworkOperator) = collect(op.pairing)
 
 function operator_support(tn::AbstractGraph, op::ITensorNetworkOperator)
     return operator_support_names(tn, inputnames(op))
diff --git a/test/test_itensornetworkoperator.jl b/test/test_itensornetworkoperator.jl
index 0fb4363c..ae6eb51c 100644
--- a/test/test_itensornetworkoperator.jl
+++ b/test/test_itensornetworkoperator.jl
@@ -77,6 +77,10 @@ end
         @test_throws ArgumentError operator(
             tn, [name(out[v]) for v in vs], [name(inp[vs[1]])]
         )
+        # An output name can be paired with only one input.
+        @test_throws ArgumentError operator(
+            tn, [name(out[vs[1]]), name(out[vs[1]])], [name(inp[vs[1]]), name(inp[vs[1]])]
+        )
 
         # An operator leg must be dangling: a link name touches two vertices.
         l = Index(2)