From 35e491eb9b72392cfcf5273b749b182d4813d574 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Thu, 9 Jul 2026 13:20:54 -0400
Subject: [PATCH 01/16] `add_edge!` and `rem_edge!` on the `ITensorNetwork`
 type now return false instead of erroring

This is inline with the `Graphs` behaviour.
---
 src/tensornetwork.jl       | 14 +++-----------
 test/test_tensornetwork.jl | 10 ++++++++--
 2 files changed, 11 insertions(+), 13 deletions(-)

diff --git a/src/tensornetwork.jl b/src/tensornetwork.jl
index ba96462..3971da1 100644
--- a/src/tensornetwork.jl
+++ b/src/tensornetwork.jl
@@ -171,17 +171,9 @@ function DataGraphs.underlying_graph_type(type::Type{<:ITensorNetwork{T, V}}) wh
     return fieldtype(type, :underlying_graph)
 end
 
-function Graphs.rem_edge!(::ITensorNetwork, _edge)
-    return throw(
-        ErrorException("removing edges from the `ITensorNetwork` type is not supported.")
-    )
-end
-
-function Graphs.add_edge!(::ITensorNetwork, _edge)
-    return throw(
-        ErrorException("Adding edges to the `ITensorNetwork` type is not supported.")
-    )
-end
+# Can't add/remove edges from `ITensorNetwork` as graph topology fixed by indices.
+Graphs.rem_edge!(::ITensorNetwork, _edge) = false
+Graphs.add_edge!(::ITensorNetwork, _edge) = false
 
 # PERF: fast lookup compared to `AbstractITensorNetwork` fallback.
 dimnamevertices(tn::ITensorNetwork, name) = tn.dimname_vertices[name]
diff --git a/test/test_tensornetwork.jl b/test/test_tensornetwork.jl
index ab333b1..bc060d6 100644
--- a/test/test_tensornetwork.jl
+++ b/test/test_tensornetwork.jl
@@ -48,8 +48,14 @@ using Test: @test, @test_throws, @testset
         @test_throws MethodError tn[e] = randn(2, 2)
         @test_throws MethodError tn[src(e) => dst(e)] = randn(2, 2)
 
-        # `rem_edge!` is intentionally unimplemented.
-        @test_throws ErrorException rem_edge!(tn, (1, 1) => (2, 1))
+        # `rem_edge!` and `add_edge!` are intentionally unimplemented; they return
+        # `false` without modifying the network.
+        @test rem_edge!(tn, (1, 1) => (2, 1)) == false
+        @test has_edge(tn, (1, 1) => (2, 1))
+        @test ne(tn) == 1
+        @test add_edge!(tn, (2, 1) => (2, 2)) == false
+        @test !has_edge(tn, (2, 1) => (2, 2))
+        @test ne(tn) == 1
 
         tn[1, 1] = randn(Index(2))
         tn[2, 1] = randn(Index(2))

From 25b18facc837be3839caf3cc487168f412d8a8ee Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Thu, 9 Jul 2026 13:36:16 -0400
Subject: [PATCH 02/16] `dimnamevertices` for `ITensorNetwork` now returns
 empty set if index not in dictionary

This is now consistant with the fallback defn of `dimnamevertices`.
would error previously.
---
 src/tensornetwork.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/tensornetwork.jl b/src/tensornetwork.jl
index 3971da1..994b650 100644
--- a/src/tensornetwork.jl
+++ b/src/tensornetwork.jl
@@ -176,7 +176,9 @@ Graphs.rem_edge!(::ITensorNetwork, _edge) = false
 Graphs.add_edge!(::ITensorNetwork, _edge) = false
 
 # PERF: fast lookup compared to `AbstractITensorNetwork` fallback.
-dimnamevertices(tn::ITensorNetwork, name) = tn.dimname_vertices[name]
+function dimnamevertices(tn::ITensorNetwork, name)
+    return get(tn.dimname_vertices, name, Set{vertextype(tn)}())
+end
 
 # PERF: fast lookup compared to `AbstractITensorNetwork` fallback.
 has_dimname(tn::ITensorNetwork, name) = haskey(tn.dimname_vertices, name)

From 5dc2aeb57fc0a7da014a2257ed023b868dd9fc58 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Thu, 9 Jul 2026 13:36:40 -0400
Subject: [PATCH 03/16] New function `supportof` that gives the vertex support
 of an operator on a tensor network.

---
 src/abstracttensornetwork.jl | 23 ++++++++++++++++++++++-
 src/apply/apply_operators.jl | 35 +++++++++++++++++++----------------
 2 files changed, 41 insertions(+), 17 deletions(-)

diff --git a/src/abstracttensornetwork.jl b/src/abstracttensornetwork.jl
index 4ea4c34..20fbcb1 100644
--- a/src/abstracttensornetwork.jl
+++ b/src/abstracttensornetwork.jl
@@ -4,7 +4,8 @@ using DataGraphs: DataGraphs, AbstractDataGraph, AbstractVertexDataGraph, edge_d
 using Dictionaries: Dictionary
 using Graphs: Graphs, AbstractEdge, AbstractGraph, add_edge!, add_vertex!, dst, edges,
     edgetype, ne, neighbors, nv, rem_edge!, src, vertices
