Add support for tensor network forms (bilinear and quadratic)

src/Forms/abstractbilinearform.jl

@@ -0,0 +1,14 @@
+abstract type AbstractForm{V,BraMap,KetMap,OperatorMap,SpaceMap} <:
+              AbstractITensorNetwork{V} end
+
+#Needed for interface
+bra_map(f::AbstractForm) = not_implemented()
+ket_map(f::AbstractForm) = not_implemented()
+operator_map(f::AbstractForm) = not_implemented()
+space_map(f::AbstractForm) = not_implemented()
+tensornetwork(f::AbstractForm) = not_implemented()
+copy(f::AbstractForm) = not_implemented()
+
+bra(f::AbstractForm) = induced_subgraph(f, collect(values(bra_map(f))))
+ket(f::AbstractForm) = induced_subgraph(f, collect(values(ket_map(f))))
+operator(f::AbstractForm) = induced_subgraph(f, collect(values(operator_map(f))))

src/Forms/bilinearform.jl

@@ -0,0 +1,42 @@
+default_index_map = prime
+default_bra_vertex_map(v) = (v, "bra")
+default_ket_vertex_map(v) = (v, "ket")
+default_operator_vertex_map(v) = (v, "operator")
+default_tno_constructor(s::IndsNetwork) = delta_network(s)
+
+struct BilinearForm{V,KetMap,BraMap,OperatorMap,SpaceMap} <:
+       AbstractForm{V,KetMap,BraMap,OperatorMap,SpaceMap}
+  tn::ITensorNetwork{V}
+  bra_vs_map::BraMap
+  ket_vs_map::KetMap
+  operator_vs_map::OperatorMap
+  dual_space_map::SpaceMap
+end
+
+#Needed for implementation
+bra_map(blf::BilinearForm) = blf.bra_vs_map
+ket_map(blf::BilinearForm) = blf.ket_vs_map
+operator_map(blf::BilinearForm) = blf.operator_vs_map
+space_map(blf::BilinearForm) = blf.dual_space_map
+tensornetwork(blf::BilinearForm) = blf.tn
+data_graph_type(::Type{<:BilinearForm}) = data_graph_type(tensornetwork(blf))
+data_graph(blf::BilinearForm) = data_graph(tensornetwork(blf))
+
+function BilinearForm(
+  bra::ITensorNetwork, operator::ITensorNetwork, ket::ITensorNetwork; kwargs...
+)
+  return Form(bra, operator, ket; form_type="Bilinear", kwargs...)
+end
+
+function BilinearForm(
+  bra::ITensorNetwork,
+  ket::ITensorNetwork;
+  dual_space_map=default_index_map,
+  tno_constructor=default_tno_constructor,
+  kwargs...,
+)
+  s = siteinds(bra)
+  operator_space = union_all_inds(s, dual_space_map(s; links=[]))
+  O = tno_constructor(operator_space)
+  return BilinearForm(bra, O, ket; kwargs...)
+end

src/Forms/construct_form.jl

@@ -0,0 +1,44 @@
+default_form_type = "Bilinear"
+
+function Form(
+  bra::ITensorNetwork,
+  operator::ITensorNetwork,
+  ket::ITensorNetwork;
+  dual_space_map=default_index_map,
+  link_space_map=default_index_map,
+  bra_vertex_map=default_bra_vertex_map,
+  ket_vertex_map=default_ket_vertex_map,
+  operator_vertex_map=default_operator_vertex_map,
+  form_type=default_form_type,
+)
+  dual_map = map(dual_space_map, Indices(externalinds(ket)))
+
+  @assert Set(externalinds(bra)) == Set(externalinds(ket))
+  @assert Set(externalinds(operator)) ==
+    union(Set(keys(dual_map)), Set(collect(values(dual_map))))
+  @assert isempty(findall(in(internalinds(bra)), internalinds(ket)))
+  @assert isempty(findall(in(internalinds(bra)), internalinds(operator)))
+  @assert isempty(findall(in(internalinds(ket)), internalinds(operator)))
+
+  bra_vs_map = map(bra_vertex_map, Indices(vertices(bra)))
+  ket_vs_map = map(ket_vertex_map, Indices(vertices(ket)))
+  operator_vs_map = map(operator_vertex_map, Indices(vertices(operator)))
+
+  bra_renamed = rename_vertices_itn(bra, bra_vs_map)
+  ket_renamed = rename_vertices_itn(ket, ket_vs_map)
+  operator_renamed = rename_vertices_itn(operator, operator_vs_map)
+  bra_renamed_dual = map_inds(dual_space_map, bra_renamed; links=[])
+
+  tn = union(union(bra_renamed_dual, operator_renamed), ket_renamed)
+
+  if form_type == "Bilinear"
+    return BilinearForm(tn, bra_vs_map, ket_vs_map, operator_vs_map, dual_map)
+  elseif form_type == "Quadratic"
+    link_map = map(link_space_map, Indices(internalinds(ket)))
+    return QuadraticForm(
+      tn, bra_vs_map, ket_vs_map, operator_vs_map, dual_map, link_space_map
+    )
+  else
+    return error("Form type not supported")
+  end
+end

src/Forms/quadraticform.jl

