diff --git a/Project.toml b/Project.toml index b5c5900..8ba3dca 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ITensorNetworksNext" uuid = "302f2e75-49f0-4526-aef7-d8ba550cb06c" -version = "0.7.1" +version = "0.8.0" authors = ["ITensor developers and contributors"] [workspace] @@ -27,7 +27,7 @@ TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" Adapt = "4.3" AlgorithmsInterface = "0.1" Combinatorics = "1" -DataGraphs = "0.4" +DataGraphs = "0.5" Dictionaries = "0.4.5" Graphs = "1.13.1" ITensorBase = "0.8" diff --git a/docs/Project.toml b/docs/Project.toml index f01c94b..bae8d7c 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -10,5 +10,5 @@ path = ".." [compat] Documenter = "1" ITensorFormatter = "0.2.27" -ITensorNetworksNext = "0.7" +ITensorNetworksNext = "0.8" Literate = "2" diff --git a/examples/Project.toml b/examples/Project.toml index ad5bc0b..2d44b8a 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -5,4 +5,4 @@ ITensorNetworksNext = "302f2e75-49f0-4526-aef7-d8ba550cb06c" path = ".." [compat] -ITensorNetworksNext = "0.7" +ITensorNetworksNext = "0.8" diff --git a/src/ITensorNetworkGenerators/delta_network.jl b/src/ITensorNetworkGenerators/delta_network.jl index 7b1370a..1721199 100644 --- a/src/ITensorNetworkGenerators/delta_network.jl +++ b/src/ITensorNetworkGenerators/delta_network.jl @@ -1,5 +1,5 @@ -using ..ITensorNetworksNext: ITensorNetwork -using Graphs: AbstractGraph +using ..ITensorNetworksNext: tensornetwork +using Graphs: AbstractGraph, vertices using ITensorBase: NamedUnitRange, name, nameddims, unnamed using NamedGraphs.GraphsExtensions: incident_edges @@ -42,10 +42,11 @@ on each vertex. Link dimensions are defined using the function `f(e)` that shoul edge `e` as an input and should output the link index on that edge. """ function delta_network(f, elt::Type, g::AbstractGraph) - return tn = ITensorNetwork(g) do v + tn = tensornetwork(vertices(g)) do v is = Tuple(f.(incident_edges(g, v))) return delta(elt, is) end + return tn end function delta_network(f, g::AbstractGraph) return delta_network(f, Float64, g) diff --git a/src/ITensorNetworkGenerators/ising_network.jl b/src/ITensorNetworkGenerators/ising_network.jl index 3cdf9a2..803adff 100644 --- a/src/ITensorNetworkGenerators/ising_network.jl +++ b/src/ITensorNetworkGenerators/ising_network.jl @@ -1,4 +1,4 @@ -using ..ITensorNetworksNext: @preserve_graph +using ..ITensorNetworksNext using Graphs: degree, dst, edges, src using ITensorBase: apply, name, operator, uniquename using LinearAlgebra: Diagonal, eigen @@ -29,20 +29,23 @@ function ising_network( ) elt = typeof(β) l̃ = Dict(e => uniquename(f(e)) for e in edges(g)) - f̃(e) = get(() -> l̃[reverse(e)], l̃, e) - tn = delta_network(f̃, elt, g) + + 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])) end - for e in edges(tn) + + for e in edges(g) v1 = src(e) v2 = dst(e) deg1 = degree(tn, v1) + deg2 = degree(tn, v2) m = sqrt_ising_bond(β; J, h, deg1, deg2) - t = operator(m, (name(f̃(e)),), (name(f(e)),)) - @preserve_graph tn[v1] = apply(t, tn[v1]) - @preserve_graph tn[v2] = apply(t, tn[v2]) + t = operator(m, (name(fp(e)),), (name(f(e)),)) + tn[v1] = apply(t, tn[v1]) + tn[v2] = apply(t, tn[v2]) end return tn end diff --git a/src/ITensorNetworksNext.jl b/src/ITensorNetworksNext.jl index 92a541f..360fb76 100644 --- a/src/ITensorNetworksNext.jl +++ b/src/ITensorNetworksNext.jl @@ -3,7 +3,7 @@ module ITensorNetworksNext 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" + "public apply_operator, apply_operators" ) ) end @@ -12,6 +12,8 @@ include("select_algorithm.jl") include("AlgorithmsInterfaceExtensions/AlgorithmsInterfaceExtensions.jl") include("abstracttensornetwork.jl") include("tensornetwork.jl") +include("normnetwork.jl") +include("normnetworkview.jl") include("ITensorNetworkGenerators/ITensorNetworkGenerators.jl") include("contract_network.jl") diff --git a/src/abstracttensornetwork.jl b/src/abstracttensornetwork.jl index 49d3a14..40bb29d 100644 --- a/src/abstracttensornetwork.jl +++ b/src/abstracttensornetwork.jl @@ -1,66 +1,54 @@ using Adapt: Adapt, adapt -using DataGraphs: DataGraphs, AbstractDataGraph, edge_data, set_vertex_data!, - underlying_graph, underlying_graph_type, vertex_data +using DataGraphs: DataGraphs, AbstractDataGraph, AbstractVertexDataGraph, edge_data, + set_vertex_data!, underlying_graph, underlying_graph_type, vertex_data 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 +using ITensorBase: dimnames, inds, name, named, nametype, prime, uniquename, unnamedtype using LinearAlgebra: LinearAlgebra using MacroTools: @capture using NamedGraphs.GraphsExtensions: directed_graph, incident_edges, rem_edges!, vertextype using NamedGraphs.OrdinalIndexing: OrdinalSuffixedInteger using NamedGraphs: NamedGraphs, NamedGraph, not_implemented, similar_graph +using TensorAlgebra: trivialrange + +abstract type AbstractITensorNetwork{T, V} <: AbstractVertexDataGraph{T, V} end -abstract type AbstractITensorNetwork{V, VD} <: AbstractDataGraph{V, VD, Nothing} end +# ====================================== Graphs.jl ======================================= # # Need to be careful about removing edges from tensor networks in case there is a bond -Graphs.rem_edge!(::AbstractITensorNetwork, edge) = not_implemented() +Graphs.rem_edge!(::AbstractITensorNetwork, _edge) = not_implemented() -# Graphs.jl overloads function Graphs.weights(graph::AbstractITensorNetwork) V = vertextype(graph) es = Tuple.(edges(graph)) ws = Dictionary{Tuple{V, V}, Float64}(es, undef) for e in edges(graph) - w = log2(dim(commoninds(graph, e))) + w = log2(dim(linkinds(graph, e))) ws[(src(e), dst(e))] = w end return ws end -# Copy -Base.copy(::AbstractITensorNetwork) = not_implemented() - -# Iteration -Base.iterate(tn::AbstractITensorNetwork, args...) = iterate(vertex_data(tn), args...) -Base.keys(tn::AbstractITensorNetwork) = vertices(tn) - -# TODO: This contrasts with the `DataGraphs.AbstractDataGraph` definition, -# where it is defined as the `vertextype`. Does that cause problems or should it be changed? -Base.eltype(tn::AbstractITensorNetwork) = eltype(vertex_data(tn)) - # Overload if needed Graphs.is_directed(::Type{<:AbstractITensorNetwork}) = false -DataGraphs.underlying_graph(::AbstractITensorNetwork) = not_implemented() -function NamedGraphs.vertex_positions(tn::AbstractITensorNetwork) - return NamedGraphs.vertex_positions(underlying_graph(tn)) -end -function NamedGraphs.ordered_vertices(tn::AbstractITensorNetwork) - return NamedGraphs.ordered_vertices(underlying_graph(tn)) +# ==================================== NamedGraphs.jl ==================================== # + +function NamedGraphs.similar_graph(::AbstractITensorNetwork, VD::Type, vertices) + return ITensorNetwork{VD}(undef, collect(vertices)) end -function Adapt.adapt_structure(to, tn::AbstractITensorNetwork) - # TODO: Define and use: - # - # @preserve_graph map_vertex_data(adapt(to), tn) - # - # or just: - # - # @preserve_graph map(adapt(to), tn) - return map_vertex_data_preserve_graph(adapt(to), tn) +# ==================================== DataGraphs.jl ===================================== # + +function DataGraphs.underlying_graph(tn::AbstractITensorNetwork) + ug = NamedGraph(vertices(tn)) + add_edges!(ug, edges(tn)) + return ug end +# ====================================== interface ======================================= # + 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 @@ -106,139 +94,40 @@ function sitenames(tn::AbstractGraph, v) return s end -function setindex_preserve_graph!(tn::AbstractGraph, value, vertex) - set_vertex_data!(tn, value, vertex) - return tn -end +# Return the non-link vertices associated with an dim name +function dimnamevertices(tn::AbstractGraph, name) + sites = vertextype(tn)[] -function is_setindex!_expr(expr::Expr) - return is_assignment_expr(expr) && is_getindex_expr(first(expr.args)) -end -is_setindex!_expr(x) = false -is_getindex_expr(expr::Expr) = (expr.head === :ref) -is_getindex_expr(x) = false -is_assignment_expr(expr::Expr) = (expr.head === :(=)) -is_assignment_expr(expr) = false - -# TODO: Define this in terms of a function mapping -# preserve_graph_function(::typeof(setindex!)) = setindex!_preserve_graph -# preserve_graph_function(::typeof(map_vertex_data)) = map_vertex_data_preserve_graph -# Also allow annotating codeblocks like `@views`. -macro preserve_graph(expr) - if !is_setindex!_expr(expr) - error( - "preserve_graph must be used with setindex! syntax (as @preserve_graph a[i,j,...] = value)" - ) + for v in vertices(tn) + if name ∈ dimnames(tn[v]) + push!(sites, v) + end end - @capture(expr, array_[indices__] = value_) - return :(setindex_preserve_graph!($(esc(array)), $(esc(value)), $(esc.(indices)...))) -end -# Update the graph of the ITensorNetwork `tn` to include -# edges that should exist based on the tensor connectivity. -function add_missing_edges!(tn::AbstractGraph) - foreach(v -> add_missing_edges!(tn, v), vertices(tn)) - return tn + return sites end -# Update the graph of the ITensorNetwork `tn` to include -# edges that should be incident to the vertex `v` -# based on the tensor connectivity. -function add_missing_edges!(tn::AbstractGraph, v) - for v′ in vertices(tn) - if v ≠ v′ - e = v => v′ - if !isempty(linknames(tn, e)) - add_edge!(tn, e) - end +function has_dimname(tn::AbstractGraph, name) + for v in vertices(tn) + if name ∈ dimnames(tn[v]) + return true end end - return tn + return false end -# Fix the edges of the ITensorNetwork `tn` to match -# the tensor connectivity. -function fix_edges!(tn::AbstractGraph) - foreach(v -> fix_edges!(tn, v), vertices(tn)) - return tn -end -# Fix the edges of the ITensorNetwork `tn` to match -# the tensor connectivity at vertex `v`. -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(linknames(tn, e)) - rem_edge!(tn, e) - end - end - add_missing_edges!(tn, v) - return tn -end +has_ind(tn::AbstractGraph, ind) = has_dimname(tn, name(ind)) -using ITensorBase: named, nametype, uniquename, unnamedtype -using TensorAlgebra: trivialrange -function insertlink!(tn, e) - add_edge!(tn, e) +function insertlink!(tn::AbstractGraph, e) T = eltype(inds(tn[src(e)])) - l = named(trivialrange(unnamedtype(T)), uniquename(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)] * conj(x) - return tn -end -using ITensorBase: 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′ = uniquename(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 + linkind = named(trivialrange(unnamedtype(T)), uniquename(nametype(T))) -function Base.setindex!(tn::AbstractITensorNetwork, value, v) - @preserve_graph tn[v] = value - fix_edges!(tn, v) - return tn -end -# Fix ambiguity error. -function Base.setindex!( - graph::AbstractITensorNetwork, - value, - vertex::OrdinalSuffixedInteger - ) - graph[vertices(graph)[vertex]] = value - return graph -end -Base.setindex!(tn::AbstractITensorNetwork, value, edge::AbstractEdge) = not_implemented() -Base.setindex!(tn::AbstractITensorNetwork, value, edge::Pair) = not_implemented() -# Fix ambiguity error. -function Base.setindex!( - tn::AbstractITensorNetwork, - value, - edge::Pair{<:OrdinalSuffixedInteger, <:OrdinalSuffixedInteger} - ) - return not_implemented() -end + x = similar(tn[src(e)], (linkind,)) + fill!(x, true) -function Base.show(io::IO, mime::MIME"text/plain", graph::AbstractITensorNetwork) - println(io, "$(typeof(graph)) with $(nv(graph)) vertices:") - show(io, mime, vertices(graph)) - println(io, "\n") - println(io, "and $(ne(graph)) edge(s):") - for e in edges(graph) - show(io, mime, e) - println(io) - end - println(io) - println(io, "with vertex data:") - show(io, mime, axes.(vertex_data(graph))) - return nothing -end + tn[src(e)] *= x + tn[dst(e)] *= conj(x) -Base.show(io::IO, graph::AbstractITensorNetwork) = show(io, MIME"text/plain"(), graph) + return tn +end diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl index b936035..f18653d 100644 --- a/src/beliefpropagation/beliefpropagation.jl +++ b/src/beliefpropagation/beliefpropagation.jl @@ -5,7 +5,8 @@ using DataGraphs: edge_data using Graphs: AbstractEdge, edges, edgetype, has_edge, vertices using ITensorBase: AbstractITensor using LinearAlgebra: norm, normalize -using NamedGraphs.GraphsExtensions: add_edges!, boundary_edges, subgraph +using NamedGraphs.GraphsExtensions: + add_edges!, boundary_edges, forest_cover_edge_sequence, subgraph using NamedGraphs.PartitionedGraphs: quotientvertices # === Top-level user entry point === @@ -65,7 +66,10 @@ an explicit `AlgorithmsInterface.StoppingCriterion`. `message_update_algorithm` controls how a single message is recomputed from its incoming neighbours. """ -function beliefpropagation( +function beliefpropagation(factors, messages; kwargs...) + return _beliefpropagation(factors, messages; kwargs...) +end +function _beliefpropagation( factors, messages; edges = default_beliefpropagation_edges(factors), stopping_criterion = nothing, diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl index 585cdf6..869ccc5 100644 --- a/src/beliefpropagation/messagecache.jl +++ b/src/beliefpropagation/messagecache.jl @@ -1,15 +1,16 @@ -using DataGraphs: DataGraphs, AbstractDataGraph, edge_data, edge_data_type, - set_vertex_data!, underlying_graph, underlying_graph_type, vertex_data, vertex_data_type +using DataGraphs: DataGraphs, AbstractDataGraph, AbstractEdgeDataGraph, edge_data, + edge_data_type, set_vertex_data!, underlying_graph, underlying_graph_type, vertex_data, + vertex_data_type using Dictionaries: Dictionary, delete!, getindices, set! using Graphs: AbstractGraph, connected_components, is_directed, is_tree using NamedGraphs.GraphsExtensions: IsDirected, boundary_edges, default_root_vertex, directed_graph, forest_cover, in_incident_edges, post_order_dfs_edges, undirected_graph, vertextype using NamedGraphs.PartitionedGraphs: QuotientEdge, QuotientView, quotient_graph -using NamedGraphs: NamedDiGraph, Vertices, convert_vertextype, ordered_vertices, - parent_graph_indices, position_graph, to_graph_index, vertex_positions +using NamedGraphs: AbstractNamedEdge, NamedDiGraph, NamedEdge, Vertices, convert_vertextype, + ordered_vertices, parent_graph_indices, position_graph, to_graph_index, vertex_positions -struct MessageCache{T, V} <: AbstractDataGraph{V, Nothing, T} +struct MessageCache{T, V} <: AbstractEdgeDataGraph{T, V} messages::Dictionary{NamedEdge{V}, T} underlying_graph::NamedDiGraph{V} function MessageCache{T, V}(::UndefInitializer, vertices) where {T, V} @@ -20,19 +21,13 @@ struct MessageCache{T, V} <: AbstractDataGraph{V, Nothing, T} end # single type parameter version of the inner constructor +function MessageCache(::UndefInitializer, vertices) + return MessageCache{Any}(undef, vertices) +end function MessageCache{T}(::UndefInitializer, vertices) where {T} return MessageCache{T, eltype(vertices)}(undef, vertices) end -# compatibility with generic key-val iterables -Base.keytype(c::MessageCache) = keytype(typeof(c)) -Base.keytype(::Type{<:MessageCache{T, V}}) where {T, V} = NamedEdge{V} - -Base.valtype(c::MessageCache) = valtype(typeof(c)) -Base.valtype(::Type{<:MessageCache{T}}) where {T} = T - -Base.keys(cache::MessageCache) = edges(cache) - MessageCache(messages) = MessageCache{valtype(messages)}(messages) function MessageCache{T}(messages) where {T} @@ -54,99 +49,72 @@ end messagecache(pairs) = MessageCache(Dict(pairs)) messagecache(f, edges) = messagecache(edge => f(edge) for edge in edges) -# ================================ NamedGraphs interface ================================= # -function NamedGraphs.add_edge!(c::MessageCache, edge) - add_edge!(c.underlying_graph, edge) +function Graphs.rem_edge!(c::MessageCache, edge) + delete!(c.messages, to_graph_index(c, edge)) + rem_edge!(c.underlying_graph, edge) return c end -function NamedGraphs.rem_edge!(c::MessageCache, edge) - delete!(c.messages, to_graph_index(c, edge)) - rem_edge!(c.underlying_graph, edge) +function Graphs.add_vertex!(c::MessageCache, vertex) + add_edge!(c.underlying_graph, vertex) return c end -# ================================= DataGraphs interface ================================= # +function Graphs.has_edge(c::MessageCache, edge::AbstractNamedEdge) + return has_edge(c.underlying_graph, edge) +end -DataGraphs.underlying_graph(cache::MessageCache) = cache.underlying_graph +# ================================ NamedGraphs interface ================================= # -DataGraphs.is_vertex_assigned(::MessageCache, _) = false -DataGraphs.is_edge_assigned(c::MessageCache, edge) = haskey(c.messages, edge) +function NamedGraphs.similar_graph(::Type{<:MessageCache}, vertices) + return MessageCache(undef, vertices) +end -function DataGraphs.get_edge_data(c::MessageCache, edge::AbstractEdge) - return c.messages[edge] +function NamedGraphs.similar_graph(cache::MessageCache, T::Type) + new_cache = similar_graph(cache, T, vertices(cache)) + add_edges!(new_cache.underlying_graph, edges(cache)) + return new_cache end -function DataGraphs.set_edge_data!(c::MessageCache, val, edge) - return set!(c.messages, edge, val) +function NamedGraphs.similar_graph(::MessageCache, ED::Type, vertices) + return MessageCache{ED}(undef, collect(vertices)) end -Base.copy(cache::MessageCache) = MessageCache(copy(cache.messages)) +# ================================= DataGraphs interface ================================= # -function Base.:(==)(cache1::MessageCache, cache2::MessageCache) - ug1 = cache1.underlying_graph - ug2 = cache2.underlying_graph +DataGraphs.underlying_graph(cache::MessageCache) = cache.underlying_graph - ms1 = cache1.messages - ms2 = cache2.messages +DataGraphs.is_vertex_assigned(::MessageCache, _) = false +DataGraphs.is_edge_assigned(c::MessageCache, edge) = haskey(c.messages, edge) - return (ug1 == ug2 && ms1 == ms2) +DataGraphs.get_edge_data(c::MessageCache, edge::AbstractEdge) = c.messages[edge] +function DataGraphs.set_edge_data!(c::MessageCache, val, edge) + has_edge(c, edge) || add_edge!(c.underlying_graph, edge) + set!(c.messages, edge, val) + return c end -function NamedGraphs.induced_subgraph_from_vertices(cache::MessageCache, subvertices) - # TODO: once we have `subgraph_edges` in `NamedGraphs`, simplify this. - underlying_subgraph, vlist = - Graphs.induced_subgraph(cache.underlying_graph, subvertices) - - assigned = v -> isassigned(cache, v) - - assigned_subedges = Iterators.filter(assigned, edges(underlying_subgraph)) - - messages = getindices(cache.messages, Indices(assigned_subedges)) - - return MessageCache(messages), vlist +function DataGraphs.insert_edge_data!(cache::MessageCache, edge, val) + add_edge!(cache.underlying_graph, edge) + insert!(cache.messages, edge, val) + return cache end -# see: copyto!(dest, src) for analogous behaviour to 2 argument method -# see: copyto!(dest, Rdest::CartesianIndices, src, Rsrc::CartesianIndices) -# for analogous behaviour to 3 argument method. -# TODO: these can be made generic for `AbtractDataGraph` in `DataGraphs.jl` -function copyto!_messagecache( - cache_dst::MessageCache, - cache_src, - inds = nothing - ) - inds = isnothing(inds) ? Indices(keys(cache_src)) : Indices(inds) - view(edge_data(cache_dst), inds) .= view(cache_src, inds) - return cache_dst -end +# =================================== Dictionaries.jl ==================================== # -function Base.copyto!( - cache_dst::MessageCache, - cache_src::AbstractDataGraph, - inds = nothing - ) - copyto!_messagecache(cache_dst, edge_data(cache_src), inds) - return cache_dst -end +Dictionaries.issettable(::MessageCache) = true +Dictionaries.isinsertable(::MessageCache) = true -function Base.copyto!( - cache_dst::MessageCache, - dictionary_src::Dictionary, - inds = nothing - ) - copyto!_messagecache(cache_dst, dictionary_src, inds) - return cache_dst +function Base.map(f, cache::MessageCache) + new_cache = similar_graph(cache, Base.promote_op(f, valtype(cache))) + map!(f, new_cache, cache) + return new_cache end -function Base.copyto!( - cache_dst::MessageCache, - dict_src::Dict, - inds = keys(dict_src) - ) - for key in inds - cache_dst[key] = dict_src[key] +function Base.map!(f, dst::MessageCache, src) + for key in keys(src) + dst[key] = f(src[key]) end - return cache_dst + return dst end # ===================================== contraction ====================================== # @@ -156,14 +124,14 @@ function incoming_messages(cache::AbstractGraph, pair::Pair) return incoming_messages(cache, edge) end function incoming_messages(cache::AbstractGraph, edge::AbstractEdge) - inds = Indices(in_incident_edges(cache, src(edge))) - return getindices(cache, filter(e -> e != reverse(edge), inds)) + dimnames = Indices(in_incident_edges(cache, src(edge))) + return getindices(cache, filter(e -> e != reverse(edge), dimnames)) end # TODO: maybe this should be defined in `DataGraphs`. function incoming_edge_data(cache::AbstractGraph, vertices) - inds = Indices(boundary_edges(cache, vertices; dir = :in)) - return getindices(cache, inds) + dimnames = Indices(boundary_edges(cache, vertices; dir = :in)) + return getindices(cache, dimnames) end function vertex_scalar(factors, messages, vertex; kwargs...) @@ -229,40 +197,3 @@ function bethe_free_energy(factors, messages) return sum(log.(numerator_terms)) - sum(log.(denominator_terms)) end - -# TODO: This needs to go in NamedGraphs.GraphsExtensions -function forest_cover_edge_sequence(gi::AbstractGraph; root_vertex = default_root_vertex) - # All we care about are the edges so the type of the graph doesnt matter - g = similar_graph(NamedGraph, vertices(gi)) - add_edges!(g, edges(gi)) - forests = forest_cover(g) - rv = edgetype(g)[] - for forest in forests - trees = [forest[Vertices(vs)] for vs in connected_components(forest)] - for tree in trees - tree_edges = post_order_dfs_edges(tree, root_vertex(tree)) - push!(rv, vcat(tree_edges, reverse(reverse.(tree_edges)))...) - end - end - return rv -end - -# ======================================= printing ======================================= # - -# TODO: This is the definition for the proposed `DataGraphs.AbstractEdgeDataGraph`. -function Base.show(io::IO, mime::MIME"text/plain", graph::MessageCache) - println(io, "$(typeof(graph)) with $(nv(graph)) vertices:") - show(io, mime, vertices(graph)) - println(io, "\n") - println(io, "and $(ne(graph)) edge(s):") - for e in edges(graph) - show(io, mime, e) - println(io) - end - println(io) - println(io, "with edge data:") - show(io, mime, edge_data(graph)) - return nothing -end - -Base.show(io::IO, graph::MessageCache) = show(io, MIME"text/plain"(), graph) diff --git a/src/beliefpropagation/normnetwork.jl b/src/beliefpropagation/normnetwork.jl index 60ef064..0237d4d 100644 --- a/src/beliefpropagation/normnetwork.jl +++ b/src/beliefpropagation/normnetwork.jl @@ -1,199 +1,56 @@ -using DataGraphs: underlying_graph using Graphs: dst, edges, edgetype, src using ITensorBase: codomainnames, domainnames, name, operator, replacedimnames, similar_operator, state, uniquename, unnamed using NamedGraphs.GraphsExtensions: all_edges, incident_edges -using Random: Random, rand!, randn! +using SplitApplyCombine: mapmany -# === Norm-network environment constructors === -# -# `*_norm_message_env(tn)` builds a `MessageCache` shaped to act as the BP environment -# for the norm network ⟨tn|tn⟩, with each entry filled per the leading verb (`identity`, -# `ones`, `randn`, `rand`). The `_env` suffix is reserved for the high-level -# environment-builder interface; the low-level `MessageCache` / `messagecache(...)` -# constructors are used internally. A parallel `*_norm_ctm_env` family is planned for -# CTMRG environments. +function message_environment(::UndefInitializer, nn::NormNetwork) + messages = mapmany(vertices(nn)) do vertex + return map(in_incident_edges(nn, vertex)) do edge + braview = BraView(nn) + ketview = KetView(nn) -""" - similar_norm_message_env(tn) -> MessageCache + ketnames = linknames(KetView(nn), edge) -Allocate a BP environment for the norm network ⟨tn|tn⟩ with **undefined** message data: -one square operator message per directed edge of `tn` (both directions on every -undirected edge). On each undirected edge the two directions share the same ket-side -names (the link axes from `tn`) and the same fresh `uniquename`-generated bra-side names, -with the codomain and domain swapped between the two directions — so `env[v1=>v2]` and -`env[v2=>v1]` contract directly with each other (matching names, dual axes) for -bond-marginal computations. Element type and backend are inherited from the factor -tensors of `tn` via `Base.similar`. + brainds = linkinds(braview, edge) + branames = name.(brainds) + braaxis = unnamed.(brainds) -Used internally by [`norm_message_env`](@ref) and the filled environment constructors -([`identity_norm_message_env`](@ref), [`ones_norm_message_env`](@ref), -[`randn_norm_message_env`](@ref), [`rand_norm_message_env`](@ref)). Use it directly to -construct environments with custom message data, e.g. by mutating each entry after -allocation. -""" -function similar_norm_message_env(tn) - pairs = [] - for e in edges(tn) - v1, v2 = src(e), dst(e) - ket_axes = linkinds(tn, e) - ket_names = name.(ket_axes) - unnamed_axes = unnamed.(ket_axes) - bra_names = uniquename.(ket_names) - # Message axes are dual to the link they contract against in the factor. - push!( - pairs, - edgetype(tn)(v1, v2) => - similar_operator(tn[v1], conj.(unnamed_axes), bra_names, ket_names) - ) - push!( - pairs, - edgetype(tn)(v2, v1) => - similar_operator(tn[v2], unnamed_axes, bra_names, ket_names) - ) - end - return messagecache(pairs) -end - -""" - norm_message_env(f, tn) -> MessageCache - -Allocate a norm-network BP environment via [`similar_norm_message_env`](@ref) and apply -`f` to each operator-message entry. Shared building block for the filled-environment -constructors. -""" -function norm_message_env(f, tn) - env = similar_norm_message_env(tn) - # TODO: replace with `map(f, env)` once `map` is defined on `MessageCache`. - foreach(e -> env[e] = f(env[e]), edges(env)) - return env -end - -""" - identity_norm_message_env(tn) -> MessageCache - -Build a norm-network BP environment with identity-operator messages on every edge — the -"uncorrelated environment" starting point for belief-propagation simple-update gauging -on ⟨tn|tn⟩. - -See also: [`ones_norm_message_env`](@ref), [`randn_norm_message_env`](@ref), -[`rand_norm_message_env`](@ref), [`similar_norm_message_env`](@ref). -""" -identity_norm_message_env(tn) = norm_message_env(one, tn) - -""" - ones_norm_message_env(tn) -> MessageCache - -Build a norm-network BP environment whose per-edge messages have every entry equal to -`1` — the rank-1 outer product of all-ones vectors on each (codomain, domain) pair. - -See also: [`identity_norm_message_env`](@ref), [`randn_norm_message_env`](@ref), -[`rand_norm_message_env`](@ref). -""" -ones_norm_message_env(tn) = norm_message_env(msg -> fill!(msg, one(eltype(msg))), tn) - -randn_norm_message_env(tn) = randn_norm_message_env(Random.default_rng(), tn) - -""" - randn_norm_message_env([rng], tn) -> MessageCache + # Message axis is conj to the tensor it points to. + message = similar_operator(ketview[vertex], braaxis, branames, ketnames) -Build a norm-network BP environment whose per-edge messages have entries drawn from a -standard normal distribution. `rng` defaults to `Random.default_rng()`. + return edge => message + end + end -See also: [`rand_norm_message_env`](@ref), [`identity_norm_message_env`](@ref), -[`ones_norm_message_env`](@ref). -""" -function randn_norm_message_env(rng::Random.AbstractRNG, tn) - return norm_message_env(msg -> randn!(rng, msg), tn) + return messagecache(messages) end -rand_norm_message_env(tn) = rand_norm_message_env(Random.default_rng(), tn) - -""" - rand_norm_message_env([rng], tn) -> MessageCache - -Build a norm-network BP environment whose per-edge messages have entries drawn from a -uniform distribution on `[0, 1)`. `rng` defaults to `Random.default_rng()`. - -See also: [`randn_norm_message_env`](@ref), [`identity_norm_message_env`](@ref), -[`ones_norm_message_env`](@ref). -""" -function rand_norm_message_env(rng::Random.AbstractRNG, tn) - return norm_message_env(msg -> rand!(rng, msg), tn) +function message_environment(f::Base.Callable, nn::NormNetwork) + return map(f, message_environment(undef, nn)) end -# === Double-layer construction and BP wrapper === - -""" - normnetwork(tn) -> norm_tn, linknames_map - -Build the double-layer norm network `⟨tn|tn⟩` together with the per-edge ket→bra name -mapping used to construct it. - -Each ket link axis on every edge is paired with a fresh `uniquename`-generated bra link -name; the bra layer at every vertex is the ket tensor with all of its incident link -names renamed accordingly. The returned `linknames_map` is keyed by both directions of -each undirected edge (the values are shared `Dict`s, so a directed edge and its reverse -look up the same `ketname => braname` table) and is the source of truth for adapting -externally-supplied messages onto the double-layer network. - -Anticipates a future `NormNetwork(tn)` struct that bundles `norm_tn` and `linknames_map` -into a single value with belief-propagation dispatch. -""" -function normnetwork(tn) - linknames_map = Dict( - e => Dict(n => uniquename(n) for n in linknames(tn, e)) - for e in edges(tn) - ) - merge!(linknames_map, Dict(reverse(e) => m for (e, m) in linknames_map)) - norm_tn = ITensorNetwork(underlying_graph(tn)) do v - t = tn[v] - ket_to_bra = Dict(p for e in incident_edges(tn, v) for p in linknames_map[e]) - return t * replacedimnames(n -> get(ket_to_bra, n, n), conj(t)) - end - return norm_tn, linknames_map -end - -""" - beliefpropagation_normnetwork(tn, messages; kwargs...) -> MessageCache - -Run belief propagation on the norm network `⟨tn|tn⟩` (treating `tn` as the ket), -starting from a pre-built operator `MessageCache` `messages` (e.g. from -[`identity_norm_message_env`](@ref) or any of the other `*_norm_message_env` -constructors). - -The norm network built by [`normnetwork`](@ref) is the source of truth for bra-link -names. Each input operator message's codomain (bra) axes are renamed to match the -norm-network's bra names before BP iterates; the converged messages are wrapped back as -operators using those same bra names on output. `kwargs` are forwarded to -`beliefpropagation`. +function beliefpropagation(nn::NormNetwork, messages; kwargs...) + renamed_messages = map(messages) do msg + if !any(name -> has_dimname(KetView(nn), name), dimnames(msg)) + error( + "provided message on does not have have any index \ + names in common with the tensor network contained in the norm." + ) + end -Anticipates a future `beliefpropagation(NormNetwork(tn), messages)` once a `NormNetwork` -wrapper type lands; until then this is the canonical entry point for BP on the norm -network. -""" -function beliefpropagation_normnetwork(tn, messages; kwargs...) - norm_tn, linknames_map = normnetwork(tn) + bramap = Dict(codomainnames(msg) .=> Base.Fix1(namemap, nn).(domainnames(msg))) - # Adapt input messages onto the norm network: rename each operator's codomain - # (bra) axes to the bra names `linknames_map` chose, paired via the operator's - # own domain (ket) → codomain (bra) bijection. - es = collect(keys(messages)) - raws = map(es) do e - msg, ket_to_bra = messages[e], linknames_map[e] - bra_rename = Dict( - cur => ket_to_bra[kn] for - (kn, cur) in zip(domainnames(msg), codomainnames(msg)) - ) - return replacedimnames(n -> get(bra_rename, n, n), state(msg)) + return replacedimnames(name -> get(bramap, name, name), state(msg)) end - raw_messages = Dict(es .=> raws) - cache = beliefpropagation(norm_tn, raw_messages; kwargs...) + cache = _beliefpropagation(nn, renamed_messages; kwargs...) # Re-wrap each converged message as an operator with codomain = bra names and # domain = ket names from the map. - return messagecache(keys(cache)) do e - return operator(cache[e], values(linknames_map[e]), keys(linknames_map[e])) + return messagecache(keys(cache)) do edge + ketnames = linknames(KetView(nn), edge) + branames = linknames(BraView(nn), edge) + return operator(cache[edge], branames, ketnames) end end diff --git a/src/normnetwork.jl b/src/normnetwork.jl new file mode 100644 index 0000000..2f814ae --- /dev/null +++ b/src/normnetwork.jl @@ -0,0 +1,87 @@ +using Dictionaries: Dictionary +using ITensorBase: LazyITensor, lazy, replacedimnames, setname, similar_operator, uniquename +using ITensorNetworksNext + +""" + struct NormNetwork{T, V, I} <: AbstractITensorNetwork{T, V} + +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} + tensornetwork::ITensorNetwork{T, V, I} + namemap::Dictionary{I, I} + function NormNetwork(tn::ITensorNetwork{T, V, I}, map::Dictionary{I, I}) where {T, V, I} + namemap = Dictionary{I, I}() + for (name, vertices) in pairs(tn.dimname_vertices) + if length(vertices) == 2 + insert!(namemap, name, map[name]) + end + end + return new{T, V, I}(tn, namemap) + end +end + +Base.eltype(::Type{<:NormNetwork{T, V, I}}) where {T, V, I} = LazyITensor{I, T} + +function NormNetwork(tn::ITensorNetwork) + return NormNetwork(tn, map(uniquename, keys(tn.dimname_vertices))) +end + +# ====================================== Graphs.jl ======================================= # + +Graphs.edges(nn::NormNetwork) = edges(nn.tensornetwork) +Graphs.vertices(nn::NormNetwork) = vertices(nn.tensornetwork) + +# ==================================== NamedGraphs.jl ==================================== # + +NamedGraphs.vertex_positions(nn::NormNetwork) = vertex_positions(nn.tensornetwork) +NamedGraphs.ordered_vertices(nn::NormNetwork) = ordered_vertices(nn.tensornetwork) +NamedGraphs.position_graph(nn::NormNetwork) = position_graph(nn.tensornetwork) + +# ==================================== DataGraphs.jl ===================================== # + +function DataGraphs.get_vertex_data(nn::NormNetwork, vertex) + A = ket(nn, vertex) + B = conjbra(nn, vertex) + # TODO: implement and use a lazy `conj` via `LazyNamedDimsArrays` here? + return lazy(A) * lazy(conj(B)) +end + +function DataGraphs.is_vertex_assigned(nn::NormNetwork, vertex) + return isassigned(nn.tensornetwork, vertex) +end +# =================================== Dictionaries.jl ==================================== # + +Dictionaries.issettable(::NormNetwork) = false +Dictionaries.isinsertable(::NormNetwork) = false + +# ====================================== interface ======================================= # + +tensornetwork(nn::NormNetwork) = nn.tensornetwork + +function namemap(nn::NormNetwork, name) + if !has_dimname(nn.tensornetwork, name) + error("index name $name not found underlying tensor network.") + end + return get(nn.namemap, name, name) +end + +indmap(nn::NormNetwork, ind) = setname(conj(ind), namemap(nn, name(ind))) + +ket(nn::NormNetwork, vertex) = nn.tensornetwork[vertex] +conjbra(nn::NormNetwork, vertex) = replacedimnames(n -> namemap(nn, n), ket(nn, vertex)) + +bra(nn::NormNetwork, vertex) = conj(conjbra(nn, vertex)) + +""" + normnetwork(tn::ITensorNetwork, [namemap]) -> NormNetwork + +Build the double-layer norm network `⟨tn|tn⟩`, represented lazily as a `NomnNetwork` object. +The optional second argument `namemap` should implement `namemap[ketdimname] = bradimname` for +every link dimension name `ketdimname` in `tn`. If this is not specified, then a name is +generated via the `ITensorBase.uniquename` function. +""" +normnetwork(tn::ITensorNetwork) = NormNetwork(tn) +normnetwork(tn::ITensorNetwork, namemap) = NormNetwork(tn, namemap) diff --git a/src/normnetworkview.jl b/src/normnetworkview.jl new file mode 100644 index 0000000..917d68a --- /dev/null +++ b/src/normnetworkview.jl @@ -0,0 +1,58 @@ +using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned +using Dictionaries: Dictionaries, isinsertable, issettable +using Graphs: Graphs, edges, vertices +using NamedGraphs: NamedGraphs, ordered_vertices, position_graph, vertex_positions + +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) + end +end + +# ==================================== NamedGraphs.jl ==================================== # + +for View in (:KetView, :BraView) + @eval begin + function NamedGraphs.vertex_positions(nnv::$View) + return vertex_positions(nnv.parent) + end + function NamedGraphs.ordered_vertices(nnv::$View) + return ordered_vertices(nnv.parent) + end + + NamedGraphs.position_graph(nnv::$View) = position_graph(nnv.parent) + end +end + +# ==================================== DataGraphs.jl ===================================== # + +DataGraphs.get_vertex_data(nn::KetView, vertex) = ket(nn.parent, vertex) +DataGraphs.get_vertex_data(nn::BraView, vertex) = bra(nn.parent, vertex) + +for View in (:KetView, :BraView) + @eval begin + function DataGraphs.is_vertex_assigned(nnv::$View, vertex) + return isassigned(nnv.parent.tensornetwork, vertex) + end + 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) + end +end diff --git a/src/tensornetwork.jl b/src/tensornetwork.jl index 996075e..4b62b83 100644 --- a/src/tensornetwork.jl +++ b/src/tensornetwork.jl @@ -2,233 +2,218 @@ using Combinatorics: combinations using DataGraphs.DataGraphsPartitionedGraphsExt using DataGraphs: DataGraphs, AbstractDataGraph, DataGraph, edge_data, get_vertices_data, vertex_data, vertex_data_type -using Dictionaries: AbstractDictionary, Indices, dictionary, set!, unset! +using Dictionaries: Dictionaries, AbstractDictionary, Indices, dictionary, set!, unset! using Graphs: AbstractSimpleGraph, rem_edge!, rem_vertex! -using ITensorBase: AbstractITensor, dimnames, lazy +using ITensorBase: + ITensorBase, AbstractITensor, dim, dimnames, dimnametype, name, unnamedtype using NamedGraphs.GraphsExtensions: GraphsExtensions, arrange_edge, arranged_edges, vertextype +using NamedGraphs.OrderedDictionaries: + OrderedDictionary, OrderedIndices, index_positions, ordered_indices using NamedGraphs.PartitionedGraphs: AbstractPartitionedGraph, PartitionedGraphs, QuotientVertex, QuotientVertexVertices, QuotientVertices, departition, partitioned_vertices, partitionedgraph, quotient_graph, quotient_graph_type, quotientvertices -using NamedGraphs: - NamedGraphs, NamedEdge, NamedGraph, Vertices, parent_graph_indices, vertextype - -function _TensorNetwork end - -struct ITensorNetwork{ - V, - VD, - UG <: AbstractGraph{V}, - Tensors <: AbstractDictionary{V, VD}, - } <: - AbstractITensorNetwork{V, VD} - underlying_graph::UG - tensors::Tensors - global @inline function _TensorNetwork( - underlying_graph::UG, tensors::Tensors - ) where {V, VD, UG <: AbstractGraph{V}, Tensors <: AbstractDictionary{V, VD}} - # This assumes the tensor connectivity matches the graph structure. - return new{V, VD, UG, Tensors}(underlying_graph, tensors) +using NamedGraphs: NamedGraphs, NamedEdge, NamedGraph, PositionGraphView, Vertices, + parent_graph_indices, vertextype +using SplitApplyCombine: mapview + +struct ITensorNetwork{T, V, I} <: AbstractITensorNetwork{T, V} + tensors::Dictionary{V, T} + dimname_vertices::Dictionary{I, Set{V}} + underlying_graph::NamedGraph{V} + function ITensorNetwork{T, V, I}(::UndefInitializer, vertices) where {T, V, I} + tensors = Dictionary{V, T}() + dimname_vertices = Dictionary{I, Set{V}}() + underlying_graph = NamedGraph(vertices) + return new{T, V, I}(tensors, dimname_vertices, underlying_graph) end end -# This assumes the tensor connectivity matches the graph structure. -function ITensorNetwork(graph::AbstractGraph, tensors) - return ITensorNetwork(graph, Dictionary(keys(tensors), values(tensors))) -end -function ITensorNetwork(graph::AbstractGraph, tensors::AbstractDictionary) - tn = _TensorNetwork(graph, tensors) - fix_links!(tn) - return tn + +function ITensorNetwork{T}(undef::UndefInitializer, vertices) where {T} + return ITensorNetwork{T, eltype(vertices)}(undef, vertices) end -function ITensorNetwork{V, VD, UG, Tensors}( - graph::UG - ) where {V, VD, UG <: AbstractGraph{V}, Tensors} - return _TensorNetwork(graph, Tensors()) +function ITensorNetwork{T, V}(undef::UndefInitializer, vertices) where {T, V} + return ITensorNetwork{T, V, dimnametype(T)}(undef, vertices) end -function Graphs.rem_vertex!(tn::ITensorNetwork, v) - delete!(tn.tensors, v) - rem_vertex!(tn.underlying_graph, v) +ITensorNetwork(tensors) = ITensorNetwork{valtype(tensors)}(tensors) +ITensorNetwork{T}(tensors) where {T} = ITensorNetwork{T, keytype(tensors)}(tensors) +function ITensorNetwork{T, V}(tensors) where {T, V} + I = dimnametype(T) + tn = ITensorNetwork{T, V, I}(undef, keys(tensors)) + copyto!(tn, tensors) return tn end -# DataGraphs interface - -DataGraphs.underlying_graph(tn::ITensorNetwork) = tn.underlying_graph +ITensorBase.dimnametype(::Type{<:ITensorNetwork{T, V, I}}) where {T, V, I} = I -DataGraphs.is_vertex_assigned(tn::ITensorNetwork, v) = haskey(tn.tensors, v) -DataGraphs.is_edge_assigned(tn::ITensorNetwork, e) = false +Graphs.vertices(tn::ITensorNetwork) = vertices(tn.underlying_graph) -DataGraphs.get_vertex_data(tn::ITensorNetwork, v) = tn.tensors[v] +function NamedGraphs.vertex_positions(graph::ITensorNetwork) + return index_positions(vertices(graph)) +end +function NamedGraphs.ordered_vertices(graph::ITensorNetwork) + return ordered_indices(vertices(graph)) +end -DataGraphs.set_vertex_data!(tn::ITensorNetwork, val, v) = set!(tn.tensors, v, val) +NamedGraphs.position_graph(graph::ITensorNetwork) = position_graph(graph.underlying_graph) -function DataGraphs.underlying_graph_type(type::Type{<:ITensorNetwork}) - return fieldtype(type, :underlying_graph) +function Base.copy(tn::ITensorNetwork{T}) where {T} + tn_dst = ITensorNetwork{T}(undef, vertices(tn)) + copyto!(tn_dst, tn) + return tn_dst end -# For a collection of tensors, return the edges implied by shared indices -# as a list of `edgetype` edges of keys/vertices. -function tensornetwork_edges(edgetype::Type, tensors) - # We need to collect the keys since in the case of `tensors::AbstractDictionary`, - # `keys(tensors)::AbstractIndices`, which is indexed by `keys(tensors)` rather - # than `1:length(keys(tensors))`, which is assumed by `combinations`. - verts = collect(keys(tensors)) - return filter( - !isnothing, map(combinations(verts, 2)) do (v1, v2) - if !isdisjoint(dimnames(tensors[v1]), dimnames(tensors[v2])) - return arrange_edge(edgetype(v1, v2)) - end - return nothing - end - ) -end -tensornetwork_edges(tensors) = tensornetwork_edges(NamedEdge, tensors) +function Graphs.rem_vertex!(tn::ITensorNetwork, vertex) + tensor = tn.tensors[vertex] -function ITensorNetwork(f::Base.Callable, graph::AbstractGraph) - return ITensorNetwork(graph, Dictionary(map(f, vertices(graph)))) -end + for name in dimnames(tensor) + + # If `ind` is associated with an edge, remove the edge. + delete_ind_edge!(tn, name) + + # Delete the vertex from that `ind`s vertex list + # (this index may still be one incident to one other vertex) + vertex_list = tn.dimname_vertices[name] + delete!(vertex_list, vertex) -# Insert trivial links for missing edges, and also check -# the vertices and edges are consistent between the graph and tensors. -function fix_links!(tn::AbstractITensorNetwork) - graph = underlying_graph(tn) - tensors = vertex_data(tn) - @assert issetequal(vertices(graph), keys(tensors)) "Graph vertices and tensor keys must match." - tn_edges = tensornetwork_edges(edgetype(graph), tensors) - tn_edges ⊆ arranged_edges(graph) || - error("The edges in the tensors do not match the graph structure.") - for e in setdiff(arranged_edges(graph), tn_edges) - insertlink!(tn, e) + # If that index is now no longer associated with any vertices, it was dangling, + # and that index should be deleted from the keys of reverse index mapping + isempty(vertex_list) && delete!(tn.dimname_vertices, name) end - return tn -end -# Determine the graph structure from the tensors. -function ITensorNetwork(tensors) - graph = NamedGraph(keys(tensors)) - add_edges!(graph, tensornetwork_edges(tensors)) - return _TensorNetwork(graph, tensors) -end + rem_vertex!(tn.underlying_graph, vertex) + delete!(tn.tensors, vertex) -function Base.copy(tn::ITensorNetwork) - return ITensorNetwork(copy(underlying_graph(tn)), copy(vertex_data(tn))) -end -ITensorNetwork(tn::ITensorNetwork) = copy(tn) -ITensorNetwork{V}(tn::ITensorNetwork{V}) where {V} = copy(tn) -function ITensorNetwork{V}(tn::ITensorNetwork) where {V} - g = convert_vertextype(V, underlying_graph(tn)) - d = dictionary(V(k) => tn[k] for k in vertices(tn)) - return ITensorNetwork(g, d) + return tn end -NamedGraphs.convert_vertextype(::Type{V}, tn::ITensorNetwork{V}) where {V} = tn -NamedGraphs.convert_vertextype(V::Type, tn::ITensorNetwork) = ITensorNetwork{V}(tn) +# Internal (unsafe) +function delete_ind_edge!(tn, ind) + vertex_list = tn.dimname_vertices[ind] -function Graphs.rem_edge!(tn::ITensorNetwork, e) - if !has_edge(underlying_graph(tn), e) - return false - end - if !isempty(linknames(tn, e)) - throw( - ArgumentError( - "cannot remove edge $e due to tensor indices existing on this edge." - ) - ) + if length(vertex_list) == 2 + src, dst = vertex_list + rem_edge!(tn.underlying_graph, src => dst) end - rem_edge!(underlying_graph(tn), e) - return true + + return tn end -function NamedGraphs.similar_graph( - type::Type{<:ITensorNetwork}, - vertices = vertextype(type)[] - ) - DT = fieldtype(type, :tensors) - empty_dict = DT() +# Internal (unsafe) +function delete_ind_vertex!(tn, ind, vertex) + vertex_list = tn.dimname_vertices[ind] - underlying_graph = similar_graph(underlying_graph_type(type), vertices) + delete!(vertex_list, vertex) + isempty(vertex_list) && delete!(tn.dimname_vertices, ind) - return _TensorNetwork(underlying_graph, empty_dict) + return tn end -function NamedGraphs.similar_graph( - graph::ITensorNetwork, - VD::Type, - ::Type{<:Nothing}, - vertices - ) - V = eltype(vertices) - empty_dict = Dictionary{V, VD}() - new_underlying_graph = similar_graph(underlying_graph(graph), vertices) +tensornetwork(f, vertices) = ITensorNetwork(Dict(v => f(v) for v in vertices)) - return _TensorNetwork(new_underlying_graph, empty_dict) -end +Graphs.is_directed(::Type{<:ITensorNetwork}) = false -function NamedGraphs.induced_subgraph_from_vertices(graph::ITensorNetwork, subvertices) - return induced_subgraph_tensornetwork(graph, subvertices) -end +# ====================================== DataGraphs ====================================== # -function induced_subgraph_tensornetwork(graph, subvertices) - underlying_subgraph, vlist = - Graphs.induced_subgraph(underlying_graph(graph), subvertices) +DataGraphs.is_vertex_assigned(tn::ITensorNetwork, vertex) = isassigned(tn.tensors, vertex) +DataGraphs.is_edge_assigned(::ITensorNetwork, _edge) = false - subgraph = ITensorNetwork(underlying_subgraph) do vertex - return graph[vertex] - end +DataGraphs.get_vertex_data(tn::ITensorNetwork, v) = tn.tensors[v] - return subgraph, vlist +function DataGraphs.insert_vertex_data!(tn::ITensorNetwork, vertex, tensor) + add_vertex!(tn.underlying_graph, vertex) + set!_tensornetwork(tn, vertex, tensor) + return tn end -## PartitionedGraphs -function PartitionedGraphs.partitioned_vertices(tn::ITensorNetwork) - return partitioned_vertices(tn.underlying_graph) +function DataGraphs.set_vertex_data!(tn::ITensorNetwork, tensor, vertex) + set!_tensornetwork(tn, vertex, tensor) + return tn end -function PartitionedGraphs.quotient_graph(tn::ITensorNetwork) - ug = quotient_graph(underlying_graph(tn)) +# "upsert" +function set!_tensornetwork(tn::ITensorNetwork, vertex, tensor) + newinds = dimnames(tensor) + + oldinds = get(mapview(dimnames, tn.tensors), vertex, Set()) - inds = Indices(parent_graph_indices(QuotientVertices(tn))) - data = map(v -> tn[QuotientVertex(v)], inds) + # 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) + continue + end + + # Now `ind` must be a new index that's not in `oldinds` - return ITensorNetwork(ug, data) + vertex_list = get!(tn.dimname_vertices, ind, Set()) + if length(vertex_list) > 1 + throw( + ArgumentError( + "index $ind can appear in at most one existing tensor, got $(length(vertex_list))." + ) + ) + end + push!(vertex_list, vertex) + + # Add an edge if the index is now shared between two vertices. + if length(vertex_list) == 2 + src, dst = vertex_list + add_edge!(tn.underlying_graph, src, dst) + end + end + + set!(tn.tensors, vertex, tensor) + + return tn end -# TODO: This method should not be required with a better interface with a better -# DataGraphsPartitionedGraphsExt interface. -function PartitionedGraphs.quotient_graph_type(type::Type{<:ITensorNetwork}) - UG = quotient_graph_type(underlying_graph_type(type)) - VD = Vector{vertex_data_type(type)} - V = vertextype(UG) - return ITensorNetwork{V, VD, UG, Dictionary{V, VD}} + +Dictionaries.isinsertable(::ITensorNetwork) = true + +function DataGraphs.underlying_graph_type(type::Type{<:ITensorNetwork{T, V}}) where {T, V} + return fieldtype(type, :underlying_graph) end -# Partition the underlying graph of the tensor network; does not affect the data. -function PartitionedGraphs.partitionedgraph(tn::ITensorNetwork, parts) - pg = partitionedgraph(underlying_graph(tn), parts) - return ITensorNetwork(pg, copy(vertex_data(tn))) +function Graphs.rem_edge!(::ITensorNetwork, _edge) + return throw( + ErrorException("removing edges from the `ITensorNetwork` type is not supported.") + ) end -PartitionedGraphs.departition(tn::ITensorNetwork) = tn -function PartitionedGraphs.departition( - tn::ITensorNetwork{<:Any, <:Any, <:AbstractPartitionedGraph} +function Graphs.add_edge!(::ITensorNetwork, _edge) + return throw( + ErrorException("Adding edges to the `ITensorNetwork` type is not supported.") ) - return ITensorNetwork(departition(underlying_graph(tn)), vertex_data(tn)) end -NamedGraphs.to_graph_index(::ITensorNetwork, vertex::QuotientVertex) = vertex -# When getting data according the quotient vertices, take a lazy contraction. -function DataGraphs.get_index_data(tn::ITensorNetwork, vertex::QuotientVertex) - data = collect(map(v -> tn[v], vertices(tn, vertex))) - return mapreduce(lazy, *, data) +# PERF: fast lookup compared to `AbstractITensorNetwork` fallback. +dimnamevertices(tn::ITensorNetwork, name) = tn.dimname_vertices[name] + +# PERF: fast lookup compared to `AbstractITensorNetwork` fallback. +has_dimname(tn::ITensorNetwork, name) = haskey(tn.dimname_vertices, name) + +function NamedGraphs.similar_graph( + T::Type{<:ITensorNetwork}, + vertices = vertextype(T)[] + ) + return T(undef, vertices) +end +function NamedGraphs.similar_graph(::ITensorNetwork, VD::Type, vertices) + return ITensorNetwork{VD}(undef, collect(vertices)) end -function DataGraphs.is_graph_index_assigned(tn::ITensorNetwork, vertex::QuotientVertex) - return isassigned(tn, Vertices(vertices(tn, vertex))) + +function NamedGraphs.convert_vertextype(V::Type, tn_src::ITensorNetwork{T}) where {T} + tn_dst = ITensorNetwork{eltype(tn_src), V}(undef, vertices(tn_src)) + copyto!(tn_dst, tn_src) + return tn_dst end -function PartitionedGraphs.quotientview(tn::ITensorNetwork) - qview = QuotientView(underlying_graph(tn)) - tensors = map(qv -> vertex_data(tn)[Indices(qv)], Indices(quotientvertices(tn))) - return ITensorNetwork(qview, tensors) +function NamedGraphs.induced_subgraph_from_vertices(tn::ITensorNetwork, subvertices) + subgraph = similar_graph(tn, subvertices) + copyto!(subgraph, tn, subvertices) + return subgraph, subvertices end diff --git a/test/Project.toml b/test/Project.toml index e4dbe5d..da0665e 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -26,12 +26,12 @@ path = ".." [compat] AlgorithmsInterface = "0.1" Aqua = "0.8.14" -DataGraphs = "0.4" +DataGraphs = "0.5" Dictionaries = "0.4.5" GradedArrays = "0.10" Graphs = "1.13.1" ITensorBase = "0.8" -ITensorNetworksNext = "0.7" +ITensorNetworksNext = "0.8" ITensorPkgSkeleton = "0.3.42" MatrixAlgebraKit = "0.6" NamedGraphs = "0.11.5" diff --git a/test/test_apply_operator.jl b/test/test_apply_operator.jl index 21dd914..5713cca 100644 --- a/test/test_apply_operator.jl +++ b/test/test_apply_operator.jl @@ -1,9 +1,8 @@ using GradedArrays: U1, gradedrange using Graphs: dst, edges, src, vertices -using ITensorBase: ITensorBase as ITB, Index, name, operator, setname, uniquename -using ITensorNetworksNext: ITensorNetwork, apply_operator, apply_operators, - beliefpropagation_normnetwork, identity_norm_message_env, insertlink!, - ones_norm_message_env +using ITensorBase: ITensorBase as ITB, Index, name, named, operator, setname, uniquename +using ITensorNetworksNext: ITensorNetwork, NormNetwork, apply_operator, apply_operators, + beliefpropagation, insertlink!, message_environment, tensornetwork using MatrixAlgebraKit: svd_trunc, truncrank using NamedGraphs.NamedGraphGenerators: named_cycle_graph, named_path_graph using NamedGraphs: NamedGraph @@ -22,13 +21,17 @@ function randn_operator(rng::AbstractRNG, elt::Type, domain_namedaxes) end function random_state(rng::AbstractRNG, elt::Type, g, site_axes; nlayers, trunc) - state = ITensorNetwork(NamedGraph(collect(vertices(g)))) do v + linkinds = Dict(e => Index(1) for e in edges(g)) + + state = tensornetwork(vertices(g)) do v return randn(rng, elt, (site_axes[v],)) end - for e in edges(g) - insertlink!(state, e) + + for edge in edges(g) + insertlink!(state, edge) end - env = identity_norm_message_env(state) + + env = message_environment(one, NormNetwork(state)) for _ in 1:nlayers, e in edges(g) gate = randn_operator(rng, elt, (site_axes[src(e)], site_axes[dst(e)])) state, env = apply_operator(gate, state, env; trunc) @@ -48,10 +51,15 @@ end g = named_cycle_graph(N) site_axes = Dict(v => Index(site_range) for v in vertices(g)) state = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4)) - env = beliefpropagation_normnetwork( - state, ones_norm_message_env(state); + + nn = NormNetwork(state) + + env = beliefpropagation( + nn, + message_environment(msg -> fill!(msg, true), nn); stopping_criterion = (; maxiter = 100, tol = 1.0e-13) ) + for gate in ( randn_operator(rng, T, (site_axes[2],)), randn_operator(rng, T, (site_axes[2], site_axes[3])), @@ -66,15 +74,21 @@ end g = named_path_graph(N) site_axes = Dict(v => Index(site_range) for v in vertices(g)) state = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4)) - env = beliefpropagation_normnetwork( - state, ones_norm_message_env(state); + + nn = NormNetwork(state) + + env = beliefpropagation( + nn, + message_environment(msg -> fill!(msg, true), nn); stopping_criterion = (; maxiter = 100, tol = 1.0e-13) ) + gate = randn_operator(rng, T, (site_axes[2], site_axes[3])) gated_full = ITB.apply(gate, prod(state)) left = [name(site_axes[v]) for v in 1:2] U, S, Vt = svd_trunc(gated_full, left; trunc = truncrank(k)) gated, _ = apply_operator(gate, state, env; trunc = truncrank(k)) + @test prod(gated) ≈ U * S * Vt rtol = eps(real(T))^(1 / 3) end @@ -83,8 +97,12 @@ end g = named_cycle_graph(N) site_axes = Dict(v => Index(site_range) for v in vertices(g)) state = random_state(rng, T, g, site_axes; nlayers = 2, trunc = truncrank(4)) - env = beliefpropagation_normnetwork( - state, ones_norm_message_env(state); + + nn = NormNetwork(state) + + env = beliefpropagation( + nn, + message_environment(msg -> fill!(msg, true), nn); stopping_criterion = (; maxiter = 100, tol = 1.0e-13) ) g1 = randn_operator(rng, T, (site_axes[2], site_axes[3])) diff --git a/test/test_basics.jl b/test/test_basics.jl index 64636f0..1edfe6c 100644 --- a/test/test_basics.jl +++ b/test/test_basics.jl @@ -1,7 +1,7 @@ using Dictionaries: Indices using Graphs: dst, edges, has_edge, ne, nv, src, vertices using ITensorBase: Index, dimnames -using ITensorNetworksNext: ITensorNetwork, linkinds, siteinds +using ITensorNetworksNext: ITensorNetwork, linkinds, siteinds, tensornetwork using NamedGraphs.GraphsExtensions: arranged_edges, incident_edges using NamedGraphs.NamedGraphGenerators: named_grid using Test: @test, @testset @@ -11,25 +11,14 @@ using Test: @test, @testset dims = (3, 3) g = named_grid(dims) s = Dict(v => Index(2) for v in vertices(g)) - tn = ITensorNetwork(g) do v + tn = tensornetwork(vertices(g)) do v return randn(s[v]) end @test nv(tn) == 9 - @test ne(tn) == ne(g) + @test ne(tn) == 0 # zero link indices @test issetequal(vertices(tn), vertices(g)) - @test issetequal(arranged_edges(tn), arranged_edges(g)) for v in vertices(tn) - @test siteinds(tn, v) == [s[v]] - end - for v1 in vertices(tn) - for v2 in vertices(tn) - v1 == v2 && continue - haslink = !isempty(linkinds(tn, v1 => v2)) - @test haslink == has_edge(tn, v1 => v2) - end - end - for e in edges(tn) - @test isone(length(only(linkinds(tn, e)))) + @test siteinds(tn, v) == (s[v],) end end @testset "Construct ITensorNetwork partition function" begin @@ -37,7 +26,7 @@ using Test: @test, @testset g = named_grid(dims) l = Dict(e => Index(2) for e in edges(g)) l = merge(l, Dict(reverse(e) => l[e] for e in edges(g))) - tn = ITensorNetwork(g) do v + tn = tensornetwork(vertices(g)) do v is = map(e -> l[e], incident_edges(g, v)) return randn(Tuple(is)) end @@ -51,8 +40,8 @@ using Test: @test, @testset for v1 in vertices(tn) for v2 in vertices(tn) v1 == v2 && continue - haslink = !isempty(linkinds(tn, v1 => v2)) - @test haslink == has_edge(tn, v1 => v2) + # haslink = !isempty(linkinds(tn, v1 => v2)) + # @test haslink == has_edge(tn, v1 => v2) end end for e in edges(tn) diff --git a/test/test_beliefpropagation.jl b/test/test_beliefpropagation.jl index 8b548d0..3ecf8ee 100644 --- a/test/test_beliefpropagation.jl +++ b/test/test_beliefpropagation.jl @@ -1,11 +1,11 @@ -using AlgorithmsInterface: AlgorithmsInterface as AI -using DataGraphs: edge_data, edge_data_type +import AlgorithmsInterface as AI +using DataGraphs: DataGraphs, DataGraph, edge_data, edge_data_type using Dictionaries: Dictionary, dictionary, set! using Graphs: AbstractGraph, dst, edges, has_edge, src, vertices using ITensorBase: ITensor, Index, inds, name, noprime, prime using ITensorNetworksNext: ITensorNetworksNext, ITensorNetwork, MessageCache, StopWhenConverged, bethe_free_energy, edge_scalar, incoming_messages, linkinds, - messagecache, region_scalar, subgraph, vertex_scalar, vertex_scalars + messagecache, region_scalar, subgraph, tensornetwork, vertex_scalar, vertex_scalars using LinearAlgebra: LinearAlgebra using NamedGraphs.GraphsExtensions: all_edges, arranged_edges, incident_edges, vertextype using NamedGraphs.NamedGraphGenerators: named_comb_tree, named_grid, named_path_graph @@ -14,12 +14,10 @@ using StableRNGs: StableRNG using Test: @test, @testset function spin_ice_tensornetwork(g) - links = Dictionary( - edges(g), - [Index(2) for e in edges(g)] - # [Index(2; tags = "edge " => "e$(src(e))_$(dst(e))") for e in edges(g)] - ) - links = merge(links, Dictionary(reverse.(edges(g)), [links[e] for e in edges(g)])) + links = DataGraph(g) + for e in edges(g) + links[e] = Index(2) + end ts = Dictionary{vertextype(g), ITensor}() for v in vertices(g) @@ -34,7 +32,7 @@ function spin_ice_tensornetwork(g) t = t_data[linkinds...] set!(ts, v, t) end - return ITensorNetwork(g, ts) + return ITensorNetwork(ts) end @testset "Belief propagation" begin @@ -42,10 +40,11 @@ end @testset "Basics" begin dims = (3, 3) g = named_grid(dims) + l = Dict(e => Index(2) for e in edges(g)) l = merge(l, Dict(reverse(e) => l[e] for e in edges(g))) - tn = ITensorNetwork(g) do v + tn = tensornetwork(vertices(g)) do v is = map(e -> l[e], incident_edges(g, v)) return randn(Tuple(is)) end @@ -85,7 +84,8 @@ end 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))) - tn = ITensorNetwork(g) do v + + tn = tensornetwork(vertices(g)) do v is = map(e -> l[e], incident_edges(g, v)) return randn(ComplexF32, Tuple(is)) end @@ -118,7 +118,8 @@ end g = named_grid((3,)) l = Dict(e => Index(2) for e in edges(g)) l = merge(l, Dict(reverse(e) => l[e] for e in edges(g))) - tn = ITensorNetwork(g) do v + + tn = tensornetwork(vertices(g)) do v is = map(e -> l[e], incident_edges(g, v)) return randn(Tuple(is)) end @@ -134,7 +135,8 @@ end g = named_grid((2,)) l = Dict(e => Index(2) for e in edges(g)) l = merge(l, Dict(reverse(e) => l[e] for e in edges(g))) - tn = ITensorNetwork(g) do v + + tn = tensornetwork(vertices(g)) do v is = map(e -> l[e], incident_edges(g, v)) return randn(Tuple(is)) end @@ -154,14 +156,16 @@ end #Chain of tensors dims = (2, 1) - g = named_grid(dims) - l = Dict(e => Index(2) for e in edges(g)) - l = merge(l, Dict(reverse(e) => l[e] for e in edges(g))) + g = DataGraph(named_grid(dims)) # graph to hold the links. + for edge in edges(g) + g[edge] = Index(2) + end - tn = ITensorNetwork(g) do v - is = map(e -> l[e], incident_edges(g, v)) - return randn(rng, T, Tuple(is)) + tensors = map(vertices(g)) do vertex + is = map(edge -> g[edge], incident_edges(g, vertex)) + return randn(T, Tuple(is)) end + tn = ITensorNetwork(tensors) messages = Dict( edge => ones(T, Tuple(linkinds(tn, edge))) for edge in all_edges(g) @@ -176,13 +180,15 @@ end #Tree of tensors dims = (4, 3) - g = named_comb_tree(dims) - l = Dict(e => Index(3) for e in edges(g)) - l = merge(l, Dict(reverse(e) => l[e] for e in edges(g))) - tn = ITensorNetwork(g) do v - is = map(e -> l[e], incident_edges(g, v)) - return randn(rng, T, Tuple(is)) + g = DataGraph(named_comb_tree(dims)) # graph to hold the links. + for edge in edges(g) + g[edge] = Index(3) + end + tensors = map(vertices(g)) do vertex + is = map(edge -> g[edge], incident_edges(g, vertex)) + return randn(T, Tuple(is)) end + tn = ITensorNetwork(tensors) messages = Dict( edge => ones(T, Tuple(linkinds(tn, edge))) for edge in all_edges(g) diff --git a/test/test_contract_network.jl b/test/test_contract_network.jl index b5383ef..a6cdb30 100644 --- a/test/test_contract_network.jl +++ b/test/test_contract_network.jl @@ -1,7 +1,7 @@ -using Graphs: edges +using Graphs: edges, vertices using ITensorBase: Greedy, Index, Optimal -using ITensorNetworksNext: - Exact, ITensorNetwork, LeftAssociative, contract_network, linkinds, siteinds +using ITensorNetworksNext: Exact, ITensorNetwork, LeftAssociative, contract_network, + linkinds, siteinds, tensornetwork using NamedGraphs.GraphsExtensions: arranged_edges, incident_edges using NamedGraphs.NamedGraphGenerators: named_grid using TensorOperations: TensorOperations @@ -28,7 +28,7 @@ using Test: @test, @testset g = named_grid(dims) l = Dict(e => Index(2) for e in edges(g)) l = merge(l, Dict(reverse(e) => l[e] for e in edges(g))) - tn = ITensorNetwork(g) do v + tn = tensornetwork(vertices(g)) do v is = map(e -> l[e], incident_edges(g, v)) return randn(Tuple(is)) end @@ -39,5 +39,8 @@ using Test: @test, @testset @test abs(z1 - z2) / abs(z1) <= 1.0e3 * eps(Float64) @test abs(z1 - z3) / abs(z1) <= 1.0e3 * eps(Float64) + + @test z1 ≈ z2 + @test z1 ≈ z3 end end diff --git a/test/test_itensornetworkgenerators.jl b/test/test_itensornetworkgenerators.jl index 3656ca6..c4b1da2 100644 --- a/test/test_itensornetworkgenerators.jl +++ b/test/test_itensornetworkgenerators.jl @@ -1,3 +1,4 @@ +using DataGraphs: DataGraph using Graphs: edges, ne, nv, vertices using ITensorBase: Index, inds using ITensorNetworksNext.ITensorNetworkGenerators: delta, delta_network, ising_network @@ -29,8 +30,13 @@ using Test: @test, @testset dims = (4,) β = 0.4 g = named_grid(dims; periodic) - ldict = Dict(e => Index(2) for e in edges(g)) - l(e) = get(() -> ldict[reverse(e)], ldict, e) + # ldict = Dict(e => Index(2) for e in edges(g)) + # l(e) = get(() -> ldict[reverse(e)], ldict, e) + ldict = DataGraph(g) + for e in edges(g) + ldict[e] = Index(2) + end + l = e -> ldict[e] tn = ising_network(l, β, g) @test nv(tn) == 4 @test ne(tn) == ne(g) diff --git a/test/test_normnetwork.jl b/test/test_normnetwork.jl new file mode 100644 index 0000000..86b088e --- /dev/null +++ b/test/test_normnetwork.jl @@ -0,0 +1,170 @@ +using Base.Broadcast: materialize +using DataGraphs: is_vertex_assigned +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, bra, conjbra, + indmap, ket, namemap, normnetwork, tensornetwork +using LinearAlgebra: norm +using NamedGraphs.GraphsExtensions: incident_edges +using NamedGraphs.NamedGraphGenerators: named_grid, named_path_graph +using NamedGraphs: NamedEdge +using Test: @test, @test_throws, @testset + +# Contract a (possibly double-layer) network into a single tensor by multiplying all +# of its vertex tensors together. For a `NormNetwork` the vertex data are lazy products +# `ket * conj(bra)`, so the result is a lazy expression that we materialize. +contract(tn) = materialize(prod(tn)) + +# 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 + +@testset "`NormNetwork`" begin + @testset "Basics" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + nn = NormNetwork(tn) + + # `normnetwork` is the public constructor and agrees with `NormNetwork`. + @test normnetwork(tn) isa NormNetwork + @test nn isa NormNetwork + + # The underlying tensor network is accessible and is the same object. + @test tensornetwork(nn) === tn + + # The norm network shares the graph structure of the underlying network. + @test issetequal(vertices(nn), vertices(tn)) + @test issetequal(edges(nn), edges(tn)) + + # `eltype` is the type of the (lazy double-layer) vertex data. + @test eltype(nn) === typeof(nn[first(vertices(nn))]) + + # Vertex data is assigned wherever the underlying network is. + @test is_vertex_assigned(nn, 1) + + # The norm network is neither settable nor insertable (it is a lazy view). + @test !issettable(nn) + @test !isinsertable(nn) + end + + @testset "ket / bra / conjbra and the name map" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + nn = NormNetwork(tn) + + # `ket` returns the underlying tensor untouched. + @test ket(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 namemap(nn, sname) == sname + @test sname in name.(inds(ket(nn, 2))) + @test sname in name.(inds(conjbra(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 namemap(nn, lname) != lname + @test lname in name.(inds(ket(nn, 2))) + @test !(lname in name.(inds(conjbra(nn, 2)))) + @test namemap(nn, lname) in name.(inds(conjbra(nn, 2))) + + # `bra` is the elementwise conjugate of `conjbra` and carries the same indices. + @test inds(bra(nn, 2)) == inds(conjbra(nn, 2)) + + # `indmap` conjugates an index and renames it according to the name map. + ind = only(i for i in inds(ket(nn, 2)) if name(i) == lname) + @test name(indmap(nn, ind)) == namemap(nn, name(ind)) + @test indmap(nn, ind) == setname(conj(ind), namemap(nn, name(ind))) + + # Querying the name map with an index name absent from the network errors. + @test_throws ErrorException namemap(nn, name(Index(2))) + end + + @testset "custom name map" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + + # A user-supplied map dictates the bra-layer name for each link. + custom = map(uniquename, keys(tn.dimname_vertices)) + nn = normnetwork(tn, custom) + + lname = name(l[NamedEdge(1 => 2)]) + @test namemap(nn, lname) == custom[lname] + @test namemap(nn, lname) in name.(inds(conjbra(nn, 2))) + end + + @testset "`KetView` / `BraView`" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + nn = NormNetwork(tn) + + kv = KetView(nn) + bv = BraView(nn) + + # Views share the graph structure of the underlying network. + @test issetequal(vertices(kv), vertices(tn)) + @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] === ket(nn, v) + @test inds(bv[v]) == inds(bra(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) + @test !issettable(bv) + @test !isinsertable(bv) + end + + @testset "contraction / physics" begin + @testset "single normalized tensor contracts to 1" begin + s = Index(3) + v = randn(s) + v = v / norm(v) + tn = ITensorNetwork(Dict(1 => v)) + nn = NormNetwork(tn) + + # ⟨ψ|ψ⟩ for a single normalized site tensor is 1. + @test contract(nn)[] ≈ 1 + end + + @testset "$T" for T in (Float64, ComplexF64) + g = named_grid((2, 2)) + tn, l, s = random_state(T, g) + + # The norm network contracts to ⟨tn|tn⟩ = ‖prod(tn)‖², a real nonnegative number. + z = contract(NormNetwork(tn))[] + @test z ≈ norm(prod(tn))^2 + @test imag(z) ≈ 0 atol = 1.0e-12 * abs(z) + @test real(z) > 0 + + # Rescaling a single tensor by 1/√z normalizes the state, so ⟨tn|tn⟩ = 1. + tn[first(vertices(tn))] = tn[first(vertices(tn))] / sqrt(real(z)) + @test contract(NormNetwork(tn))[] ≈ 1 + + # The contracted norm does not depend on the chosen bra-layer name map. + custom = map(uniquename, keys(tn.dimname_vertices)) + @test contract(normnetwork(tn, custom))[] ≈ contract(NormNetwork(tn))[] + end + end +end diff --git a/test/test_tensornetwork.jl b/test/test_tensornetwork.jl index 4efe4cb..8d044a5 100644 --- a/test/test_tensornetwork.jl +++ b/test/test_tensornetwork.jl @@ -1,9 +1,10 @@ -using DataGraphs: assigned_edge_data, assigned_vertex_data, underlying_graph, vertex_data +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_vertex!, src, vertices -using ITensorBase: Index, LazyITensor -using ITensorNetworksNext: - ITensorNetwork, fix_edges!, linkaxes, linkinds, linknames, siteaxes, siteinds, sitenames + is_directed, ne, nv, rem_edge!, rem_vertex!, src, vertices +using ITensorBase: Index, LazyITensor, inds +using ITensorNetworksNext: ITensorNetwork, has_ind, linkaxes, linkinds, linknames, siteaxes, + siteinds, sitenames, tensornetwork using NamedGraphs.GraphsExtensions: incident_edges, subgraph, vertextype using NamedGraphs.NamedGraphGenerators: named_grid, named_path_graph using NamedGraphs.PartitionedGraphs: AbstractPartitionedGraph, QuotientVertex, departition, @@ -15,9 +16,8 @@ using Test: @test, @test_throws, @testset @testset "`ITensorNetwork`" begin @testset "Basics" begin g = named_grid((2, 2)) - s = Dict(v => Index(2) for v in vertices(g)) - tn = ITensorNetwork(g) do v - return randn(s[v]) + tn = tensornetwork(vertices(g)) do _ + return randn(Index(2)) end # `iterate` works (delegates to `vertex_data`). @@ -36,23 +36,65 @@ using Test: @test, @test_throws, @testset s_default = sprint(show, tn) @test occursin("vertices", s_default) + # No link indices so should have no edges + @test ne(tn) == 0 + + j = Index(2) + tn[1, 1] = randn(j) + tn[2, 1] = randn(j) + + # No link indices so should have no edges + @test ne(tn) == 1 + @test has_edge(tn, (1, 1) => (2, 1)) + # `setindex!` for edges is intentionally unimplemented. e = first(edges(tn)) - @test_throws ErrorException tn[e] = randn(2, 2) - @test_throws ErrorException tn[src(e) => dst(e)] = randn(2, 2) + @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)) + + tn[1, 1] = randn(Index(2)) + tn[2, 1] = randn(Index(2)) + + @test !has_edge(tn, (1, 1) => (2, 1)) + @test ne(tn) == 0 rem_vertex!(tn, (2, 2)) @test !has_vertex(tn, (2, 2)) - add_vertex!(tn, (2, 2)) + insert!(tn, (2, 2), randn(Index(2))) @test has_vertex(tn, (2, 2)) - @test !isassigned(tn, (2, 2)) - - # Test `fix_edges!` removes edges where there is no link index - t = randn(s[(2, 2)]) - tn[(2, 2)] = t - add_edge!(tn.underlying_graph, (1, 2) => (2, 2)) - fix_edges!(tn, (2, 2)) - @test !has_edge(tn, (1, 2) => (2, 2)) + @test isassigned(tn, (2, 2)) + + links = DataGraph(named_path_graph(4)) + i = Index(2) + j = Index(3) + k = Index(4) + links[1 => 2] = i + links[2 => 3] = j + links[3 => 4] = k + + tn = tensornetwork(vertices(links)) do v + indices = map(e -> getindex(links, e), incident_edges(links, v)) + return randn(Tuple(indices)) + end + + @test has_ind(tn, i) + @test has_ind(tn, j) + @test has_ind(tn, k) + + ip = Index(2) + + tn[1] = tn[1] * randn((i, ip)) + @test has_ind(tn, ip) + @test has_ind(tn, j) + @test has_ind(tn, k) + + @test inds(tn[1]) == (ip,) + @test inds(tn[2]) == (i, j) + @test inds(tn[3]) == (j, k) + @test inds(tn[4]) == (k,) end @testset "link and site functions" begin @@ -60,7 +102,7 @@ using Test: @test, @test_throws, @testset 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 = ITensorNetwork(g) do v + tn = tensornetwork(vertices(g)) do v is = map(e -> l[e], incident_edges(g, v)) return randn((s[v], is...)) end @@ -82,9 +124,11 @@ using Test: @test, @test_throws, @testset @testset "`subgraph`" begin g = named_grid((3,)) + l = Dict(e => Index(2) for e in edges(g)) l = merge(l, Dict(reverse(e) => l[e] for e in edges(g))) - tn = ITensorNetwork(g) do v + + tn = tensornetwork(vertices(g)) do v is = map(e -> l[e], incident_edges(g, v)) return randn(Tuple(is)) end @@ -100,7 +144,7 @@ using Test: @test, @test_throws, @testset dims = (3, 3) g = named_grid(dims) s = Dict(v => Index(2) for v in vertices(g)) - tn = ITensorNetwork(g) do v + tn = tensornetwork(vertices(g)) do v return randn(s[v]) end @@ -131,105 +175,4 @@ using Test: @test, @test_throws, @testset @test vertextype(ctn) == Tuple{Float64, Float64} @test collect(vertex_data(ctn)) == collect(vertex_data(tn)) end - - @testset "`PartitionedGraphs`" begin - dims = (3, 3) - g = named_grid(dims) - s = Dict(v => Index(2) for v in vertices(g)) - tn = ITensorNetwork(g) do v - return randn(s[v]) - end - - # Row partition: each partition is one row of the grid. - row_parts = [[(i, j) for i in 1:dims[1]] for j in 1:dims[2]] - - @testset "default `partitioned_vertices`" begin - # By default the entire underlying graph is one partition. - pvs = partitioned_vertices(tn) - @test length(pvs) == 1 - @test issetequal(only(pvs), vertices(tn)) - end - - @testset "default `quotientvertices`" begin - qvs = collect(quotientvertices(tn)) - @test length(qvs) == 1 - @test only(qvs) isa QuotientVertex - end - - @testset "`tn[QuotientVertex(...)]` (default)" begin - qv = only(collect(quotientvertices(tn))) - data = tn[qv] - @test data isa LazyITensor - end - - @testset "`quotient_graph` (default partitioning)" begin - qtn = quotient_graph(tn) - @test qtn isa ITensorNetwork - @test nv(qtn) == 1 - @test ne(qtn) == 0 - v = only(collect(vertices(qtn))) - @test qtn[v] isa LazyITensor - end - - @testset "`quotient_graph_type`" begin - QT = quotient_graph_type(typeof(tn)) - @test QT <: ITensorNetwork - qtn = quotient_graph(tn) - @test vertextype(qtn) === vertextype(QT) - end - - @testset "`partitionedgraph(tn, parts)`" begin - ptn = partitionedgraph(tn, row_parts) - @test ptn isa ITensorNetwork - # The set of underlying vertices/edges is preserved. - @test issetequal(vertices(ptn), vertices(tn)) - @test issetequal(edges(ptn), edges(tn)) - @test nv(ptn) == nv(tn) - @test ne(ptn) == ne(tn) - # Vertex data is copied, not aliased. - @test collect(vertex_data(ptn)) == collect(vertex_data(tn)) - @test vertex_data(ptn) !== vertex_data(tn) - end - - @testset "`partitioned_vertices` of partitioned tn" begin - ptn = partitionedgraph(tn, row_parts) - pvs = partitioned_vertices(ptn) - @test length(pvs) == dims[2] - for part in pvs - @test length(part) == dims[1] - end - @test issetequal(reduce(vcat, pvs), vertices(tn)) - end - - @testset "`tn[QuotientVertex(...)]` (partitioned)" begin - ptn = partitionedgraph(tn, row_parts) - for qv in quotientvertices(ptn) - @test ptn[qv] isa LazyITensor - end - end - - @testset "`quotient_graph` of partitioned tn" begin - ptn = partitionedgraph(tn, row_parts) - qtn = quotient_graph(ptn) - @test qtn isa ITensorNetwork - @test nv(qtn) == dims[2] - # The row-partitioned grid quotients to a path graph of length `dims[2]`. - @test ne(qtn) == dims[2] - 1 - for v in vertices(qtn) - @test qtn[v] isa LazyITensor - end - end - - @testset "`departition`" begin - # `departition` on a non-partitioned tn returns itself. - @test departition(tn) === tn - - # `departition` on a partitioned tn unwraps one layer of partitioning. - ptn = partitionedgraph(tn, row_parts) - dtn = departition(ptn) - @test dtn isa ITensorNetwork - @test issetequal(vertices(dtn), vertices(tn)) - @test issetequal(edges(dtn), edges(tn)) - end - end end