Add ITensorNetworkOperator and QuadraticFormNetwork data structures

Project.toml

@@ -1,6 +1,6 @@
 name = "ITensorNetworksNext"
 uuid = "302f2e75-49f0-4526-aef7-d8ba550cb06c"
-version = "0.11.1"
+version = "0.11.2"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
@@ -23,6 +23,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
+UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [weakdeps]
 OMEinsumContractionOrders = "6f22d1fd-8eed-4bb7-9776-e7d684900715"
@@ -48,4 +49,5 @@ Random = "1.10"
 SimpleTraits = "0.9.5"
 SplitApplyCombine = "1.2.3"
 TensorAlgebra = "0.23.1"
+UUIDs = "1.10"
 julia = "1.10"

src/ITensorNetworkGenerators/ising_network.jl

@@ -1,6 +1,6 @@
 using ..ITensorNetworksNext
 using Graphs: degree, dst, edges, src
-using ITensorBase: NamedTensor, name, uniquename
+using ITensorBase: NamedTensor, inds, name, uniquename
 using LinearAlgebra: Diagonal, eigen
 using NamedGraphs: vertextype
 
@@ -33,7 +33,7 @@ function ising_network(
     fp(e) = get(() -> l̃[reverse(e)], l̃, e)
     tn = delta_network(fp, elt, g)
     for v in sz_vertices
-        tn[v] = diagonaltensor(elt[1, -1], axes(tn[v]))
+        tn[v] = diagonaltensor(elt[1, -1], Tuple(inds(tn, v)))
     end
 
     for e in edges(g)

src/ITensorNetworksNext.jl

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

src/abstracttensornetwork.jl

@@ -3,8 +3,8 @@ using DataGraphs: DataGraphs, AbstractDataGraph, AbstractVertexDataGraph, edge_d
     set_vertex_data!, underlying_graph, underlying_graph_type, vertex_data
 using Graphs: Graphs, AbstractEdge, AbstractGraph, add_edge!, add_vertex!, dst, edges,
     edgetype, ne, neighbors, nv, rem_edge!, src, vertices
-using ITensorBase: ITensorOperator, NamedUnitRange, inds, inputnames, name, names, nametype,
-    prime, uniquename, unnamedtype
+using ITensorBase: ITensorBase, ITensorOperator, NamedUnitRange, inds, inputnames, name,
+    names, nametype, prime, uniquename, unnamedtype
 using LinearAlgebra: LinearAlgebra
 using MacroTools: @capture
 using NamedGraphs:
@@ -37,57 +37,63 @@ end
 
 # ====================================== interface ======================================= #
 
+ITensorBase.inds(tn::AbstractITensorNetwork, v) = inds(tn[v])
+ITensorBase.names(tn::AbstractITensorNetwork, v) = name.(inds(tn, v))
+
 linkinds(tn::AbstractGraph, edge::Pair) = linkinds(tn, edgetype(tn)(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]
+    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)
     ln = linknames(tn, edge)
-    return [ax for ax in axes(tn[src(edge)]) if name(ax) in ln]
+    return [ax for ax in inds(tn, src(edge)) if name(ax) in ln]
 end
 function linknames(tn::AbstractGraph, edge::Pair)
     return linknames(tn, edgetype(tn)(edge))
 end
 function linknames(tn::AbstractGraph, edge::AbstractEdge)
-    return names(tn[src(edge)]) ∩ names(tn[dst(edge)])
+    return names(tn, src(edge)) ∩ names(tn, dst(edge))
+end
+function linknames(tn::AbstractGraph)
+    return unique(Iterators.flatten(linknames(tn, e) for e in edges(tn)))
 end
 
 function siteinds(tn::AbstractGraph, v)
-    s = inds(tn[v])
+    s = inds(tn, v)
     for v′ in neighbors(tn, v)
-        s = setdiff(s, inds(tn[v′]))
+        s = setdiff(s, inds(tn, v′))
     end
     return s
 end
 function siteaxes(tn::AbstractGraph, v)
-    s = axes(tn[v])
+    s = Tuple(inds(tn, v))
     for v′ in neighbors(tn, v)
-        s = setdiff(s, axes(tn[v′]))
+        s = setdiff(s, inds(tn, v′))
     end
     return s
 end
 function sitenames(tn::AbstractGraph, v)
-    s = names(tn[v])
+    s = names(tn, v)
     for v′ in neighbors(tn, v)
-        s = setdiff(s, names(tn[v′]))
+        s = setdiff(s, names(tn, v′))
     end
     return s
 end
 
 # Return the non-link vertices associated with an dim name
 function dimnamevertices(tn::AbstractGraph, name)
-    sites = vertextype(tn)[]
+    sites = Set{vertextype(tn)}()
 
     for v in vertices(tn)
-        if name ∈ names(tn[v])
+        if name ∈ names(tn, v)
             push!(sites, v)
         end
     end
@@ -97,7 +103,7 @@ end
 
 function has_dimname(tn::AbstractGraph, name)
     for v in vertices(tn)
-        if name ∈ names(tn[v])
+        if name ∈ names(tn, v)
             return true
         end
     end
@@ -107,7 +113,7 @@ end
 has_ind(tn::AbstractGraph, ind) = has_dimname(tn, name(ind))
 
 function insertlink!(tn::AbstractGraph, e)
-    T = eltype(inds(tn[src(e)]))
+    T = eltype(inds(tn, src(e)))
 
     linkind = NamedUnitRange(trivialrange(unnamedtype(T)), uniquename(nametype(T)))
 
@@ -121,9 +127,15 @@ function insertlink!(tn::AbstractGraph, e)
 end
 
 function operator_support(tn::AbstractGraph, op::ITensorOperator)
+    return operator_support_names(tn, inputnames(op))
+end
+
+# Shared with the `ITensorNetworkOperator` method, which is defined alongside that type
+# because it is not yet known at this point in the include order.
+function operator_support_names(tn::AbstractGraph, opnames)
     support = Indices{vertextype(tn)}()
 
-    for name in inputnames(op)
+    for name in opnames
         vertices = dimnamevertices(tn, name)
 
         if length(vertices) > 1

src/beliefpropagation/beliefpropagation.jl

@@ -246,6 +246,12 @@ factor_tensors(factors, vertex) = [factors[vertex]]
 function factor_tensors(factors::NormNetwork, vertex)
     return [kettensor(factors, vertex), bratensor(factors, vertex)]
 end
+function factor_tensors(factors::QuadraticFormNetwork, vertex)
+    return [
+        kettensor(factors, vertex), operatortensor(factors, vertex),
+        bratensor(factors, vertex),
+    ]
+end
 
 # Contract the incoming messages into the source factor to form the (unnormalized) new message on
 # `edge`.
@@ -272,9 +278,9 @@ end
 # the message is positive semidefinite and its trace is a positive normalization.
 function message_update!(algorithm::SimpleMessageUpdate, cache, factors::NormNetwork, edge)
     new_tensor = updated_message(algorithm, cache, factors, edge)
-    new_message = operator(
-        new_tensor, linknames(BraView(factors), edge), linknames(KetView(factors), edge)
-    )
+    branames = linknames(branetwork(factors), edge)
+    ketnames = linknames(ketnetwork(factors), edge)
+    new_message = operator(new_tensor, branames, ketnames)
     if algorithm.normalize
         message_norm = tr(new_message)
         iszero(message_norm) || (new_message /= message_norm)

src/beliefpropagation/messagecache.jl

@@ -191,16 +191,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