@@ -0,0 +1,70 @@
+struct QuadraticForm{V,KetMap,BraMap,OperatorMap,VectorSpaceMap,LinkSpaceMap} <:
+       AbstractForm{V,KetMap,BraMap,OperatorMap,VectorSpaceMap}
+  tn::ITensorNetwork{V}
+  bra_vs_map::BraMap
+  ket_vs_map::KetMap
+  operator_vs_map::OperatorMap
+  dual_space_map::VectorSpaceMap
+  link_space_map::LinkSpaceMap
+end
+
+#Needed for implementation
+bra_map(qf::QuadraticForm) = qf.bra_vs_map
+ket_map(qf::QuadraticForm) = qf.ket_vs_map
+operator_map(qf::QuadraticForm) = qf.operator_vs_map
+space_map(qf::QuadraticForm) = qf.dual_space_map
+tensornetwork(qf::QuadraticForm) = qf.tn
+data_graph_type(::Type{<:QuadraticForm}) = data_graph_type(tensornetwork(qf))
+data_graph(qf::QuadraticForm) = data_graph(tensornetwork(qf))
+
+function QuadraticForm(
+  ket::ITensorNetwork, operator::ITensorNetwork; link_space_map=default_index_map, kwargs...
+)
+  bra = link_space_map(dag(ket); sites=[])
+  return Form(bra, operator, ket; form_type="Quadratic", link_space_map, kwargs...)
+end
+
+function QuadraticForm(
+  ket::ITensorNetwork;
+  dual_space_map=default_index_map,
+  tno_constructor=default_tno_constructor,
+  kwargs...,
+)
+  s = siteinds(ket)
+  operator_space = union_all_inds(s, dual_space_map(s; links=[]))
+  operator = tno_constructor(operator_space)
+  return QuadraticForm(ket, operator; kwargs...)
+end
+
+function bra_ket_vertices(qf::QuadraticForm, state_vertices::Vector)
+  ket_vertices = [ket_map(qf)[sv] for sv in state_vertices]
+  bra_vertices = [bra_map(qf)[sv] for sv in state_vertices]
+  return bra_vertices, ket_vertices
+end
+
+function gradient(
+  qf::QuadraticForm, state_vertices::Vector; alg="Exact", (cache!)=nothing, kwargs...
+)
+  return gradient(Algorithm(alg), qf, state_vertices; (cache!), kwargs...)
+end
+
+function gradient(
+  ::Algorithm"Exact",
+  qf::QuadraticForm,
+  state_vertices::Vector;
+  (cache!)=nothing,
+  contractor_kwargs...,
+)
+  qf_bra_vertices, qf_ket_vertices = bra_ket_vertices(qf, state_vertices)
+  tn_vertices = setdiff(vertices(qf), vcat(qf_bra_vertices, qf_ket_vertices))
+  reduced_tn, _ = induced_subgraph(tensornetwork(qf), tn_vertices)
+  return contract(reduced_tn, contractor_kwargs...)
+end
+
+function update_qf!(qf::QuadraticForm, state_vertex, state::ITensor)
+  setindex_preserve_graph!(tensornetwork(qf)[ket_map(state_vertex)], state, state_vertex)
+  #Change to be consistent with dual space map and link space map. How can we do this elegantly??
+  return setindex_preserve_graph!(
+    tensornetwork(qf)[bra_map(state_vertex)], prime(dag(state)), state_vertex
+  )
+end

src/ITensorNetworks.jl

@@ -98,6 +98,10 @@ include("renameitensornetwork.jl")
 include("boundarymps.jl")
 include(joinpath("beliefpropagation", "beliefpropagation.jl"))
 include(joinpath("beliefpropagation", "beliefpropagation_schedule.jl"))
+include(joinpath("Forms", "abstractbilinearform.jl"))
+include(joinpath("Forms", "bilinearform.jl"))
+include(joinpath("Forms", "quadraticform.jl"))
+include(joinpath("Forms", "construct_form.jl"))
 include("contraction_tree_to_graph.jl")
 include("gauging.jl")
 include("utils.jl")

src/exports.jl

@@ -71,6 +71,8 @@ export AbstractITensorNetwork,
   mps,
   ortho_center,
   set_ortho_center,
+  QuadraticForm,
+  BilinearForm,
   TreeTensorNetwork,
   TTN,
   random_ttn,

test/test_forms.jl

@@ -0,0 +1,31 @@
+using ITensors
+using Graphs
+using NamedGraphs
+using ITensorNetworks
+using ITensorNetworks: delta_network, bra_ket_vertices, gradient
+using Test
+using Random
+
+@testset "Forms" begin
+  g = named_grid((1, 4))
+  s = siteinds("S=1/2", g)
+  s_operator = union_all_inds(s, prime(s; links=[]))
+  χ = 2
+  Random.seed!(1234)
+  ψket = randomITensorNetwork(s; link_space=χ)
+  ψbra = randomITensorNetwork(s; link_space=χ)
+  A = delta_network(s_operator)
+
+  blf = BilinearForm(ψbra, A, ψket)
+  @test nv(blf) == nv(ψket) + nv(ψbra) + nv(A)
+  @test isempty(externalinds(blf))
+
+  qf = QuadraticForm(ψbra, A)
+  @test nv(qf) == nv(ψket) + nv(ψbra) + nv(A)
+  @test isempty(externalinds(qf))
+  v = (1, 1)
+
+  qf_gradient_v = gradient(qf, [v]; alg="Exact")
+  z_qf = ITensors.contract([qf_gradient_v, dag(prime(ψbra[v])), ψbra[v]])[]
+  @test z_qf ≈ ITensors.contract(qf)[]
+end