Add apply_operator(s) for gate application with BP simple-update backend (#114)

## Summary

Adds `apply_operator` and `apply_operators` for applying gates to a
`TensorNetwork`, with a belief-propagation simple-update backend
(`BPApplyGate`) and a `MessageCache`-based BP environment for the
double-layer norm network `⟨tn|tn⟩`. Built on the operator-construction
primitives in ITensor/TensorAlgebra.jl#177 and
ITensor/NamedDimsArrays.jl#229.

Also renames `insert_trivial_link!` to `insertlink!` and routes its
underlying-range construction through `TensorAlgebra.trivialrange`. Adds
`randlinknames`, which copies a tensor network, renaming every shared
link with a fresh `randname`, used to contract networks such as `⟨ψ|ψ⟩`.

---------

Co-authored-by: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: mtfishman

Project.toml

@@ -1,6 +1,6 @@
 name = "ITensorNetworksNext"
 uuid = "302f2e75-49f0-4526-aef7-d8ba550cb06c"
-version = "0.4.3"
+version = "0.4.4"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
@@ -19,8 +19,10 @@ FunctionImplementations = "7c7cc465-9c6a-495f-bdd1-f42428e86d0c"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde"
 NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
@@ -47,11 +49,13 @@ FunctionImplementations = "0.4.1"
 Graphs = "1.13.1"
 LinearAlgebra = "1.10"
 MacroTools = "0.5.16"
-NamedDimsArrays = "0.14.3, 0.15"
+MatrixAlgebraKit = "0.6"
+NamedDimsArrays = "0.15.8"
 NamedGraphs = "0.11"
+Random = "1.10"
 SimpleTraits = "0.9.5"
 SplitApplyCombine = "1.2.3"
-TensorAlgebra = "0.9.2"
+TensorAlgebra = "0.9.5"
 TensorOperations = "5.3.1"
 TermInterface = "2"
 TypeParameterAccessors = "0.4.4"

src/AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl

@@ -1,6 +1,6 @@
 module AlgorithmsInterfaceExtensions
 
-import AlgorithmsInterface as AI
+using AlgorithmsInterface: AlgorithmsInterface as AI
 
 # ============================ NestedAlgorithm =============================================
 
@@ -52,44 +52,6 @@ function Base.propertynames(state::NestedState)
     return (fieldnames(typeof(state))..., :iterate)
 end
 
-# ============================ select_algorithm / default_algorithm ========================
-
-# Like `MatrixAlgebraKit.select_algorithm` / `default_algorithm`, but
-# selection-relevant inputs are packed into an `args` tuple so the value
-# and type domains stay disjoint: `(1.2,)` vs `Tuple{Float64}`. Strategy
-# types subtype `AbstractAlgorithm` so the passthrough overload is generic.
-abstract type AbstractAlgorithm end
-
-function default_algorithm(f, ::Type{Args}; kwargs...) where {Args <: Tuple}
-    return throw(MethodError(default_algorithm, (f, Args)))
-end
-function default_algorithm(f, args::Tuple; kwargs...)
-    return default_algorithm(f, typeof(args); kwargs...)
-end
-
-function select_algorithm(f, alg, args::Tuple; kwargs...)
-    return select_algorithm(f, alg, typeof(args); kwargs...)
-end
-function select_algorithm(f, ::Nothing, ::Type{Args}; kwargs...) where {Args <: Tuple}
-    return default_algorithm(f, Args; kwargs...)
-end
-function select_algorithm(f, alg::NamedTuple, ::Type{Args}; kwargs...) where {Args <: Tuple}
-    isempty(kwargs) || throw(
-        ArgumentError(
-            "Additional keyword arguments are not allowed when `alg` is a `NamedTuple`."
-        )
-    )
-    return default_algorithm(f, Args; alg...)
-end
-function select_algorithm(f, alg::AbstractAlgorithm, ::Type{<:Tuple}; kwargs...)
-    isempty(kwargs) || throw(
-        ArgumentError(
-            "Additional keyword arguments are not allowed when `alg` is an `AbstractAlgorithm` instance."
-        )
-    )
-    return alg
-end
-
 # ============================ StopWhenConverged ===========================================
 
 # Stopping criterion that fires once `iterate_diff(iterate, previous_iterate) < tol`.

src/ITensorNetworksNext.jl

@@ -1,11 +1,14 @@
 module ITensorNetworksNext
 
-# Imported as a name only so that the `[compat] TensorAlgebra = "0.9.2"` floor
-# (needed for a `bipermutedimsopadd!` fix in `TensorAlgebra` 0.9.2 that affects
-# `NamedDimsArrays`-mediated tensor multiplication) isn't reported as a stale
-# dependency by Aqua.
-using TensorAlgebra: TensorAlgebra
+if VERSION >= v"1.11.0-DEV.469"
+    eval(
+        Meta.parse(
+            "public apply_operator, apply_operators, beliefpropagation_normnetwork, identity_norm_message_env, normnetwork, norm_message_env, ones_norm_message_env, rand_norm_message_env, randn_norm_message_env, similar_norm_message_env"
+        )
+    )
+end
 
+include("select_algorithm.jl")
 include("AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl")
 include("LazyNamedDimsArrays/LazyNamedDimsArrays.jl")
 include("abstracttensornetwork.jl")
@@ -15,5 +18,8 @@ include("contract_network.jl")
 
 include("beliefpropagation/messagecache.jl")
 include("beliefpropagation/beliefpropagation.jl")
+include("beliefpropagation/normnetwork.jl")
+
+include("apply/apply_operators.jl")
 
 end

src/abstracttensornetwork.jl

@@ -63,13 +63,20 @@ function Adapt.adapt_structure(to, tn::AbstractTensorNetwork)
 end
 
 linkinds(tn::AbstractGraph, edge::Pair) = linkinds(tn, edgetype(tn)(edge))
-linkinds(tn::AbstractGraph, edge::AbstractEdge) = inds(tn[src(edge)]) ∩ inds(tn[dst(edge)])
+# Pick the link indices from the `src` side, identified by name match with `dst`.
+# A range-strict intersection (`inds(src) ∩ inds(dst)`) would drop graded links
+# whose two endpoints carry dual-related ranges.
+function linkinds(tn::AbstractGraph, edge::AbstractEdge)
+    ln = linknames(tn, edge)
+    return [i for i in inds(tn[src(edge)]) if name(i) in ln]
+end
 
 function linkaxes(tn::AbstractGraph, edge::Pair)
     return linkaxes(tn, edgetype(tn)(edge))
 end
 function linkaxes(tn::AbstractGraph, edge::AbstractEdge)
-    return axes(tn[src(edge)]) ∩ axes(tn[dst(edge)])
+    ln = linknames(tn, edge)
+    return [ax for ax in axes(tn[src(edge)]) if name(ax) in ln]
 end
 function linknames(tn::AbstractGraph, edge::Pair)
     return linknames(tn, edgetype(tn)(edge))
@@ -143,7 +150,7 @@ function add_missing_edges!(tn::AbstractGraph, v)
     for v′ in vertices(tn)
         if v ≠ v′
             e = v => v′
-            if !isempty(linkinds(tn, e))
+            if !isempty(linknames(tn, e))
                 add_edge!(tn, e)
             end
         end
@@ -162,37 +169,40 @@ end
 function fix_edges!(tn::AbstractGraph, v)
     for e in incident_edges(tn, v)
         # Remove an edge if there is no index on that edge.
-        if isempty(linkinds(tn, e))
+        if isempty(linknames(tn, e))
             rem_edge!(tn, e)
         end
     end
     add_missing_edges!(tn, v)
     return tn
 end
 
-# Customization point.
-using NamedDimsArrays: AbstractNamedUnitRange, namedunitrange, nametype, randname
-function trivial_unitrange(type::Type{<:AbstractUnitRange})
-    return Base.oneto(one(eltype(type)))
-end
-function rand_trivial_namedunitrange(
-        ::Type{<:AbstractNamedUnitRange{<:Any, R, N}}
-    ) where {R, N}
-    return namedunitrange(trivial_unitrange(R), randname(N))
-end
-
-dag(x) = x
-
-function insert_trivial_link!(tn, e)
+using NamedDimsArrays: denamedtype, named, nametype, randname
+using TensorAlgebra: trivialrange
+function insertlink!(tn, e)
     add_edge!(tn, e)
-    l = rand_trivial_namedunitrange(eltype(inds(tn[src(e)])))
-    x = similar(tn[src(e)], (l,))
-    x[1] = 1
+    T = eltype(inds(tn[src(e)]))
+    l = named(trivialrange(denamedtype(T)), randname(nametype(T)))
+    x = fill!(similar(tn[src(e)], (l,)), one(eltype(tn[src(e)])))
     @preserve_graph tn[src(e)] = tn[src(e)] * x
-    @preserve_graph tn[dst(e)] = tn[dst(e)] * dag(x)
+    @preserve_graph tn[dst(e)] = tn[dst(e)] * conj(x)
     return tn
 end
 
+using NamedDimsArrays: replacedimnames
+function randlinknames(tn)
+    new_tn = copy(tn)
+    for e in edges(new_tn)
+        u, v = src(e), dst(e)
+        for n in intersect(dimnames(new_tn[u]), dimnames(new_tn[v]))
+            n′ = randname(n)
+            new_tn[u] = replacedimnames(new_tn[u], n => n′)
+            new_tn[v] = replacedimnames(new_tn[v], n => n′)
+        end
+    end
+    return new_tn
+end
+
 function Base.setindex!(tn::AbstractTensorNetwork, value, v)
     @preserve_graph tn[v] = value
     fix_edges!(tn, v)