-using ITensorBase: dimnames, inds, name, named, nametype, prime, uniquename, unnamedtype
+using ITensorBase: ITensorOperator, dimnames, domainnames, inds, name, named, nametype,
+    prime, uniquename, unnamedtype
 using LinearAlgebra: LinearAlgebra
 using MacroTools: @capture
 using NamedGraphs:
@@ -130,3 +131,23 @@ function insertlink!(tn::AbstractGraph, e)
 
     return tn
 end
+
+function supportof(tn::AbstractGraph, op::ITensorOperator)
+    support = Base.Generator(domainnames(op)) do name
+        vertices = dimnamevertices(tn, name)
+
+        length(vertices) == 1 && return only(vertices)
+
+        if length(vertices) == 0
+            throw(ArgumentError("operator dim name $name not found in tensor network."))
+        elseif length(vertices) > 1
+            throw(
+                ArgumentError(
+                    "operator dim name $name associated with multiple vertices in tensor network."
+                )
+            )
+        end
+    end
+
+    return Set(support)
+end
diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index f824b9d..5bbfbf4 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -2,9 +2,9 @@ using .AlgorithmsInterfaceExtensions: AlgorithmsInterfaceExtensions as AIE
 using AlgorithmsInterface: AlgorithmsInterface as AI
 using Base: @kwdef
 using Graphs: dst, src, vertices
-using ITensorBase: ITensorBase as ITB, AbstractITensor, dimnames, inputnames, operator,
-    outputnames, replacedimnames
-using LinearAlgebra: norm
+using ITensorBase:
+    ITensorBase, AbstractITensor, apply, dimnames, domainnames, operator, replacedimnames
+using LinearAlgebra: norm, normalize!
 using MatrixAlgebraKit: eigh_full, project_hermitian, qr_compact, svd_trunc
 using NamedGraphs: boundary_edges
 using TensorAlgebra.MatrixAlgebra: invsqrth_safe, sqrth_safe
@@ -239,12 +239,15 @@ function apply_gate_bp!(
         dest::AbstractITensorNetwork, op::AbstractITensor,
         state::AbstractITensorNetwork, env; kwargs...
     )
