diff --git a/Project.toml b/Project.toml index c614f208..df9a94f6 100644 --- a/Project.toml +++ b/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 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" diff --git a/src/ITensorNetworkGenerators/ising_network.jl b/src/ITensorNetworkGenerators/ising_network.jl index 702b8787..e2ebe8b3 100644 --- a/src/ITensorNetworkGenerators/ising_network.jl +++ b/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) diff --git a/src/ITensorNetworksNext.jl b/src/ITensorNetworksNext.jl index 0298a85f..838d7339 100644 --- a/src/ITensorNetworksNext.jl +++ b/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") diff --git a/src/abstracttensornetwork.jl b/src/abstracttensornetwork.jl index eabe39ba..197e0c87 100644 --- a/src/abstracttensornetwork.jl +++ b/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,13 +37,16 @@ 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) @@ -51,43 +54,46 @@ function linkaxes(tn::AbstractGraph, edge::Pair) 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 diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl index fa387253..fbb8f625 100644 --- a/src/beliefpropagation/beliefpropagation.jl +++ b/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) diff --git a/src/beliefpropagation/messagecache.jl b/src/beliefpropagation/messagecache.jl index ec91cd58..ebbfd759 100644 --- a/src/beliefpropagation/messagecache.jl +++ b/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 diff --git a/src/bilinearforms/abstractbilinearformnetwork.jl b/src/bilinearforms/abstractbilinearformnetwork.jl new file mode 100644 index 00000000..0d4441b0 --- /dev/null +++ b/src/bilinearforms/abstractbilinearformnetwork.jl @@ -0,0 +1,189 @@ +using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned +using Dictionaries: Dictionaries, isinsertable, issettable +using Graphs: Graphs, edges, vertices +using ITensorBase: ITensorBase, conj, inds, nametype, rename +using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex + +""" + abstract type AbstractBilinearFormNetwork{T, V, I} <: AbstractITensorNetwork{T, V} + +Supertype of the lazy multi-layer networks built from a ket layer of type +`ITensorNetwork{T, V, I}` and a ket→bra index name mapping. + +A subtype implements [`ketnetwork`](@ref) and [`branamemap`](@ref), from which the graph +structure, [`braname`](@ref), [`kettensor`](@ref) and [`bratensor`](@ref) are defined, and, +where it has an operator layer, [`operatornetwork`](@ref) and [`operatortensor`](@ref). The layers as whole networks are returned by [`ketnetwork`](@ref), +[`branetwork`](@ref) and [`operatornetwork`](@ref). +""" +abstract type AbstractBilinearFormNetwork{T, V, I} <: AbstractITensorNetwork{T, V} end + +# The checks on a ket→bra name map shared by the bilinear-form constructors; `layers` are the +# networks other than `ket` whose index names a bra name must not reuse. +function check_branamemap(ket, braname, layers...) + if !issubset(linknames(ket), keys(braname)) + error("every link name of the ket layer must have a bra name.") + end + if !allunique(braname) + error("the bra names of a bilinear-form network must be distinct.") + end + for name in braname + if has_dimname(ket, name) || any(layer -> has_dimname(layer, name), layers) + error("bra name $name is already an index name of another layer.") + end + end + return nothing +end + +# ====================================== Graphs.jl ======================================= # + +Graphs.edges(bn::AbstractBilinearFormNetwork) = edges(ketnetwork(bn)) +Graphs.vertices(bn::AbstractBilinearFormNetwork) = vertices(ketnetwork(bn)) + +# ==================================== NamedGraphs.jl ==================================== # + +function NamedGraphs.encoded_vertex(bn::AbstractBilinearFormNetwork, vertex) + return encoded_vertex(ketnetwork(bn), vertex) +end +function NamedGraphs.decoded_vertex(bn::AbstractBilinearFormNetwork, code::Integer) + return decoded_vertex(ketnetwork(bn), code) +end +NamedGraphs.encoded_graph(bn::AbstractBilinearFormNetwork) = encoded_graph(ketnetwork(bn)) + +# ==================================== DataGraphs.jl ===================================== # + +function DataGraphs.is_vertex_assigned(bn::AbstractBilinearFormNetwork, vertex) + return isassigned(ketnetwork(bn), vertex) +end + +# =================================== Dictionaries.jl ==================================== # + +Dictionaries.issettable(::AbstractBilinearFormNetwork) = false +Dictionaries.isinsertable(::AbstractBilinearFormNetwork) = false + +# ====================================== interface ======================================= # + +# The name type of the indices, the same as the ket network's. +function ITensorBase.nametype( + ::Type{<:AbstractBilinearFormNetwork{T, V, I}} + ) where {T, V, I} + return I +end +ITensorBase.nametype(bn::AbstractBilinearFormNetwork) = nametype(typeof(bn)) + +""" + braname(bn::AbstractBilinearFormNetwork, name) + +The bra-layer index name corresponding to the ket-layer index name `name`. +""" +function braname(bn::AbstractBilinearFormNetwork, name) + if !has_dimname(ketnetwork(bn), name) + error("index name $name not found underlying tensor network.") + end + # A name absent from the map has no separate bra copy and maps to itself: a site index of a + # norm network, or a site index a quadratic form's operator does not act on. + return get(branamemap(bn), name, name) +end + +""" + branamemap(bn::AbstractBilinearFormNetwork) + +The ket→bra name map, holding a bra name for each ket index name that has a separate bra copy. +""" +function branamemap end + +""" + kettensor(bn::AbstractBilinearFormNetwork, vertex) + +The ket-layer tensor at `vertex`. +""" +kettensor(bn::AbstractBilinearFormNetwork, vertex) = ketnetwork(bn)[vertex] + +""" + operatortensor(bn::AbstractBilinearFormNetwork, vertex) + +The operator-layer tensor at `vertex`, with its index names renamed so that its input legs +meet the ket layer and its output legs meet the bra layer. +""" +function operatortensor end + +function conj_bratensor(bn::AbstractBilinearFormNetwork, vertex) + return rename(n -> braname(bn, n), kettensor(bn, vertex)) +end + +""" + bratensor(bn::AbstractBilinearFormNetwork, vertex) + +The bra-layer tensor at `vertex`. +""" +bratensor(bn::AbstractBilinearFormNetwork, vertex) = conj(conj_bratensor(bn, vertex)) + +""" + ketnetwork(bn::AbstractBilinearFormNetwork) + +The ket-layer network of `bn`. +""" +function ketnetwork end + +""" + operatornetwork(bn::AbstractBilinearFormNetwork) + +The operator-layer network of `bn`, for a subtype that has an operator layer. +""" +function operatornetwork end + +""" + branetwork(bn::AbstractBilinearFormNetwork) + +The bra-layer network of `bn`. Unless a subtype stores its bra layer as a network, this is a +`BraView`, whose tensors are built by [`bratensor`](@ref) when accessed. +""" +branetwork(bn::AbstractBilinearFormNetwork) = BraView(bn) + +""" + struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <: AbstractITensorNetwork{T, V} + +The bra layer of the bilinear-form network `parent(view)`, with each vertex tensor built by +[`bratensor`](@ref) when accessed. Its graph structure and mutability are those of the parent. +""" +struct BraView{T, V, I, P <: AbstractBilinearFormNetwork{T, V, I}} <: + AbstractITensorNetwork{T, V} + parent::P + function BraView(parent::AbstractBilinearFormNetwork{T, V, I}) where {T, V, I} + return new{T, V, I, typeof(parent)}(parent) + end +end + +Base.parent(nnv::BraView) = nnv.parent + +# ==================================== DataGraphs.jl ===================================== # + +DataGraphs.get_vertex_data(nnv::BraView, vertex) = bratensor(parent(nnv), vertex) +function DataGraphs.is_vertex_assigned(nnv::BraView, vertex) + return is_vertex_assigned(parent(nnv), vertex) +end + +# ==================================== ITensorBase.jl ==================================== # + +function ITensorBase.inds(nnv::BraView, vertex) + return conj.(inds(conj_bratensor(parent(nnv), vertex))) +end + +# ====================================== Graphs.jl ======================================= # + +Graphs.edges(nnv::BraView) = edges(parent(nnv)) +Graphs.vertices(nnv::BraView) = vertices(parent(nnv)) + +# ==================================== NamedGraphs.jl ==================================== # + +function NamedGraphs.encoded_vertex(nnv::BraView, vertex) + return encoded_vertex(parent(nnv), vertex) +end +function NamedGraphs.decoded_vertex(nnv::BraView, code::Integer) + return decoded_vertex(parent(nnv), code) +end +NamedGraphs.encoded_graph(nnv::BraView) = encoded_graph(parent(nnv)) + +# =================================== Dictionaries.jl ==================================== # + +Dictionaries.issettable(nnv::BraView) = issettable(parent(nnv)) +Dictionaries.isinsertable(nnv::BraView) = isinsertable(parent(nnv)) diff --git a/src/bilinearforms/normnetwork.jl b/src/bilinearforms/normnetwork.jl new file mode 100644 index 00000000..8410ae49 --- /dev/null +++ b/src/bilinearforms/normnetwork.jl @@ -0,0 +1,85 @@ +using Dictionaries: Dictionary, Indices, dictionary +using ITensorBase: similar_operator +using ITensorNetworksNext + +""" + struct NormNetwork{T, V, I} <: AbstractBilinearFormNetwork{T, V, I} + +Lazy wrapper representing the norm `⟨tn|tn⟩` of `tn::ITensorNetwork{T, V, I}`, +together with a per-edge ket→bra name mapping that, for each index in the ket layer, defines +the name of the corresponding index in the bra layer. +""" +struct NormNetwork{T, V, I} <: AbstractBilinearFormNetwork{T, V, I} + ket::ITensorNetwork{T, V, I} + braname::Dictionary{I, I} + function NormNetwork( + ket::ITensorNetwork{T, V, I}, + map::Dictionary{I, I} + ) where {T, V, I} + check_input(NormNetwork, ket, map) + return new{T, V, I}(ket, map) + end +end + +function NormNetwork( + f::Function, tn::ITensorNetwork, linknames = ITensorNetworksNext.linknames(tn) + ) + return NormNetwork(tn, narrow_map(f, Indices(linknames))) +end +function NormNetwork(tn::ITensorNetwork, linknames = ITensorNetworksNext.linknames(tn)) + return NormNetwork(name -> derivedname(name, "bra"), tn, linknames) +end + +function check_input(::Type{<:NormNetwork}, ket, braname) + return check_branamemap(ket, braname) +end + +function Base.eltype(::Type{<:NormNetwork}) + return error( + "`eltype` of a `NormNetwork` is not defined, since the double-layer tensor at a vertex has no representation of its own. Use `kettensor` and `bratensor` to reach the individual layers." + ) +end + +# ==================================== DataGraphs.jl ===================================== # + +function DataGraphs.get_vertex_data(nn::NormNetwork, vertex) + return error( + "Indexing a `NormNetwork` is not defined, since the double-layer tensor at a vertex has no representation of its own. Use `kettensor` and `bratensor` to reach the individual layers." + ) +end + +# ====================================== interface ======================================= # + +ketnetwork(nn::NormNetwork) = nn.ket +branamemap(nn::NormNetwork) = nn.braname + +""" + flatten_network(nn::NormNetwork) -> ITensorNetwork + +Expand a norm network into a plain tensor network carrying one vertex per layer, so that vertex +`v` of `nn` becomes the two vertices `(v, :ket)` and `(v, :bra)` and `nv` doubles. A vertex's +two layers stay adjacent in the vertex order, ket first. +""" +function flatten_network(nn::NormNetwork) + return ITensorNetwork( + dictionary( + Iterators.flatten( + ((v, :ket) => kettensor(nn, v), (v, :bra) => bratensor(nn, v)) + for v in vertices(nn) + ) + ) + ) +end + +""" + normnetwork([f::Function], tn::ITensorNetwork, [linknames]) -> NormNetwork + normnetwork(tn::ITensorNetwork, braname::Dictionary) -> NormNetwork + +Build the double-layer norm network `⟨tn|tn⟩`, represented lazily as a `NormNetwork` object. +Each ket index name in `linknames` (by default `linknames(tn)`) gets the bra name `f(name)`, +where `f` defaults to `name -> derivedname(name, "bra")`. Alternatively pass the ket→bra name +map `braname` directly. Either way the names mapped must include every link name of `tn` and +may include others, and an error is thrown if a link name is missing, two bra names coincide +or a bra name is already an index name of `tn`. +""" +normnetwork(args...) = NormNetwork(args...) diff --git a/src/bilinearforms/quadraticformnetwork.jl b/src/bilinearforms/quadraticformnetwork.jl new file mode 100644 index 00000000..169b0b18 --- /dev/null +++ b/src/bilinearforms/quadraticformnetwork.jl @@ -0,0 +1,124 @@ +using Dictionaries: Dictionary, Indices +using ITensorBase: inputnames, outputnames, rename, state +using ITensorNetworksNext + +""" + struct QuadraticFormNetwork{T, V, I, O} <: AbstractBilinearFormNetwork{T, V, I} + +Lazy wrapper representing the quadratic form `⟨tn|op|tn⟩` of `tn::ITensorNetwork{T, V, I}` +sandwiched around the operator layer `op::ITensorNetworkOperator`, together with a per-index +ket→bra name mapping that, for each index in the ket layer, defines the name of the +corresponding index in the bra layer. + +The operator layer must be defined on every vertex of `tn`. Its input names are ket index +names, and each paired output name is renamed to the bra layer, so the operator's remaining +index names must be distinct from every index name of `tn`. +""" +struct QuadraticFormNetwork{T, V, I, O <: ITensorNetworkOperator} <: + AbstractBilinearFormNetwork{T, V, I} + ket::ITensorNetwork{T, V, I} + operator::O + braname::Dictionary{I, I} + function QuadraticFormNetwork( + ket::ITensorNetwork{T, V, I}, + operator::ITensorNetworkOperator, + map::Dictionary{I, I} + ) where {T, V, I} + check_input(QuadraticFormNetwork, ket, operator, map) + return new{T, V, I, typeof(operator)}(ket, operator, map) + end +end + +function QuadraticFormNetwork( + f::Function, + ket::ITensorNetwork, + operator::ITensorNetworkOperator, + linknames = ITensorNetworksNext.linknames(ket) + ) + names = union(linknames, inputnames(operator)) + return QuadraticFormNetwork(ket, operator, narrow_map(f, Indices(names))) +end +function QuadraticFormNetwork( + ket::ITensorNetwork, + operator::ITensorNetworkOperator, + linknames = ITensorNetworksNext.linknames(ket) + ) + return QuadraticFormNetwork(name -> derivedname(name, "bra"), ket, operator, linknames) +end + +function check_input(::Type{<:QuadraticFormNetwork}, ket, operator, braname) + if !issetequal(vertices(operator), vertices(ket)) + error("the operator layer must be defined on every vertex of the ket layer.") + end + inputs = inputnames(operator) + for input in inputs + if !has_dimname(ket, input) + error("operator input name $input is not an index name of the ket layer.") + end + if dimnamevertices(ket, input) != dimnamevertices(state(operator), input) + error("operator input name $input must be a ket site index on the same vertex.") + end + end + for vertex in vertices(operator), name in names(state(operator), vertex) + if name ∉ inputs && has_dimname(ket, name) + error("operator index name $name is not an input but is a ket index name.") + end + end + if !issubset(inputs, keys(braname)) + error("every operator input name must have a bra name.") + end + return check_branamemap(ket, braname, state(operator)) +end + +function Base.eltype(::Type{<:QuadraticFormNetwork}) + return error( + "`eltype` of a `QuadraticFormNetwork` is not defined, since the triple-layer tensor \ + at a vertex has no representation of its own. Use `kettensor`, `operatortensor` and \ + `bratensor` to reach the individual layers." + ) +end + +# ==================================== DataGraphs.jl ===================================== # + +function DataGraphs.get_vertex_data(qf::QuadraticFormNetwork, vertex) + return error( + "Indexing a `QuadraticFormNetwork` is not defined, since the triple-layer tensor at \ + a vertex has no representation of its own. Use `kettensor`, `operatortensor` and \ + `bratensor` to reach the individual layers." + ) +end + +function DataGraphs.is_vertex_assigned(qf::QuadraticFormNetwork, vertex) + return isassigned(ketnetwork(qf), vertex) && isassigned(operatornetwork(qf), vertex) +end + +# ====================================== interface ======================================= # + +ketnetwork(qf::QuadraticFormNetwork) = qf.ket +branamemap(qf::QuadraticFormNetwork) = qf.braname +operatornetwork(qf::QuadraticFormNetwork) = qf.operator + +# Each output name is renamed to the bra name of the input it is paired with, so the output legs +# meet the bra layer and the input legs meet the ket layer. +function operatortensor(qf::QuadraticFormNetwork, vertex) + op = qf.operator[vertex] + replacements = [ + output => braname(qf, input) for + (output, input) in zip(outputnames(op), inputnames(op)) + ] + return rename(state(op), replacements...) +end + +""" + quadraticformnetwork([f::Function], tn::ITensorNetwork, op::ITensorNetworkOperator, [linknames]) + quadraticformnetwork(tn::ITensorNetwork, op::ITensorNetworkOperator, braname::Dictionary) + +Build the triple-layer network `⟨tn|op|tn⟩`, represented lazily as a `QuadraticFormNetwork` +object. Each ket index name in `linknames` (by default `linknames(tn)`) and each input name of +`op` gets the bra name `f(name)`, where `f` defaults to `name -> derivedname(name, "bra")`. +Alternatively pass the ket→bra name map `braname` directly. Either way the names mapped must +include every link name of `tn` and every input name of `op` and may include others, and an +error is thrown if one is missing, two bra names coincide or a bra name is already an index +name of `tn` or `op`. +""" +quadraticformnetwork(args...) = QuadraticFormNetwork(args...) diff --git a/src/itensornetworkoperator.jl b/src/itensornetworkoperator.jl new file mode 100644 index 00000000..58e3e844 --- /dev/null +++ b/src/itensornetworkoperator.jl @@ -0,0 +1,119 @@ +using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned +using Dictionaries: Dictionaries, Dictionary +using Graphs: Graphs, AbstractGraph, edges, vertices +using ITensorBase: ITensorBase, NamedTensorOperator, inputnames, names, nametype, operator, + outputnames, state +using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex + +""" + struct ITensorNetworkOperator{T, V, I, P} <: AbstractITensorNetwork{T, V} + +The network equivalent of `ITensorBase.ITensorOperator`: a tensor network of type +`P <: AbstractITensorNetwork{T, V}` together with a pairing of its dangling index names, +stored as a map from each output name to its input name. Applying the operator contracts over +the input names and leaves the output names. + +The output and input of each pair must sit on the same vertex. Indexing returns the vertex +tensor wrapped as an `ITensorOperator` carrying the pairs at that vertex. +""" +struct ITensorNetworkOperator{T, V, I, P <: AbstractITensorNetwork{T, V}} <: + AbstractITensorNetwork{T, V} + parent::P + pairing::Dictionary{I, I} + function ITensorNetworkOperator( + parent::AbstractITensorNetwork{T, V}, outputnames, inputnames + ) where {T, V} + I = nametype(T) + outputnames = collect(I, outputnames) + inputnames = collect(I, inputnames) + if length(outputnames) != length(inputnames) + throw( + ArgumentError( + "Operator `outputnames` and `inputnames` must have equal length " * + "(positional pairing), got $(length(outputnames)) and " * + "$(length(inputnames))." + ) + ) + end + if !allunique(vcat(outputnames, inputnames)) + throw( + ArgumentError("each operator output and input name must appear only once.") + ) + end + for opname in Iterators.flatten((outputnames, inputnames)) + nvertices = length(dimnamevertices(parent, opname)) + if nvertices != 1 + throw( + ArgumentError( + "operator dim name $opname is associated with $nvertices vertices " * + "in the tensor network; an operator leg must be a dangling index." + ) + ) + end + end + for (output, input) in zip(outputnames, inputnames) + if only(dimnamevertices(parent, output)) != only(dimnamevertices(parent, input)) + throw( + ArgumentError( + "operator output $output and its paired input $input must sit on " * + "the same vertex." + ) + ) + end + end + pairing = Dictionary(outputnames, inputnames) + return new{T, V, I, typeof(parent)}(parent, pairing) + end +end + +function ITensorBase.operator(tn::AbstractITensorNetwork, outputnames, inputnames) + return ITensorNetworkOperator(tn, outputnames, inputnames) +end + +function Base.eltype(::Type{<:ITensorNetworkOperator{T, V, I}}) where {T, V, I} + return NamedTensorOperator{I, T} +end + +# ====================================== Graphs.jl ======================================= # + +Graphs.edges(op::ITensorNetworkOperator) = edges(state(op)) +Graphs.vertices(op::ITensorNetworkOperator) = vertices(state(op)) + +# ==================================== NamedGraphs.jl ==================================== # + +function NamedGraphs.encoded_vertex(op::ITensorNetworkOperator, vertex) + return encoded_vertex(state(op), vertex) +end +function NamedGraphs.decoded_vertex(op::ITensorNetworkOperator, code::Integer) + return decoded_vertex(state(op), code) +end +NamedGraphs.encoded_graph(op::ITensorNetworkOperator) = encoded_graph(state(op)) + +# ==================================== DataGraphs.jl ===================================== # + +# Both names of a pair sit on one vertex, so the pairs whose output is on `vertex` are its pairing. +function DataGraphs.get_vertex_data(op::ITensorNetworkOperator, vertex) + tensor = state(op)[vertex] + outputs = filter(name -> haskey(op.pairing, name), names(tensor)) + return operator(tensor, outputs, [op.pairing[output] for output in outputs]) +end + +function DataGraphs.is_vertex_assigned(op::ITensorNetworkOperator, vertex) + return is_vertex_assigned(state(op), vertex) +end + +# =================================== Dictionaries.jl ==================================== # + +Dictionaries.issettable(::ITensorNetworkOperator) = false +Dictionaries.isinsertable(::ITensorNetworkOperator) = false + +# ====================================== interface ======================================= # + +ITensorBase.state(op::ITensorNetworkOperator) = op.parent +Base.parent(op::ITensorNetworkOperator) = state(op) +ITensorBase.outputnames(op::ITensorNetworkOperator) = collect(keys(op.pairing)) +ITensorBase.inputnames(op::ITensorNetworkOperator) = collect(op.pairing) + +function operator_support(tn::AbstractGraph, op::ITensorNetworkOperator) + return operator_support_names(tn, inputnames(op)) +end diff --git a/src/normnetwork.jl b/src/normnetwork.jl deleted file mode 100644 index d735b22f..00000000 --- a/src/normnetwork.jl +++ /dev/null @@ -1,113 +0,0 @@ -using Dictionaries: Dictionary, dictionary -using ITensorBase: rename, 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} - ket::ITensorNetwork{T, V, I} - braname::Dictionary{I, I} - function NormNetwork( - ket::ITensorNetwork{T, V, I}, - map::Dictionary{I, I} - ) where {T, V, I} - braname = Dictionary{I, I}() - for (name, vertices) in pairs(ket.dimname_vertices) - if length(vertices) == 2 - insert!(braname, name, map[name]) - end - end - return new{T, V, I}(ket, braname) - end -end - -function Base.eltype(::Type{<:NormNetwork}) - return error( - "`eltype` of a `NormNetwork` is not defined, since the double-layer tensor at a vertex has no representation of its own. Use `kettensor` and `bratensor` to reach the individual layers." - ) -end - -function NormNetwork(tn::ITensorNetwork) - return NormNetwork(tn, map(uniquename, keys(tn.dimname_vertices))) -end - -# ====================================== Graphs.jl ======================================= # - -Graphs.edges(nn::NormNetwork) = edges(nn.ket) -Graphs.vertices(nn::NormNetwork) = vertices(nn.ket) - -# ==================================== NamedGraphs.jl ==================================== # - -NamedGraphs.encoded_vertex(nn::NormNetwork, vertex) = encoded_vertex(nn.ket, vertex) -NamedGraphs.decoded_vertex(nn::NormNetwork, code::Integer) = decoded_vertex(nn.ket, code) -NamedGraphs.encoded_graph(nn::NormNetwork) = encoded_graph(nn.ket) - -# ==================================== DataGraphs.jl ===================================== # - -function DataGraphs.get_vertex_data(nn::NormNetwork, vertex) - return error( - "Indexing a `NormNetwork` is not defined, since the double-layer tensor at a vertex has no representation of its own. Use `kettensor` and `bratensor` to reach the individual layers." - ) -end - -function DataGraphs.is_vertex_assigned(nn::NormNetwork, vertex) - return isassigned(nn.ket, vertex) -end -# =================================== Dictionaries.jl ==================================== # - -Dictionaries.issettable(::NormNetwork) = false -Dictionaries.isinsertable(::NormNetwork) = false - -# ====================================== interface ======================================= # - -function braname(nn::NormNetwork, name) - if !has_dimname(nn.ket, name) - error("index name $name not found underlying tensor network.") - end - # The indices not stored in `nn.braname` are precisely the site indices, which - # get mapped to themselves. - return get(nn.braname, name, name) -end - -indmap(nn::NormNetwork, ind) = setname(conj(ind), braname(nn, name(ind))) - -kettensor(nn::NormNetwork, vertex) = nn.ket[vertex] -function conj_bratensor(nn::NormNetwork, vertex) - return rename(n -> braname(nn, n), kettensor(nn, vertex)) -end - -bratensor(nn::NormNetwork, vertex) = conj(conj_bratensor(nn, vertex)) - -""" - flatten_network(nn::NormNetwork) -> ITensorNetwork - -Expand a norm network into a plain tensor network carrying one vertex per layer, so that vertex -`v` of `nn` becomes the two vertices `(v, :ket)` and `(v, :bra)` and `nv` doubles. A vertex's -two layers stay adjacent in the vertex order, ket first. -""" -function flatten_network(nn::NormNetwork) - return ITensorNetwork( - dictionary( - Iterators.flatten( - ((v, :ket) => kettensor(nn, v), (v, :bra) => bratensor(nn, v)) - for v in vertices(nn) - ) - ) - ) -end - -""" - normnetwork(tn::ITensorNetwork, [braname]) -> NormNetwork - -Build the double-layer norm network `⟨tn|tn⟩`, represented lazily as a `NomnNetwork` object. -The optional second argument `braname` should implement `braname[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, braname) = NormNetwork(tn, braname) diff --git a/src/normnetworkview.jl b/src/normnetworkview.jl deleted file mode 100644 index 75c57926..00000000 --- a/src/normnetworkview.jl +++ /dev/null @@ -1,58 +0,0 @@ -using DataGraphs: DataGraphs, get_vertex_data, is_vertex_assigned -using Dictionaries: Dictionaries, isinsertable, issettable -using Graphs: Graphs, edges, vertices -using NamedGraphs: NamedGraphs, decoded_vertex, encoded_graph, encoded_vertex - -struct KetView{T, V, I} <: 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.encoded_vertex(nnv::$View, vertex) - return encoded_vertex(nnv.parent, vertex) - end - function NamedGraphs.decoded_vertex(nnv::$View, code::Integer) - return decoded_vertex(nnv.parent, code) - end - - NamedGraphs.encoded_graph(nnv::$View) = encoded_graph(nnv.parent) - end -end - -# ==================================== DataGraphs.jl ===================================== # - -DataGraphs.get_vertex_data(nn::KetView, vertex) = kettensor(nn.parent, vertex) -DataGraphs.get_vertex_data(nn::BraView, vertex) = bratensor(nn.parent, vertex) - -for View in (:KetView, :BraView) - @eval begin - function DataGraphs.is_vertex_assigned(nnv::$View, vertex) - return isassigned(nnv.parent.ket, vertex) - end - 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 57ad5b2b..1c972cb2 100644 --- a/src/tensornetwork.jl +++ b/src/tensornetwork.jl @@ -146,7 +146,7 @@ function DataGraphs.set_vertex_data!(tn::ITensorNetwork, tensor, vertex) end function update_tensornetwork_metadata!(tn, vertex, tensor) - oldnames = isassigned(tn, vertex) ? names(tn[vertex]) : Set{nametype(tn)}() + oldnames = isassigned(tn, vertex) ? names(tn, vertex) : Set{nametype(tn)}() newnames = names(tensor) update_tensornetwork_metadata!(tn, vertex, oldnames, newnames) @@ -197,12 +197,19 @@ Graphs.add_edge!(::ITensorNetwork, _edge) = false # PERF: fast lookup compared to `AbstractITensorNetwork` fallback. function dimnamevertices(tn::ITensorNetwork, name) - return get(tn.dimname_vertices, name, Set{vertextype(tn)}()) + return copy(get(tn.dimname_vertices, name, Set{vertextype(tn)}())) end # PERF: fast lookup compared to `AbstractITensorNetwork` fallback. has_dimname(tn::ITensorNetwork, name) = haskey(tn.dimname_vertices, name) +# PERF: fast lookup compared to `AbstractITensorNetwork` fallback. +function linknames(tn::ITensorNetwork) + return [ + name for (name, vertices) in pairs(tn.dimname_vertices) if length(vertices) == 2 + ] +end + function NamedGraphs.similar_graph( T::Type{<:ITensorNetwork}, vertices = vertextype(T)[] diff --git a/src/utils.jl b/src/utils.jl index abf666e4..5113f317 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,6 +1,21 @@ using Dictionaries: AbstractIndices, Dictionary +using ITensorBase: ITensorBase, IndexName, settags +using UUIDs: uuid5 # `map` does something clever to figure out the element type of the output when not # inferable, but the `map` overload on `Dictionary` does not, so we fix this here. narrow_map(f, v) = map(f, v) narrow_map(f, v::AbstractIndices) = Dictionary(v, [f(x) for x in v]) + +""" + derivedname(name::IndexName, key) -> IndexName + +Return a name derived deterministically from `name` and `key`: the identifier is replaced by +the name-based UUID `uuid5(uuid(name), string(key))`, `key` is added as a bare-label tag, and +the other tags and the prime level are kept. Every process derives the same name from the same +`name` and `key`, and the result differs from `name` and from any independently minted name. +""" +function derivedname(name::IndexName, key) + derived = ITensorBase.setuuid(name, uuid5(ITensorBase.uuid(name), string(key))) + return settags(derived, string(key)) +end diff --git a/test/test_itensornetworkoperator.jl b/test/test_itensornetworkoperator.jl new file mode 100644 index 00000000..a2918297 --- /dev/null +++ b/test/test_itensornetworkoperator.jl @@ -0,0 +1,115 @@ +using DataGraphs: is_vertex_assigned +using Dictionaries: isinsertable, issettable +using Graphs: edges, vertices +using ITensorBase: ITensor, Index, inputnames, name, operator, outputnames, state +using ITensorNetworksNext: + ITensorNetwork, ITensorNetworkOperator, operator_support, tensornetwork +using NamedGraphs: incident_edges, named_path_graph +using Test: @test, @test_throws, @testset + +# A bondless operator network on `g`, with one input/output pair per vertex. +function local_operator(g; d = 2) + inp = Dict(v => Index(d) for v in vertices(g)) + out = Dict(v => Index(d) for v in vertices(g)) + tn = tensornetwork(vertices(g)) do v + return randn((out[v], inp[v])) + end + vs = collect(vertices(g)) + op = operator(tn, [name(out[v]) for v in vs], [name(inp[v]) for v in vs]) + return op, inp, out +end + +@testset "`ITensorNetworkOperator`" begin + @testset "Basics" begin + g = named_path_graph(3) + op, inp, out = local_operator(g) + + @test op isa ITensorNetworkOperator + + # The operator shares the graph structure of the network it wraps. + @test issetequal(vertices(op), vertices(state(op))) + @test issetequal(edges(op), edges(state(op))) + @test parent(op) === state(op) + + # The pairing is positional: `outputnames[i]` goes with `inputnames[i]`. + vs = collect(vertices(g)) + @test outputnames(op) == [name(out[v]) for v in vs] + @test inputnames(op) == [name(inp[v]) for v in vs] + + @test is_vertex_assigned(op, first(vs)) + + # It is a lazy wrapper, so it is neither settable nor insertable. + @test !issettable(op) + @test !isinsertable(op) + end + + @testset "`getindex` wraps the vertex tensor" begin + g = named_path_graph(3) + op, inp, out = local_operator(g) + + # Both halves of each pair sit on the same vertex here, so the wrapper carries them. + for v in vertices(g) + @test outputnames(op[v]) == [name(out[v])] + @test inputnames(op[v]) == [name(inp[v])] + @test state(op[v]) === state(op)[v] + end + + # `eltype` is the type of the wrapped vertex data. + @test eltype(op) === typeof(op[first(vertices(g))]) + end + + @testset "a pair must sit on one vertex" begin + # A swap pairs the output at vertex 1 with the input at vertex 2, and vice versa. + a, b = Index(2), Index(2) + a′, b′ = Index(2), Index(2) + tn = ITensorNetwork(Dict(1 => randn((a′, a)), 2 => randn((b′, b)))) + @test_throws ArgumentError operator(tn, [name(a′), name(b′)], [name(b), name(a)]) + end + + @testset "constructor validation" begin + g = named_path_graph(2) + op, inp, out = local_operator(g) + tn = state(op) + vs = collect(vertices(g)) + + # The pairing is a bijection, so the two name lists must have equal length. + @test_throws ArgumentError operator(tn, [name(out[vs[1]])], []) + @test_throws ArgumentError operator( + tn, [name(out[v]) for v in vs], [name(inp[vs[1]])] + ) + # An output name can be paired with only one input. + @test_throws ArgumentError operator( + tn, [name(out[vs[1]]), name(out[vs[1]])], [name(inp[vs[1]]), name(inp[vs[1]])] + ) + # An input name can be paired with only one output, and no name is both. + @test_throws ArgumentError operator( + tn, [name(out[vs[1]]), name(out[vs[2]])], [name(inp[vs[1]]), name(inp[vs[1]])] + ) + @test_throws ArgumentError operator( + tn, [name(out[vs[1]]), name(inp[vs[1]])], [name(inp[vs[1]]), name(out[vs[2]])] + ) + + # An operator leg must be dangling: a link name touches two vertices. + l = Index(2) + linked = tensornetwork(vertices(g)) do v + return randn((inp[v], l)) + end + @test_throws ArgumentError operator(linked, [name(l)], [name(inp[vs[1]])]) + end + + @testset "`operator_support`" begin + g = named_path_graph(3) + s = Dict(v => Index(2) for v in vertices(g)) + state_tn = tensornetwork(vertices(g)) do v + is = map(e -> Index(2), incident_edges(g, v)) + return randn((s[v], is...)) + end + + # An operator acting on sites 1 and 3 is supported on exactly those vertices. + out1, out3 = Index(2), Index(2) + optn = ITensorNetwork(Dict(1 => randn((out1, s[1])), 3 => randn((out3, s[3])))) + op = operator(optn, [name(out1), name(out3)], [name(s[1]), name(s[3])]) + + @test issetequal(operator_support(state_tn, op), Set([1, 3])) + end +end diff --git a/test/test_normnetwork.jl b/test/test_normnetwork.jl index 4528b220..b5109561 100644 --- a/test/test_normnetwork.jl +++ b/test/test_normnetwork.jl @@ -1,10 +1,11 @@ using DataGraphs: is_vertex_assigned -using Dictionaries: isinsertable, issettable +using Dictionaries: Indices, isinsertable, issettable using Graphs: edges, vertices -using ITensorBase: ITensor, Index, IndexName, conj, inds, name, setname, uniquename -using ITensorNetworksNext: BraView, Greedy, ITensorNetwork, KetView, NormNetwork, braname, - bratensor, conj_bratensor, contraction_tree, flatten_network, indmap, kettensor, - normnetwork, prod_tensors, tensornetwork +using ITensorBase: + ITensor, Index, IndexName, inds, name, nametype, plev, prime, setname, tags, uniquename +using ITensorNetworksNext: BraView, Greedy, ITensorNetwork, NormNetwork, braname, + branetwork, bratensor, conj_bratensor, contraction_tree, derivedname, flatten_network, + ketnetwork, kettensor, linknames, normnetwork, prod_tensors, tensornetwork using LinearAlgebra: norm using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph using Test: @test, @test_throws, @testset @@ -44,6 +45,7 @@ contract_norm(nn) = prod_tensors(flatten_network(nn), Greedy())[] # `eltype` is defined. @test_throws ErrorException nn[1] @test_throws ErrorException eltype(nn) + @test nametype(nn) === nametype(tn) === IndexName # Vertex data is assigned wherever the underlying network is. @test is_vertex_assigned(nn, 1) @@ -79,11 +81,6 @@ contract_norm(nn) = prod_tensors(flatten_network(nn), Greedy())[] # `bra` is the elementwise conjugate of `conj_bratensor` and carries the same indices. @test inds(bratensor(nn, 2)) == inds(conj_bratensor(nn, 2)) - # `indmap` conjugates an index and renames it according to the name map. - ind = only(i for i in inds(kettensor(nn, 2)) if name(i) == lname) - @test name(indmap(nn, ind)) == braname(nn, name(ind)) - @test indmap(nn, ind) == setname(conj(ind), braname(nn, name(ind))) - # Querying the name map with an index name absent from the network errors. @test_throws ErrorException braname(nn, name(Index(2))) end @@ -93,7 +90,7 @@ contract_norm(nn) = prod_tensors(flatten_network(nn), Greedy())[] 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)) + custom = map(uniquename, Indices(linknames(tn))) nn = normnetwork(tn, custom) lname = name(l[NamedEdge(1 => 2)]) @@ -101,32 +98,71 @@ contract_norm(nn) = prod_tensors(flatten_network(nn), Greedy())[] @test braname(nn, lname) in name.(inds(conj_bratensor(nn, 2))) end - @testset "`KetView` / `BraView`" begin + @testset "`derivedname`" begin + n = prime(name(Index(2; tags = "Link"))) + + # The derived name depends only on the name and key, keeps the decoration and is tagged + # with the key. + @test derivedname(n, :bra) == derivedname(n, :bra) + @test derivedname(n, :bra) != n + @test derivedname(n, :bra) != derivedname(n, :other) + @test derivedname(n, :bra) != prime(n) + @test tags(derivedname(n, :bra)) == merge(tags(n), Dict("bra" => "")) + @test plev(derivedname(n, :bra)) == plev(n) + end + + @testset "default and function name maps" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + lname = name(l[NamedEdge(1 => 2)]) + + # Two independently constructed norm networks agree on every bra name. + @test braname(NormNetwork(tn), lname) == derivedname(lname, :bra) + @test braname(normnetwork(tn), lname) == braname(NormNetwork(tn), lname) + @test braname(normnetwork(prime, tn), lname) == prime(lname) + + # A bra name that is already an index name of the network is rejected. + taken = setname(Index(2), derivedname(lname, :bra)) + tn[1] = tn[1] * ITensor(ones(2), (taken,)) + @test_throws ErrorException NormNetwork(tn) + end + + @testset "bond to a vertex outside the network" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + + # The part of `tn` on vertices 1 and 2, as one rank of a distributed network holds it. + part = ITensorNetwork(Dict(v => tn[v] for v in (1, 2))) + nn = normnetwork(part, map(uniquename, Indices(linknames(tn)))) + + # The bond to vertex 3 is not a link of `part`, but it is in the map, so it gets its + # own bra name and the two layers at vertex 2 share only the site index. + lname = name(l[NamedEdge(2 => 3)]) + @test braname(nn, lname) != lname + shared = + intersect(name.(inds(kettensor(nn, 2))), name.(inds(conj_bratensor(nn, 2)))) + @test issetequal(shared, [name(s[2])]) + end + + @testset "`ketnetwork` / `branetwork`" begin g = named_path_graph(3) tn, l, s = random_state(Float64, g) nn = NormNetwork(tn) - kv = KetView(nn) - bv = BraView(nn) + # The ket layer is the network the norm network was built from. + @test ketnetwork(nn) === tn - # Views share the graph structure of the underlying network. - @test issetequal(vertices(kv), vertices(tn)) + # The bra layer is not stored, so it is a view sharing the ket layer's graph structure. + bv = branetwork(nn) + @test bv isa BraView @test issetequal(vertices(bv), vertices(tn)) - @test issetequal(edges(kv), edges(tn)) @test issetequal(edges(bv), edges(tn)) - - # The ket view exposes the bare ket tensors; the bra view exposes the bra tensors. for v in vertices(tn) - @test kv[v] === kettensor(nn, v) @test inds(bv[v]) == inds(bratensor(nn, v)) end - - @test is_vertex_assigned(kv, 1) @test is_vertex_assigned(bv, 1) - # Views inherit the (non-)mutability of their parent norm network. - @test !issettable(kv) - @test !isinsertable(kv) + # The view inherits the (non-)mutability of its parent norm network. @test !issettable(bv) @test !isinsertable(bv) end @@ -181,7 +217,7 @@ contract_norm(nn) = prod_tensors(flatten_network(nn), Greedy())[] @test contract_norm(NormNetwork(tn)) ≈ 1 # The contracted norm does not depend on the chosen bra-layer name map. - custom = map(uniquename, keys(tn.dimname_vertices)) + custom = map(uniquename, Indices(linknames(tn))) @test contract_norm(normnetwork(tn, custom)) ≈ contract_norm(NormNetwork(tn)) end end diff --git a/test/test_quadraticformnetwork.jl b/test/test_quadraticformnetwork.jl new file mode 100644 index 00000000..cfcbdc5d --- /dev/null +++ b/test/test_quadraticformnetwork.jl @@ -0,0 +1,258 @@ +using DataGraphs: is_vertex_assigned +using Dictionaries: isinsertable, issettable +using Graphs: edges, vertices +using ITensorBase: ITensor, Index, IndexName, conj, inds, inputnames, name, names, nametype, + operator, outputnames, rename, setname, uniquename +using ITensorNetworksNext: BraView, Greedy, ITensorNetwork, NormNetwork, + QuadraticFormNetwork, braname, branetwork, bratensor, conj_bratensor, derivedname, + factor_tensors, ketnetwork, kettensor, operatornetwork, operatortensor, prod_tensors, + quadraticformnetwork, tensornetwork +using LinearAlgebra: I, norm +using NamedGraphs: NamedEdge, incident_edges, named_grid, named_path_graph +using Test: @test, @test_throws, @testset + +# Build a random `ITensorNetwork` state on the graph `g` with site dimension `d` and +# bond dimension `χ`. +function random_state(::Type{T}, g; d = 2, χ = 2) where {T} + l = Dict(e => Index(χ) for e in edges(g)) + l = merge(l, Dict(reverse(e) => l[e] for e in edges(g))) + s = Dict(v => Index(d) for v in vertices(g)) + tn = tensornetwork(vertices(g)) do v + is = map(e -> l[e], incident_edges(g, v)) + return randn(T, (s[v], is...)) + end + return tn, l, s +end + +# Build a bondless operator layer on the vertices of `g`, with `f(v)` supplying the matrix +# acting on the site index `s[v]`. +function product_operator(f, g, s; d = 2) + out = Dict(v => Index(d) for v in vertices(g)) + tn = tensornetwork(vertices(g)) do v + return ITensor(f(v), (out[v], s[v])) + end + vs = collect(vertices(g)) + return operator(tn, [name(out[v]) for v in vs], [name(s[v]) for v in vs]) +end + +identity_operator(g, s; d = 2) = product_operator(v -> Matrix(1.0I, d, d), g, s; d) + +function contract(bn) + return prod_tensors(reduce(vcat, factor_tensors(bn, v) for v in vertices(bn)), Greedy()) +end + +@testset "`QuadraticFormNetwork`" begin + @testset "Basics" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + op = identity_operator(g, s) + qf = QuadraticFormNetwork(tn, op) + + # `quadraticformnetwork` is the public constructor and agrees with the type. + @test quadraticformnetwork(tn, op) isa QuadraticFormNetwork + @test qf isa QuadraticFormNetwork + @test nametype(qf) === nametype(tn) === IndexName + + # The quadratic form shares the graph structure of the ket layer. + @test issetequal(vertices(qf), vertices(tn)) + @test issetequal(edges(qf), edges(tn)) + + # The triple-layer tensor at a vertex has no representation, so neither indexing nor + # `eltype` is defined. + @test_throws ErrorException qf[1] + @test_throws ErrorException eltype(qf) + + # Vertex data is assigned wherever both layers are. + @test is_vertex_assigned(qf, 1) + + # The quadratic form is neither settable nor insertable (it is a lazy view). + @test !issettable(qf) + @test !isinsertable(qf) + + # The operator layer must cover every vertex of the ket layer. + g4 = named_path_graph(4) + _, _, s4 = random_state(Float64, g4) + @test_throws ErrorException QuadraticFormNetwork(tn, identity_operator(g4, s4)) + end + + @testset "operator input outside the ket network" begin + g = named_path_graph(2) + tn, l, s = random_state(Float64, g) + stray = Index(2) + out1, out2 = Index(2), Index(2) + optn = ITensorNetwork(Dict(1 => randn((out1, stray)), 2 => randn((out2, s[2])))) + op = operator(optn, [name(out1), name(out2)], [name(stray), name(s[2])]) + @test_throws ErrorException QuadraticFormNetwork(tn, op) + end + + @testset "operator input on a different vertex of the ket network" begin + g = named_path_graph(2) + tn, l, s = random_state(Float64, g) + out1, out2 = Index(2), Index(2) + optn = ITensorNetwork(Dict(1 => randn((out1, s[2])), 2 => randn((out2, s[1])))) + op = operator(optn, [name(out1), name(out2)], [name(s[2]), name(s[1])]) + @test_throws ErrorException QuadraticFormNetwork(tn, op) + end + + @testset "operator input equal to a ket link name" begin + g = named_path_graph(2) + tn, l, s = random_state(Float64, g) + e = only(edges(g)) + out1, out2 = Index(2), Index(2) + optn = ITensorNetwork(Dict(1 => randn((out1, l[e])), 2 => randn((out2, s[2])))) + op = operator(optn, [name(out1), name(out2)], [name(l[e]), name(s[2])]) + @test_throws ErrorException QuadraticFormNetwork(tn, op) + end + + @testset "operator link name equal to a ket link name" begin + g = named_path_graph(2) + tn, l, s = random_state(Float64, g) + e = only(edges(g)) + out1, out2 = Index(2), Index(2) + optn = ITensorNetwork( + Dict(1 => randn((out1, s[1], l[e])), 2 => randn((out2, s[2], l[e]))) + ) + op = operator(optn, [name(out1), name(out2)], [name(s[1]), name(s[2])]) + @test_throws ErrorException QuadraticFormNetwork(tn, op) + end + + @testset "a default bra name equal to an operator name" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + + # Each operator output takes the name the default map would derive for its input. + out = + Dict(v => setname(Index(2), derivedname(name(s[v]), :bra)) for v in vertices(g)) + optn = tensornetwork(v -> ITensor(Matrix(1.0I, 2, 2), (out[v], s[v])), vertices(g)) + vs = collect(vertices(g)) + op = operator(optn, [name(out[v]) for v in vs], [name(s[v]) for v in vs]) + + @test_throws ErrorException QuadraticFormNetwork(tn, op) + end + + @testset "layer tensors and the name map" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + op = identity_operator(g, s) + qf = QuadraticFormNetwork(tn, op) + + # `kettensor` returns the underlying tensor untouched. + @test kettensor(qf, 2) === tn[2] + + # Unlike the norm network, the site indices *are* renamed in the bra layer: the + # operator sits between the two layers, so they no longer contract directly. + sname = name(s[2]) + @test braname(qf, sname) != sname + @test sname in name.(inds(kettensor(qf, 2))) + @test !(sname in name.(inds(conj_bratensor(qf, 2)))) + @test braname(qf, sname) in name.(inds(conj_bratensor(qf, 2))) + + # Link indices are shared by two tensors, so they are renamed in the bra layer to + # keep the two layers' bonds distinct. + lname = name(l[NamedEdge(1 => 2)]) + @test braname(qf, lname) != lname + @test lname in name.(inds(kettensor(qf, 2))) + @test !(lname in name.(inds(conj_bratensor(qf, 2)))) + @test braname(qf, lname) in name.(inds(conj_bratensor(qf, 2))) + + # The operator's input name meets the ket and its output name is renamed to meet + # the bra. + o = operatortensor(qf, 2) + @test sname in names(o) + @test braname(qf, sname) in names(o) + @test inputnames(op) == [name(s[v]) for v in vertices(g)] + + # `bratensor` is the elementwise conjugate of `conj_bratensor` and carries the same + # indices. + @test inds(bratensor(qf, 2)) == inds(conj_bratensor(qf, 2)) + + # Querying the name map with an index name absent from the ket layer errors. + @test_throws ErrorException braname(qf, name(Index(2))) + + # `factor_tensors` returns the three layers as separate operands. + @test factor_tensors(qf, 2) == + [kettensor(qf, 2), operatortensor(qf, 2), bratensor(qf, 2)] + end + + @testset "custom name map" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + op = identity_operator(g, s) + + # A user-supplied map dictates the bra-layer name for each renamed index. + custom = map(uniquename, keys(tn.dimname_vertices)) + qf = quadraticformnetwork(tn, op, custom) + + lname = name(l[NamedEdge(1 => 2)]) + @test braname(qf, lname) == custom[lname] + @test braname(qf, lname) in name.(inds(conj_bratensor(qf, 2))) + + sname = name(s[2]) + @test braname(qf, sname) == custom[sname] + end + + @testset "`ketnetwork` / `operatornetwork` / `branetwork`" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + op = identity_operator(g, s) + qf = QuadraticFormNetwork(tn, op) + + # The ket and operator layers are the networks the quadratic form was built from. + @test ketnetwork(qf) === tn + @test operatornetwork(qf) === op + + # The bra layer is not stored, so it is a view sharing the ket layer's graph structure. + bv = branetwork(qf) + @test bv isa BraView + @test issetequal(vertices(bv), vertices(tn)) + @test issetequal(edges(bv), edges(tn)) + @test !issettable(bv) + @test !isinsertable(bv) + @test is_vertex_assigned(bv, 1) + for v in vertices(tn) + @test inds(bv[v]) == inds(bratensor(qf, v)) + end + end + + @testset "contraction / physics" begin + @testset "identity operator layer reproduces the norm network" begin + g = named_grid((2, 2)) + tn, l, s = random_state(Float64, g) + op = identity_operator(g, s) + + @test contract(QuadraticFormNetwork(tn, op))[] ≈ + contract(NormNetwork(tn))[] + end + + @testset "$T" for T in (Float64, ComplexF64) + g = named_grid((2, 2)) + tn, l, s = random_state(T, g) + + # A product of on-site matrices, contracted densely for the reference value. + mats = Dict(v => randn(T, 2, 2) for v in vertices(g)) + op = product_operator(v -> mats[v], g, s) + qf = QuadraticFormNetwork(tn, op) + + # ⟨ψ|O|ψ⟩ built by applying each on-site matrix to the dense ket and + # overlapping with the dense bra. + psi = prod(tn) + ket = psi + for v in vertices(g) + out = Index(2) + gate = ITensor(mats[v], (out, s[v])) + ket = rename(gate * ket, name(out) => name(s[v])) + end + @test contract(qf)[] ≈ (conj(psi) * ket)[] + end + + @testset "scaled identity" begin + g = named_path_graph(3) + tn, l, s = random_state(Float64, g) + op = product_operator(v -> 2.0 * Matrix(1.0I, 2, 2), g, s) + + # A factor of 2 at each of the three vertices scales the norm by 2³. + @test contract(QuadraticFormNetwork(tn, op))[] ≈ + 8 * contract(NormNetwork(tn))[] + end + end +end diff --git a/test/test_tensornetwork.jl b/test/test_tensornetwork.jl index 0beaad01..7a2bee94 100644 --- a/test/test_tensornetwork.jl +++ b/test/test_tensornetwork.jl @@ -1,7 +1,7 @@ using DataGraphs: DataGraph, assigned_edge_data, assigned_vertex_data, underlying_graph, vertex_data -using Graphs: add_edge!, add_vertex!, dst, edges, edgetype, has_edge, has_vertex, - is_directed, ne, nv, rem_edge!, rem_vertex!, src, vertices +using Graphs: AbstractGraph, add_edge!, add_vertex!, dst, edges, edgetype, has_edge, + has_vertex, is_directed, ne, nv, rem_edge!, rem_vertex!, src, vertices using ITensorBase: Index, inds, operator using ITensorNetworksNext: ITensorNetwork, has_ind, linkaxes, linkinds, linknames, operator_support, siteaxes, siteinds, sitenames, tensornetwork @@ -119,6 +119,8 @@ using Test: @test, @test_throws, @testset @test linknames(tn, 1 => 2) == [l[E(1 => 2)].name] @test linknames(tn, E(1 => 2)) == [l[E(1 => 2)].name] + @test issetequal(linknames(tn), [l[e].name for e in edges(g)]) + @test issetequal(invoke(linknames, Tuple{AbstractGraph}, tn), linknames(tn)) @test siteinds(tn, 1) == [s[1]] @test siteaxes(tn, 2) == [s[2]]