-    op_in = inputnames(op)
-    vs = [v for v in vertices(state) if !isempty(intersect(op_in, sitenames(state, v)))]
-    isempty(vs) && throw(
+    vertices = supportof(state, op)
+
+    isempty(vertices) && throw(
         ArgumentError("operator shares no indices with the tensor network")
     )
-    return apply_gate_bp_nsite!(Val(length(vs)), dest, op, state, env, vs; kwargs...)
+
+    N = Val(length(vertices))
+
+    return apply_gate_bp_nsite!(N, dest, op, state, env, vertices; kwargs...)
 end
 
 function apply_gate_bp_nsite!(
@@ -256,11 +259,11 @@ end
 
 function apply_gate_bp_nsite!(
         ::Val{1}, dest::AbstractITensorNetwork, op::AbstractITensor,
-        state::AbstractITensorNetwork, env, vs;
+        state::AbstractITensorNetwork, env, vertices;
         normalize, kwargs...
     )
-    v = only(vs)
-    ψv = ITB.apply(op, state[v])
+    vertex = only(vertices)
+    ψv = apply(op, state[vertex])
     if normalize
         gauges = [
             message_root(env[e])
@@ -268,13 +271,13 @@ function apply_gate_bp_nsite!(
         ]
         ψv /= norm(prod([[ψv]; gauges]))
     end
-    dest[v] = ψv
+    dest[vertex] = ψv
     return dest
 end
 
 function apply_gate_bp_nsite!(
         ::Val{2}, dest::AbstractITensorNetwork, op::AbstractITensor,
-        state::AbstractITensorNetwork, env, vs;
+        state::AbstractITensorNetwork, env, vertices;
         trunc, normalize
     )
     v1, v2 = vs
@@ -289,11 +292,11 @@ function apply_gate_bp_nsite!(
 
     Q_v1, R_v1 = qr_compact(ψ_v1, setdiff(dimnames(ψ_v1), dimnames(ψ_v2), dimnames(op)))
     Q_v2, R_v2 = qr_compact(ψ_v2, setdiff(dimnames(ψ_v2), dimnames(ψ_v1), dimnames(op)))
-    op_R_v1v2 = ITB.apply(op, R_v1 * R_v2)
+    op_R_v1v2 = apply(op, R_v1 * R_v2)
     U_v1, S, U_v2 = svd_trunc(op_R_v1v2, setdiff(dimnames(R_v1), dimnames(R_v2)); trunc)
-    if normalize
-        S = S / norm(S)
-    end
+
+    normalize && normalize!(S)
+
     name_v1, name_v2 = dimnames(S)
     sqrt_S = sqrth_safe(S, (name_v1,), (name_v2,); atol = 0, rtol = 0)
     R_v1 = replacedimnames(U_v1 * sqrt_S, name_v2 => name_v1)

From 9e57b58d3c747caccbc83b773cd11d8aab0eacd2 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Thu, 9 Jul 2026 15:14:47 -0400
Subject: [PATCH 04/16] Refactor `ITensorNetwork` topology functions; add index
 precheck for setting tensors

---
 src/tensornetwork.jl | 48 +++++++++++++++++++++++++++++++-------------
 1 file changed, 34 insertions(+), 14 deletions(-)

diff --git a/src/tensornetwork.jl b/src/tensornetwork.jl
index 994b650..964b682 100644
--- a/src/tensornetwork.jl
+++ b/src/tensornetwork.jl
@@ -116,38 +116,60 @@ DataGraphs.is_edge_assigned(::ITensorNetwork, _edge) = false
 
 DataGraphs.get_vertex_data(tn::ITensorNetwork, v) = tn.tensors[v]
 
+function check_incoming_dimnames(tn, tensor, vertex)
+    for name in dimnames(tensor)
+        vertices = get(tn.dimname_vertices, name, Set())
+        if length(setdiff(vertices, Set([vertex]))) > 1
+            throw(
+                ArgumentError(
+                    "index $name can appear in at most one existing tensor"
+                )
+            )
+        end
+    end
+    return nothing
+end
+
 function DataGraphs.insert_vertex_data!(tn::ITensorNetwork, vertex, tensor)
+    check_incoming_dimnames(tn, tensor, vertex)
     add_vertex!(tn.underlying_graph, vertex)
-    set!_tensornetwork(tn, vertex, tensor)
+    update_tensornetwork_metadata!(tn, vertex, tensor)
+    insert!(tn.tensors, vertex, tensor)
     return tn
 end
 
 function DataGraphs.set_vertex_data!(tn::ITensorNetwork, tensor, vertex)
-    set!_tensornetwork(tn, vertex, tensor)
+    check_incoming_dimnames(tn, tensor, vertex)
+    update_tensornetwork_metadata!(tn, vertex, tensor)
+    set!(tn.tensors, vertex, tensor)
     return tn
 end
 
-# "upsert"
-function set!_tensornetwork(tn::ITensorNetwork, vertex, tensor)
-    newinds = dimnames(tensor)
+function update_tensornetwork_metadata!(tn, vertex, tensor)
+    oldnames = isassigned(tn, vertex) ? dimnames(tn[vertex]) : Set()
+    newnames = dimnames(tensor)
 
-    oldinds = get(mapview(dimnames, tn.tensors), vertex, Set())
+    update_tensornetwork_metadata!(tn, vertex, oldnames, newnames)
 
+    return tn
+end
+
+function update_tensornetwork_metadata!(tn, vertex, oldinds, newinds)
     # Only have to deal with the indices that aren't shared.
-    for ind in symdiff(oldinds, newinds)
-        if ind in oldinds
-            delete_ind_edge!(tn, ind)
-            delete_ind_vertex!(tn, ind, vertex)
+    for name in symdiff(oldinds, newinds)
+        if name in oldinds
+            delete_ind_edge!(tn, name)
+            delete_ind_vertex!(tn, name, vertex)
             continue
         end
 
         # Now `ind` must be a new index that's not in `oldinds`
 
-        vertex_list = get!(tn.dimname_vertices, ind, Set())
+        vertex_list = get!(tn.dimname_vertices, name, Set())
         if length(vertex_list) > 1
             throw(
                 ArgumentError(
-                    "index $ind can appear in at most one existing tensor, got $(length(vertex_list))."
+                    "index $name can appear in at most one existing tensor, got $(length(vertex_list))."
                 )
             )
         end
@@ -160,8 +182,6 @@ function set!_tensornetwork(tn::ITensorNetwork, vertex, tensor)
         end
     end
 
-    set!(tn.tensors, vertex, tensor)
-
     return tn
 end
 

From fc2b71fe84c0a148f0065c3947207a601840d688 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 20 Jul 2026 15:50:42 -0400
Subject: [PATCH 05/16] Refactor sum of log scalars into own function.

Avoids some minor code duplication.
---
 src/beliefpropagation/messagecache.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index fc6bc61..8100a86 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/src/beliefpropagation/messagecache.jl
@@ -175,23 +175,23 @@ function region_scalar(factors, messages, region)
     return mapreduce(vertex -> vertex_scalar(factors, messages, vertex), *, region)
 end
 
+function sum_log_scalars(terms)
+    if any(t -> real(t) < 0, terms)
+        terms = complex.(terms)
+    end
+    return sum(log.(terms))
+end
+
 # We need a graph structure here, so assume `factors` is a graph.
 function bethe_free_energy(factors, messages)
     numerator_terms = vertex_scalars(factors, messages)
     denominator_terms = edge_scalars(messages)
 
-    if any(t -> real(t) < 0, numerator_terms)
-        numerator_terms = complex.(numerator_terms)
-    end
-    if any(t -> real(t) < 0, denominator_terms)
-        denominator_terms = complex.(denominator_terms)
-    end
-
     if any(iszero, denominator_terms)
         return -Inf
     end
 
-    return sum(log.(numerator_terms)) - sum(log.(denominator_terms))
+    return sum_log_scalars(numerator_terms) - sum_log_scalars(denominator_terms)
 end
 
 # ===================================== NormNetwork ====================================== #

From 0c59060719af01ebe51ace94c759930815d90465 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 20 Jul 2026 15:51:21 -0400
Subject: [PATCH 06/16] The `finalize_substate` function now dispatches on
 subsolve rather than solve.

---
 .../AlgorithmsInterfaceExtensions.jl                          | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl b/src/AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl
index da7bf76..917fef9 100644
--- a/src/AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl
+++ b/src/AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl
@@ -2,6 +2,8 @@ module AlgorithmsInterfaceExtensions
 
 using AlgorithmsInterface: AlgorithmsInterface as AI
 
+abstract type NestedProblem <: AI.Problem end
+
 # ============================ NestedAlgorithm =============================================
 
 abstract type NestedAlgorithm <: AI.Algorithm end
@@ -27,7 +29,7 @@ end
 function AI.step!(problem::AI.Problem, algorithm::NestedAlgorithm, state::AI.State)
     subproblem, subalgorithm, substate = initialize_subsolve(problem, algorithm, state)
     AI.solve!(subproblem, subalgorithm, substate)
-    finalize_substate!(problem, algorithm, state, substate)
+    finalize_substate!(subproblem, subalgorithm, substate, state)
     return state
 end
 

From 01b76eed8e7c9dac489efaf6ac081b03213131b6 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Tue, 21 Jul 2026 17:52:49 -0400
Subject: [PATCH 07/16] Upgrade to ITensorBase v0.13

Fix imports in `apply_operators.jl`
---
 src/abstracttensornetwork.jl | 4 ++--
 src/apply/apply_operators.jl | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/abstracttensornetwork.jl b/src/abstracttensornetwork.jl
index 20fbcb1..e9b5c6b 100644
--- a/src/abstracttensornetwork.jl
+++ b/src/abstracttensornetwork.jl
@@ -4,7 +4,7 @@ using DataGraphs: DataGraphs, AbstractDataGraph, AbstractVertexDataGraph, edge_d
 using Dictionaries: Dictionary
 using Graphs: Graphs, AbstractEdge, AbstractGraph, add_edge!, add_vertex!, dst, edges,
     edgetype, ne, neighbors, nv, rem_edge!, src, vertices
-using ITensorBase: ITensorOperator, dimnames, domainnames, inds, name, named, nametype,
+using ITensorBase: ITensorOperator, dimnames, inputnames, inds, name, named, nametype,
     prime, uniquename, unnamedtype
 using LinearAlgebra: LinearAlgebra
 using MacroTools: @capture
@@ -133,7 +133,7 @@ function insertlink!(tn::AbstractGraph, e)
 end
 
 function supportof(tn::AbstractGraph, op::ITensorOperator)
-    support = Base.Generator(domainnames(op)) do name
+    support = Base.Generator(inputnames(op)) do name
         vertices = dimnamevertices(tn, name)
 
         length(vertices) == 1 && return only(vertices)
diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index 5bbfbf4..372f996 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -2,8 +2,8 @@ using .AlgorithmsInterfaceExtensions: AlgorithmsInterfaceExtensions as AIE
 using AlgorithmsInterface: AlgorithmsInterface as AI
 using Base: @kwdef
 using Graphs: dst, src, vertices
-using ITensorBase:
-    ITensorBase, AbstractITensor, apply, dimnames, domainnames, operator, replacedimnames
+using ITensorBase: ITensorBase as ITB, AbstractITensor, apply, dimnames, inputnames, operator,
+    outputnames, replacedimnames
 using LinearAlgebra: norm, normalize!
 using MatrixAlgebraKit: eigh_full, project_hermitian, qr_compact, svd_trunc
 using NamedGraphs: boundary_edges

From 511703ed09699afac5aee83d5cdf81e79ca10c29 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Tue, 8 Sep 2026 09:04:34 -0400
Subject: [PATCH 08/16] Fix rename of `vs` to `vertices` in function body.

---
 src/apply/apply_operators.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index 372f996..9edf9ba 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -267,7 +267,7 @@ function apply_gate_bp_nsite!(
     if normalize
         gauges = [
             message_root(env[e])
-                for e in boundary_edges(state, vs; dir = :in)
+                for e in boundary_edges(state, vertices; dir = :in)
         ]
         ψv /= norm(prod([[ψv]; gauges]))
     end
@@ -280,8 +280,8 @@ function apply_gate_bp_nsite!(
         state::AbstractITensorNetwork, env, vertices;
         trunc, normalize
     )
-    v1, v2 = vs
-    edges_in = boundary_edges(state, vs; dir = :in)
+    v1, v2 = vertices
+    edges_in = boundary_edges(state, vertices; dir = :in)
     siv_v1 = [message_gauge(env[e]) for e in edges_in if dst(e) == v1]
     siv_v2 = [message_gauge(env[e]) for e in edges_in if dst(e) == v2]
     gauges_v1, inv_gauges_v1 = first.(siv_v1), conj.(last.(siv_v1))

From d9db4865340f5fe67795cacfe6ff47307cd18789 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Tue, 8 Sep 2026 10:35:18 -0400
Subject: [PATCH 09/16] Function `supportof` can now return an empty set; add
 tests.

---
 src/abstracttensornetwork.jl | 16 ++++++++--------
 test/test_tensornetwork.jl   | 32 ++++++++++++++++++++++++++++++--
 2 files changed, 38 insertions(+), 10 deletions(-)

diff --git a/src/abstracttensornetwork.jl b/src/abstracttensornetwork.jl
index e9b5c6b..275068e 100644
--- a/src/abstracttensornetwork.jl
+++ b/src/abstracttensornetwork.jl
@@ -4,7 +4,7 @@ using DataGraphs: DataGraphs, AbstractDataGraph, AbstractVertexDataGraph, edge_d
 using Dictionaries: Dictionary
 using Graphs: Graphs, AbstractEdge, AbstractGraph, add_edge!, add_vertex!, dst, edges,
     edgetype, ne, neighbors, nv, rem_edge!, src, vertices
-using ITensorBase: ITensorOperator, dimnames, inputnames, inds, name, named, nametype,
+using ITensorBase: ITensorOperator, dimnames, inds, inputnames, name, named, nametype,
     prime, uniquename, unnamedtype
 using LinearAlgebra: LinearAlgebra
 using MacroTools: @capture
@@ -133,21 +133,21 @@ function insertlink!(tn::AbstractGraph, e)
 end
 
 function supportof(tn::AbstractGraph, op::ITensorOperator)
-    support = Base.Generator(inputnames(op)) do name
-        vertices = dimnamevertices(tn, name)
+    support = Set{vertextype(tn)}()
 
-        length(vertices) == 1 && return only(vertices)
+    for name in inputnames(op)
+        vertices = dimnamevertices(tn, name)
 
-        if length(vertices) == 0
-            throw(ArgumentError("operator dim name $name not found in tensor network."))
-        elseif length(vertices) > 1
+        if length(vertices) > 1
             throw(
                 ArgumentError(
                     "operator dim name $name associated with multiple vertices in tensor network."
                 )
             )
         end
+
+        union!(support, vertices)
     end
 
-    return Set(support)
+    return support
 end
diff --git a/test/test_tensornetwork.jl b/test/test_tensornetwork.jl
index bc060d6..29c600a 100644
--- a/test/test_tensornetwork.jl
+++ b/test/test_tensornetwork.jl
@@ -2,9 +2,9 @@ using DataGraphs:
     DataGraph, assigned_edge_data, assigned_vertex_data, underlying_graph, vertex_data
 using Graphs: add_edge!, add_vertex!, dst, edges, edgetype, has_edge, has_vertex,
     is_directed, ne, nv, rem_edge!, rem_vertex!, src, vertices
-using ITensorBase: Index, LazyITensor, inds
+using ITensorBase: Index, LazyITensor, inds, operator
 using ITensorNetworksNext: ITensorNetwork, has_ind, linkaxes, linkinds, linknames, siteaxes,
-    siteinds, sitenames, tensornetwork
+    siteinds, sitenames, supportof, tensornetwork
 using NamedGraphs: convert_vertextype, incident_edges, named_grid, named_path_graph,
     similar_graph, subgraph, vertextype
 using Test: @test, @test_throws, @testset
@@ -125,6 +125,34 @@ using Test: @test, @test_throws, @testset
         @test sitenames(tn, 3) == [s[3].name]
     end
 
+    @testset "`supportof`" begin
+        g = named_path_graph(3)
+        l = Dict(e => Index(2) for e in edges(g))
+        l = merge(l, Dict(reverse(e) => l[e] for e in edges(g)))
+        s = Dict(v => Index(2) for v in vertices(g))
+        tn = tensornetwork(vertices(g)) do v
+            is = map(e -> l[e], incident_edges(g, v))
+            return randn((s[v], is...))
+        end
+
+        o1 = operator(randn(2, 2), (Index(2),), (s[2],))
+        @test supportof(tn, o1) == Set([2])
+
+        o12 = operator(randn(2, 2, 2, 2), (Index(2), Index(2)), (s[1], s[2]))
+        @test supportof(tn, o12) == Set([1, 2])
+
+        # An input name that no tensor in the network carries contributes no vertex.
+        o_absent = operator(randn(2, 2), (Index(2),), (Index(2),))
+        @test isempty(supportof(tn, o_absent))
+
+        o_partial = operator(randn(2, 2, 2, 2), (Index(2), Index(2)), (s[3], Index(2)))
+        @test supportof(tn, o_partial) == Set([3])
+
+        # A link index is carried by both endpoints of its edge.
+        o_link = operator(randn(2, 2), (Index(2),), (l[first(edges(g))],))
+        @test_throws ArgumentError supportof(tn, o_link)
+    end
+
     @testset "`subgraph`" begin
         g = named_grid((3,))
 

From 7ca185255959043a331640e00510dcf77eeb39c3 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 14 Sep 2026 13:32:51 -0400
Subject: [PATCH 10/16] Rename function `supportof` to `operator_support`.

---
 src/abstracttensornetwork.jl |  2 +-
 src/apply/apply_operators.jl |  2 +-
 test/test_tensornetwork.jl   | 14 +++++++-------
 3 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/abstracttensornetwork.jl b/src/abstracttensornetwork.jl
index 275068e..3a314df 100644
--- a/src/abstracttensornetwork.jl
+++ b/src/abstracttensornetwork.jl
@@ -132,7 +132,7 @@ function insertlink!(tn::AbstractGraph, e)
     return tn
 end
 
-function supportof(tn::AbstractGraph, op::ITensorOperator)
+function operator_support(tn::AbstractGraph, op::ITensorOperator)
     support = Set{vertextype(tn)}()
 
     for name in inputnames(op)
diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index 9edf9ba..dac3dd4 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -239,7 +239,7 @@ function apply_gate_bp!(
         dest::AbstractITensorNetwork, op::AbstractITensor,
         state::AbstractITensorNetwork, env; kwargs...
     )
-    vertices = supportof(state, op)
+    vertices = operator_support(state, op)
 
     isempty(vertices) && throw(
         ArgumentError("operator shares no indices with the tensor network")
diff --git a/test/test_tensornetwork.jl b/test/test_tensornetwork.jl
index 29c600a..3443f4a 100644
--- a/test/test_tensornetwork.jl
+++ b/test/test_tensornetwork.jl
@@ -4,7 +4,7 @@ using Graphs: add_edge!, add_vertex!, dst, edges, edgetype, has_edge, has_vertex
     is_directed, ne, nv, rem_edge!, rem_vertex!, src, vertices
 using ITensorBase: Index, LazyITensor, inds, operator
 using ITensorNetworksNext: ITensorNetwork, has_ind, linkaxes, linkinds, linknames, siteaxes,
-    siteinds, sitenames, supportof, tensornetwork
+    siteinds, sitenames, operator_support, tensornetwork
 using NamedGraphs: convert_vertextype, incident_edges, named_grid, named_path_graph,
     similar_graph, subgraph, vertextype
 using Test: @test, @test_throws, @testset
@@ -125,7 +125,7 @@ using Test: @test, @test_throws, @testset
         @test sitenames(tn, 3) == [s[3].name]
     end
 
-    @testset "`supportof`" begin
+    @testset "`operator_support`" begin
         g = named_path_graph(3)
         l = Dict(e => Index(2) for e in edges(g))
         l = merge(l, Dict(reverse(e) => l[e] for e in edges(g)))
@@ -136,21 +136,21 @@ using Test: @test, @test_throws, @testset
         end
 
         o1 = operator(randn(2, 2), (Index(2),), (s[2],))
-        @test supportof(tn, o1) == Set([2])
+        @test operator_support(tn, o1) == Set([2])
 
         o12 = operator(randn(2, 2, 2, 2), (Index(2), Index(2)), (s[1], s[2]))
-        @test supportof(tn, o12) == Set([1, 2])
+        @test operator_support(tn, o12) == Set([1, 2])
 
         # An input name that no tensor in the network carries contributes no vertex.
         o_absent = operator(randn(2, 2), (Index(2),), (Index(2),))
-        @test isempty(supportof(tn, o_absent))
+        @test isempty(operator_support(tn, o_absent))
 
         o_partial = operator(randn(2, 2, 2, 2), (Index(2), Index(2)), (s[3], Index(2)))
-        @test supportof(tn, o_partial) == Set([3])
+        @test operator_support(tn, o_partial) == Set([3])
 
         # A link index is carried by both endpoints of its edge.
         o_link = operator(randn(2, 2), (Index(2),), (l[first(edges(g))],))
-        @test_throws ArgumentError supportof(tn, o_link)
+        @test_throws ArgumentError operator_support(tn, o_link)
     end
 
     @testset "`subgraph`" begin

From 5b277f368e0deacb0d156bb00bb9d8028f7a3317 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 14 Sep 2026 13:34:52 -0400
Subject: [PATCH 11/16] Fix variable names that refer to `inds` instead of the
 correct `names.`

---
 src/tensornetwork.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/tensornetwork.jl b/src/tensornetwork.jl
index 964b682..9eb2218 100644
--- a/src/tensornetwork.jl
+++ b/src/tensornetwork.jl
@@ -154,16 +154,16 @@ function update_tensornetwork_metadata!(tn, vertex, tensor)
     return tn
 end
 
-function update_tensornetwork_metadata!(tn, vertex, oldinds, newinds)
+function update_tensornetwork_metadata!(tn, vertex, oldnames, newnames)
     # Only have to deal with the indices that aren't shared.
-    for name in symdiff(oldinds, newinds)
-        if name in oldinds
+    for name in symdiff(oldnames, newnames)
+        if name in oldnames
             delete_ind_edge!(tn, name)
             delete_ind_vertex!(tn, name, vertex)
             continue
         end
 
-        # Now `ind` must be a new index that's not in `oldinds`
+        # Now `name` must be a new index that's not in `oldinds`
 
         vertex_list = get!(tn.dimname_vertices, name, Set())
         if length(vertex_list) > 1

From 112b78c3cd0e0908a748e6e21eedb30efb37f055 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 14 Sep 2026 13:38:57 -0400
Subject: [PATCH 12/16] Refactor and rename `sum_log_scalars` to `sumlog`.

---
 src/beliefpropagation/messagecache.jl | 22 ++++++++++++++++------
 1 file changed, 16 insertions(+), 6 deletions(-)

diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index 8100a86..d9f3b28 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/src/beliefpropagation/messagecache.jl
@@ -175,11 +175,21 @@ function region_scalar(factors, messages, region)
     return mapreduce(vertex -> vertex_scalar(factors, messages, vertex), *, region)
 end
 
-function sum_log_scalars(terms)
-    if any(t -> real(t) < 0, terms)
-        terms = complex.(terms)
-    end
-    return sum(log.(terms))
+# (log|∏terms|, sign(∏terms))
+function sumlogabs(terms)
+
+    T = typeof(first(terms))
+
+    return mapreduce(
+        t -> (log(abs(t)), sign(t)),
+        ((d1, s1), (d2, s2)) -> (d1 + d2, s1 * s2),
+        terms; init = (zero(float(real(T))), one(T))
+    )
+end
+
+function sumlog(terms)
+    d, s = sumlogabs(terms)
+    return s isa Real && s > 0 ? d : d + log(complex(s))
 end
 
 # We need a graph structure here, so assume `factors` is a graph.
@@ -191,7 +201,7 @@ function bethe_free_energy(factors, messages)
         return -Inf
     end
 
-    return sum_log_scalars(numerator_terms) - sum_log_scalars(denominator_terms)
+    return sumlog(numerator_terms) - sumlog(denominator_terms)
 end
 
 # ===================================== NormNetwork ====================================== #

From 0bc2b8f9e32413f9c0c5977b90b5230534d64001 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 14 Sep 2026 13:42:24 -0400
Subject: [PATCH 13/16] Replace `check_incoming_dimnames` with `check_input`
 and function dispatch

Convention from `MatrixAlgebraKit`.
---
 src/tensornetwork.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/tensornetwork.jl b/src/tensornetwork.jl
index 9eb2218..1cb3cbc 100644
--- a/src/tensornetwork.jl
+++ b/src/tensornetwork.jl
@@ -116,7 +116,7 @@ DataGraphs.is_edge_assigned(::ITensorNetwork, _edge) = false
 
 DataGraphs.get_vertex_data(tn::ITensorNetwork, v) = tn.tensors[v]
 
-function check_incoming_dimnames(tn, tensor, vertex)
+function check_input(::typeof(set_vertex_data!), tn, tensor, vertex)
     for name in dimnames(tensor)
         vertices = get(tn.dimname_vertices, name, Set())
         if length(setdiff(vertices, Set([vertex]))) > 1
@@ -131,7 +131,7 @@ function check_incoming_dimnames(tn, tensor, vertex)
 end
 
 function DataGraphs.insert_vertex_data!(tn::ITensorNetwork, vertex, tensor)
-    check_incoming_dimnames(tn, tensor, vertex)
+    check_input(set_vertex_data!, tn, tensor, vertex)
     add_vertex!(tn.underlying_graph, vertex)
     update_tensornetwork_metadata!(tn, vertex, tensor)
     insert!(tn.tensors, vertex, tensor)
@@ -139,7 +139,7 @@ function DataGraphs.insert_vertex_data!(tn::ITensorNetwork, vertex, tensor)
 end
 
 function DataGraphs.set_vertex_data!(tn::ITensorNetwork, tensor, vertex)
-    check_incoming_dimnames(tn, tensor, vertex)
+    check_input(set_vertex_data!, tn, tensor, vertex)
     update_tensornetwork_metadata!(tn, vertex, tensor)
     set!(tn.tensors, vertex, tensor)
     return tn

From 4482a48a9633b7a1d61b72eac9aedb8379f43868 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 14 Sep 2026 13:42:52 -0400
Subject: [PATCH 14/16] Formatting

---
 src/apply/apply_operators.jl | 4 ++--
 test/test_tensornetwork.jl   | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl
index dac3dd4..f82528d 100644
--- a/src/apply/apply_operators.jl
+++ b/src/apply/apply_operators.jl
@@ -2,8 +2,8 @@ using .AlgorithmsInterfaceExtensions: AlgorithmsInterfaceExtensions as AIE
 using AlgorithmsInterface: AlgorithmsInterface as AI
 using Base: @kwdef
 using Graphs: dst, src, vertices
-using ITensorBase: ITensorBase as ITB, AbstractITensor, apply, dimnames, inputnames, operator,
-    outputnames, replacedimnames
+using ITensorBase: ITensorBase as ITB, AbstractITensor, apply, dimnames, inputnames,
+    operator, outputnames, replacedimnames
 using LinearAlgebra: norm, normalize!
 using MatrixAlgebraKit: eigh_full, project_hermitian, qr_compact, svd_trunc
 using NamedGraphs: boundary_edges
diff --git a/test/test_tensornetwork.jl b/test/test_tensornetwork.jl
index 3443f4a..b9e318b 100644
--- a/test/test_tensornetwork.jl
+++ b/test/test_tensornetwork.jl
@@ -3,8 +3,8 @@ using DataGraphs:
 using Graphs: add_edge!, add_vertex!, dst, edges, edgetype, has_edge, has_vertex,
     is_directed, ne, nv, rem_edge!, rem_vertex!, src, vertices
 using ITensorBase: Index, LazyITensor, inds, operator
-using ITensorNetworksNext: ITensorNetwork, has_ind, linkaxes, linkinds, linknames, siteaxes,
-    siteinds, sitenames, operator_support, tensornetwork
+using ITensorNetworksNext: ITensorNetwork, has_ind, linkaxes, linkinds, linknames,
+    operator_support, siteaxes, siteinds, sitenames, tensornetwork
 using NamedGraphs: convert_vertextype, incident_edges, named_grid, named_path_graph,
     similar_graph, subgraph, vertextype
 using Test: @test, @test_throws, @testset

From 48534f9423179d2d3aa81b367b2bc88fbe6db749 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 14 Sep 2026 14:32:40 -0400
Subject: [PATCH 15/16] `finalize_substate!` now takes both solve and subsolve
 objects.

---
 .../AlgorithmsInterfaceExtensions.jl                      | 8 ++++++--
 src/beliefpropagation/messagecache.jl                     | 1 -
 test/test_algorithmsinterfaceextensions.jl                | 2 +-
 3 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl b/src/AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl
index 917fef9..7863604 100644
--- a/src/AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl
+++ b/src/AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl
@@ -20,7 +20,8 @@ function initialize_subsolve(
 end
 
 function finalize_substate!(
-        problem::AI.Problem, algorithm::AI.Algorithm, state::AI.State, substate::AI.State
+        _problem::AI.Problem, _algorithm::AI.Algorithm, state::AI.State,
+        _subproblem::AI.Problem, _subalgorithm::AI.Algorithm, substate::AI.State
     )
     state.iterate = substate.iterate
     return state
@@ -29,7 +30,10 @@ end
 function AI.step!(problem::AI.Problem, algorithm::NestedAlgorithm, state::AI.State)
     subproblem, subalgorithm, substate = initialize_subsolve(problem, algorithm, state)
     AI.solve!(subproblem, subalgorithm, substate)
-    finalize_substate!(subproblem, subalgorithm, substate, state)
+    finalize_substate!(
+        problem, algorithm, state,
+        subproblem, subalgorithm, substate
+    )
     return state
 end
 
diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index d9f3b28..998cd18 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/src/beliefpropagation/messagecache.jl
@@ -177,7 +177,6 @@ end
 
 # (log|∏terms|, sign(∏terms))
 function sumlogabs(terms)
-
     T = typeof(first(terms))
 
     return mapreduce(
diff --git a/test/test_algorithmsinterfaceextensions.jl b/test/test_algorithmsinterfaceextensions.jl
index 290ebb8..e6b816a 100644
--- a/test/test_algorithmsinterfaceextensions.jl
+++ b/test/test_algorithmsinterfaceextensions.jl
@@ -116,7 +116,7 @@ end
         # `finalize_substate!` copies the substate's iterate back into the
         # parent state.
         substate = AI.initialize_state(problem, algorithm; iterate = [42.0])
-        AIE.finalize_substate!(problem, algorithm, state, substate)
+        AIE.finalize_substate!(problem, algorithm, state, problem, algorithm, substate)
         @test state.iterate == [42.0]
     end
 

From 1ed3f96bfea5829be4b5d014b858b8ea39c8c329 Mon Sep 17 00:00:00 2001
From: Jack Dunham <jdunham@flatironinstitute.org>
Date: Mon, 14 Sep 2026 17:52:17 -0400
Subject: [PATCH 16/16] Fix type stabilty in `vertex/edge_scalars`; `sumlogabs`
 now infers eltype from param

---
 src/beliefpropagation/messagecache.jl | 19 ++++++++-----------
 1 file changed, 8 insertions(+), 11 deletions(-)

diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl
index 998cd18..91f6aa2 100644
--- a/src/beliefpropagation/messagecache.jl
+++ b/src/beliefpropagation/messagecache.jl
@@ -142,8 +142,10 @@ vertex_scalars(factors, messages) = vertex_scalars(factors, messages, keys(facto
 function vertex_scalars(factors::AbstractGraph, messages)
     return vertex_scalars(factors, messages, vertices(factors))
 end
+# `vertex_scalar` reads a number out of an `ITensor`, whose array field is untyped, so `map` would
+# give element type `Any`; collecting the values instead picks up the type they actually have.
 function vertex_scalars(factors, messages, vertices)
-    return map(v -> vertex_scalar(factors, messages, v), vertices)
+    return [vertex_scalar(factors, messages, v) for v in vertices]
 end
 
 function edge_scalar(cache, edge)
@@ -153,22 +155,17 @@ end
 edge_scalars(cache) = edge_scalars(cache, keys(cache))
 
 function edge_scalars(cache, edges)
-    processed = Set{eltype(edges)}()
-
-    T = Base.promote_op(edge_scalar, typeof(cache), eltype(edges))
-
-    scalars = T[]
+    unique_edges = Indices{eltype(edges)}()
 
     # Ignore repeated edges and their reverses.
     for e in edges
-        if e in processed || reverse(e) in processed
+        if e in unique_edges || reverse(e) in unique_edges
             continue
         end
-        push!(processed, e)
-        push!(scalars, edge_scalar(cache, e))
+        insert!(unique_edges, e)
     end
 
-    return scalars
+    return [edge_scalar(cache, e) for e in unique_edges]
 end
 
 function region_scalar(factors, messages, region)
@@ -177,7 +174,7 @@ end
 
 # (log|∏terms|, sign(∏terms))
 function sumlogabs(terms)
-    T = typeof(first(terms))
+    T = eltype(terms)
 
     return mapreduce(
         t -> (log(abs(t)), sign(t)),