From a6250acf6118b9ed6b0d4e983c7d00878faee7c9 Mon Sep 17 00:00:00 2001 From: Joseph Tindall <51231103+JoeyT1994@users.noreply.github.com> Date: Fri, 26 Jan 2024 10:50:59 -0500 Subject: [PATCH] Use `NamedGraphs.PartitionedGraph` type to rewrite BP code (#126) --- examples/apply/apply_bp/apply_bp.jl | 170 --------- examples/apply/apply_bp/apply_bp_run.jl | 72 ---- examples/belief_propagation/bpexample.jl | 113 ------ examples/belief_propagation/bpsequences.jl | 81 ----- examples/belief_propagation/sqrt_bp.jl | 78 ---- examples/boundary.jl | 5 +- examples/distances.jl | 2 +- .../dynamics/heavy_hex_ising_real_tebd.jl | 134 ------- examples/gauging/gauging_itns.jl | 144 -------- examples/group_partition.jl | 17 - examples/partition/kahypar_vs_metis.jl | 44 --- examples/partition/partitioning.jl | 31 -- src/ITensorNetworks.jl | 11 +- src/apply.jl | 6 +- .../approx_itensornetwork.jl | 2 +- .../binary_tree_partition.jl | 8 +- src/approx_itensornetwork/partition.jl | 99 +++++ src/beliefpropagation/beliefpropagation.jl | 199 ++++------ .../beliefpropagation_schedule.jl | 3 + .../sqrt_beliefpropagation.jl | 149 -------- src/boundarymps.jl | 2 +- src/contract.jl | 27 ++ src/exports.jl | 1 + src/gauging.jl | 87 ++--- src/itensornetwork.jl | 9 +- src/partition.jl | 341 ------------------ src/requires/kahypar.jl | 33 -- src/requires/metis.jl | 33 -- test/test_apply.jl | 31 +- test/test_belief_propagation.jl | 111 ++---- test/test_binary_tree_partition.jl | 11 +- test/test_contract_deltas.jl | 2 +- test/test_examples/test_apply_bp.jl | 28 -- test/test_examples/test_examples.jl | 10 +- test/test_examples/test_sqrt_bp.jl | 10 - test/test_gauging.jl | 28 +- 36 files changed, 327 insertions(+), 1805 deletions(-) delete mode 100644 examples/apply/apply_bp/apply_bp.jl delete mode 100644 examples/apply/apply_bp/apply_bp_run.jl delete mode 100644 examples/belief_propagation/bpexample.jl delete mode 100644 examples/belief_propagation/bpsequences.jl delete mode 100644 examples/belief_propagation/sqrt_bp.jl delete mode 100644 examples/dynamics/heavy_hex_ising_real_tebd.jl delete mode 100644 examples/gauging/gauging_itns.jl delete mode 100644 examples/group_partition.jl delete mode 100644 examples/partition/kahypar_vs_metis.jl delete mode 100644 examples/partition/partitioning.jl rename src/{ => approx_itensornetwork}/binary_tree_partition.jl (96%) create mode 100644 src/approx_itensornetwork/partition.jl delete mode 100644 src/beliefpropagation/sqrt_beliefpropagation.jl delete mode 100644 src/partition.jl delete mode 100644 src/requires/kahypar.jl delete mode 100644 src/requires/metis.jl delete mode 100644 test/test_examples/test_apply_bp.jl delete mode 100644 test/test_examples/test_sqrt_bp.jl diff --git a/examples/apply/apply_bp/apply_bp.jl b/examples/apply/apply_bp/apply_bp.jl deleted file mode 100644 index 3b02b923..00000000 --- a/examples/apply/apply_bp/apply_bp.jl +++ /dev/null @@ -1,170 +0,0 @@ -using ITensorNetworks -using ITensorNetworks: - approx_network_region, - belief_propagation, - get_environment, - contract_inner, - find_subgraph, - message_tensors, - neighbor_vertices, - nested_graph_leaf_vertices, - symmetric_gauge, - vidal_gauge, - vidal_to_symmetric_gauge, - norm_network -using Test -using Compat -using Dictionaries -using ITensors -using Metis -using NamedGraphs -using Observers -using Random -using LinearAlgebra -using SplitApplyCombine -using OMEinsumContractionOrders - -function expect_bp(opname, v, ψ, mts) - s = siteinds(ψ) - ψψ = norm_network(ψ) - numerator_network = approx_network_region( - ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork(ITensor[apply(op(opname, s[v]), ψ[v])]) - ) - denominator_network = approx_network_region(ψψ, mts, [(v, 1)]) - return contract(numerator_network)[] / contract(denominator_network)[] -end - -function vertex_array(ψ, v, v⃗ⱼ) - indsᵥ = unioninds((linkinds(ψ, v => vⱼ) for vⱼ in v⃗ⱼ)...) - indsᵥ = unioninds(siteinds(ψ, v), indsᵥ) - ψᵥ = ψ[v] - ψᵥ /= norm(ψᵥ) - return array(permute(ψᵥ, indsᵥ)) -end - -function simple_update_bp( - os, - ψ::ITensorNetwork; - maxdim, - variational_optimization_only=false, - regauge=false, - reduced=true, -) - println("Simple update, BP") - ψψ = norm_network(ψ) - mts = message_tensors(partition(ψψ, group(v -> v[1], vertices(ψψ)))) - mts = belief_propagation( - ψψ, mts; contract_kwargs=(; alg="exact"), niters=50, target_precision=1e-5 - ) - for layer in eachindex(os) - @show layer - o⃗ = os[layer] - for o in o⃗ - v⃗ = neighbor_vertices(ψ, o) - for e in edges(mts) - @assert order(only(mts[e])) == 2 - @assert order(only(mts[reverse(e)])) == 2 - end - - @assert length(v⃗) == 2 - v1, v2 = v⃗ - - s1 = find_subgraph((v1, 1), mts) - s2 = find_subgraph((v2, 1), mts) - envs = get_environment(ψψ, mts, [(v1, 1), (v1, 2), (v2, 1), (v2, 2)]) - obs = observer() - # TODO: Make a version of `apply` that accepts message tensors, - # and computes the environment and does the message tensor update of the bond internally. - ψ = apply( - o, - ψ; - envs, - (observer!)=obs, - maxdim, - normalize=true, - variational_optimization_only, - nfullupdatesweeps=20, - symmetrize=true, - reduced, - ) - S = only(obs.singular_values) - S /= norm(S) - - # Update message tensor - ψψ = norm_network(ψ) - mts[s1] = ITensorNetwork(dictionary([(v1, 1) => ψψ[v1, 1], (v1, 2) => ψψ[v1, 2]])) - mts[s2] = ITensorNetwork(dictionary([(v2, 1) => ψψ[v2, 1], (v2, 2) => ψψ[v2, 2]])) - mts[s1 => s2] = ITensorNetwork(obs.singular_values) - mts[s2 => s1] = ITensorNetwork(obs.singular_values) - end - if regauge - println("regauge") - mts = belief_propagation( - ψψ, mts; contract_kwargs=(; alg="exact"), niters=50, target_precision=1e-5 - ) - end - end - return ψ, mts -end - -function simple_update_vidal(os, ψ::ITensorNetwork; maxdim, regauge=false) - println("Simple update, Vidal gauge") - ψψ = norm_network(ψ) - mts = message_tensors(partition(ψψ, group(v -> v[1], vertices(ψψ)))) - mts = belief_propagation( - ψψ, mts; contract_kwargs=(; alg="exact"), niters=50, target_precision=1e-5 - ) - ψ, bond_tensors = vidal_gauge(ψ, mts) - for layer in eachindex(os) - @show layer - o⃗ = os[layer] - for o in o⃗ - v⃗ = neighbor_vertices(ψ, o) - ψ, bond_tensors = apply(o, ψ, bond_tensors; maxdim, normalize=true) - end - if regauge - println("regauge") - ψ_symmetric, mts = vidal_to_symmetric_gauge(ψ, bond_tensors) - ψψ = norm_network(ψ_symmetric) - mts = belief_propagation( - ψψ, mts; contract_kwargs=(; alg="exact"), niters=50, target_precision=1e-5 - ) - ψ, bond_tensors = vidal_gauge(ψ_symmetric, mts) - end - end - return ψ, bond_tensors -end - -function main(; - seed=5623, - graph, - opname, - dims, - χ, - nlayers, - variational_optimization_only=false, - regauge=false, - reduced=true, -) - Random.seed!(seed) - n = prod(dims) - g = graph(dims) - s = siteinds("S=1/2", g) - ψ = randomITensorNetwork(s; link_space=χ) - es = edges(g) - os = [ - [op(opname, s[src(e)]..., s[dst(e)]...; eltype=Float64) for e in es] for _ in 1:nlayers - ] - - # BP SU - ψ_bp, mts = simple_update_bp( - os, ψ; maxdim=χ, variational_optimization_only, regauge, reduced - ) - # ψ_bp, mts = vidal_to_symmetric_gauge(vidal_gauge(ψ_bp, mts)...) - - # Vidal SU - ψ_vidal, bond_tensors = simple_update_vidal(os, ψ; maxdim=χ, regauge) - ψ_vidal, mts_vidal = vidal_to_symmetric_gauge(ψ_vidal, bond_tensors) - - return ψ_bp, mts, ψ_vidal, mts_vidal -end diff --git a/examples/apply/apply_bp/apply_bp_run.jl b/examples/apply/apply_bp/apply_bp_run.jl deleted file mode 100644 index 723e70ae..00000000 --- a/examples/apply/apply_bp/apply_bp_run.jl +++ /dev/null @@ -1,72 +0,0 @@ -include("apply_bp.jl") - -# opname = "Id" -opname = "RandomUnitary" - -@show opname - -# graph = named_comb_tree -graph = named_grid - -dims = (6, 6) - -ψ_bp, mts_bp, ψ_vidal, mts_vidal = main(; - seed=1234, - opname, - graph, - dims, - χ=2, - nlayers=2, - variational_optimization_only=false, - regauge=false, - reduced=true, -) - -v = dims .÷ 2 - -sz_bp = @show expect_bp("Sz", v, ψ_bp, mts_bp) -sz_vidal = @show expect_bp("Sz", v, ψ_vidal, mts_vidal) -@show abs(sz_bp - sz_vidal) / abs(sz_vidal) - -# Run BP again -mts_bp = belief_propagation( - norm_network(ψ_bp), - mts_bp; - contract_kwargs=(; alg="exact"), - niters=50, - target_precision=1e-5, -) -mts_vidal = belief_propagation( - norm_network(ψ_vidal), - mts_vidal; - contract_kwargs=(; alg="exact"), - niters=50, - target_precision=1e-5, -) - -sz_bp = @show expect_bp("Sz", v, ψ_bp, mts_bp) -sz_vidal = @show expect_bp("Sz", v, ψ_vidal, mts_vidal) -@show abs(sz_bp - sz_vidal) / abs(sz_vidal) - -ψ_symmetric, _ = symmetric_gauge(ψ_bp) - -v⃗ⱼ = [v .+ (1, 0), v .- (1, 0), v .+ (0, 1), v .- (0, 1)] -ψ_bp_v = vertex_array(ψ_bp, v, v⃗ⱼ) -ψ_vidal_v = vertex_array(ψ_vidal, v, v⃗ⱼ) -ψ_symmetric_v = vertex_array(ψ_symmetric, v, v⃗ⱼ) - -@show norm(abs.(ψ_bp_v) - abs.(ψ_vidal_v)) -@show norm(abs.(ψ_bp_v) - abs.(ψ_symmetric_v)) -@show norm(abs.(ψ_vidal_v) - abs.(ψ_symmetric_v)) - -# inner_ψ_net = inner_network(ψ_bp, ψ_vidal) -# norm_ψ_bp_net = norm_network(ψ_bp) -# norm_ψ_vidal_net = norm_network(ψ_vidal) -# seq = contraction_sequence(inner_ψ_net; alg="sa_bipartite") -# @disable_warn_order begin -# inner_ψ = contract(inner_ψ_net; sequence=seq)[] -# norm_sqr_ψ_bp = contract(norm_ψ_bp_net; sequence=seq)[] -# norm_sqr_ψ_vidal = contract(norm_ψ_vidal_net; sequence=seq)[] -# end -# @show 1 - inner_ψ / (sqrt(norm_sqr_ψ_bp) * sqrt(norm_sqr_ψ_vidal)) -# @show log(inner_ψ) - (log(norm_sqr_ψ_bp) / 2 + log(norm_sqr_ψ_vidal) / 2) diff --git a/examples/belief_propagation/bpexample.jl b/examples/belief_propagation/bpexample.jl deleted file mode 100644 index a0af8818..00000000 --- a/examples/belief_propagation/bpexample.jl +++ /dev/null @@ -1,113 +0,0 @@ -using Compat -using ITensors -using Metis -using ITensorNetworks -using Random -using SplitApplyCombine -using NamedGraphs - -using ITensorNetworks: - belief_propagation, - approx_network_region, - contract_inner, - message_tensors, - nested_graph_leaf_vertices - -function main() - n = 4 - g_dims = (n, n) - g = named_grid(g_dims) - s = siteinds("S=1/2", g) - chi = 2 - - Random.seed!(5467) - - #bra - ψ = randomITensorNetwork(s; link_space=chi) - #bra-ket (but not actually contracted) - ψψ = ψ ⊗ prime(dag(ψ); sites=[]) - - #Site to take expectation value on - v = (1, 1) - - #Now do Simple Belief Propagation to Measure Sz on Site v - mts = message_tensors( - ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ)))) - ) - - mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"), niters=20) - - numerator_network = approx_network_region( - ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork([apply(op("Sz", s[v]), ψ[v])]) - ) - denominator_network = approx_network_region(ψψ, mts, [(v, 1)]) - sz_bp = contract(numerator_network)[] / contract(denominator_network)[] - - println( - "Simple Belief Propagation Gives Sz on Site " * string(v) * " as " * string(sz_bp) - ) - - #Now do General Belief Propagation to Measure Sz on Site v - nsites = 4 - Zp = partition( - partition(ψψ, group(v -> v[1], vertices(ψψ))); nvertices_per_partition=nsites - ) - Zpp = partition(ψψ; subgraph_vertices=nested_graph_leaf_vertices(Zp)) - mts = message_tensors(Zpp) - mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"), niters=20) - numerator_network = approx_network_region( - ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork([apply(op("Sz", s[v]), ψ[v])]) - ) - denominator_network = approx_network_region(ψψ, mts, [(v, 1)]) - sz_bp = contract(numerator_network)[] / contract(denominator_network)[] - - println( - "General Belief Propagation (4-site subgraphs) Gives Sz on Site " * - string(v) * - " as " * - string(sz_bp), - ) - - #Now do General Belief Propagation with Matrix Product State Message Tensors Measure Sz on Site v - ψψ = flatten_networks(ψ, dag(ψ); combine_linkinds=false, map_bra_linkinds=prime) - Oψ = copy(ψ) - Oψ[v] = apply(op("Sz", s[v]), ψ[v]) - ψOψ = flatten_networks(ψ, dag(Oψ); combine_linkinds=false, map_bra_linkinds=prime) - - combiners = linkinds_combiners(ψψ) - ψψ = combine_linkinds(ψψ, combiners) - ψOψ = combine_linkinds(ψOψ, combiners) - - Z = partition(ψψ, group(v -> v[1], vertices(ψψ))) - maxdim = 8 - mts = message_tensors(Z) - - mts = belief_propagation( - ψψ, - mts; - contract_kwargs=(; - alg="density_matrix", - output_structure=path_graph_structure, - maxdim, - contraction_sequence_alg="optimal", - ), - ) - - numerator_network = approx_network_region(ψψ, mts, [v]; verts_tn=ITensorNetwork(ψOψ[v])) - denominator_network = approx_network_region(ψψ, mts, [v]) - sz_bp = contract(numerator_network)[] / contract(denominator_network)[] - - println( - "General Belief Propagation with Column Partitioning and MPS Message Tensors (Max dim 8) Gives Sz on Site " * - string(v) * - " as " * - string(sz_bp), - ) - - #Now do it exactly - sz_exact = contract(ψOψ)[] / contract(ψψ)[] - - return println("The exact value of Sz on Site " * string(v) * " is " * string(sz_exact)) -end - -main() diff --git a/examples/belief_propagation/bpsequences.jl b/examples/belief_propagation/bpsequences.jl deleted file mode 100644 index 456e96c7..00000000 --- a/examples/belief_propagation/bpsequences.jl +++ /dev/null @@ -1,81 +0,0 @@ -using Compat -using ITensors -using Metis -using ITensorNetworks -using Random -using SplitApplyCombine -using Graphs -using NamedGraphs - -using ITensorNetworks: - belief_propagation, - approx_network_region, - contract_inner, - message_tensors, - nested_graph_leaf_vertices, - edge_sequence - -function main() - g_labels = [ - "Comb Tree", - "100 Site Random Regular Graph z = 3", - "6x6 Square Grid", - "4x4 Hexagonal Lattice", - ] - gs = [ - named_comb_tree((6, 6)), - NamedGraph(Graphs.random_regular_graph(100, 3)), - named_grid((6, 6)), - NamedGraphs.hexagonal_lattice_graph(4, 4), - ] - χs = [4, 4, 2, 3] - - for (i, g) in enumerate(gs) - Random.seed!(5467) - g_label = g_labels[i] - χ = χs[i] - s = siteinds("S=1/2", g) - ψ = randomITensorNetwork(s; link_space=χ) - ψψ = ψ ⊗ prime(dag(ψ); sites=[]) - - #Initial message tensors for BP - mts_init = message_tensors( - ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ)))) - ) - - println("\nFirst testing out a $g_label. Random network with bond dim $χ") - - #Now test out various sequences - print("Parallel updates (sequence is irrelevant): ") - belief_propagation( - ψψ, - mts_init; - contract_kwargs=(; alg="exact"), - target_precision=1e-10, - niters=100, - edges=edge_sequence(mts_init; alg="parallel"), - verbose=true, - ) - print("Sequential updates (sequence is default edge list of the message tensors): ") - belief_propagation( - ψψ, - mts_init; - contract_kwargs=(; alg="exact"), - target_precision=1e-10, - niters=100, - edges=[e for e in edges(mts_init)], - verbose=true, - ) - print("Sequential updates (sequence is our custom sequence finder): ") - belief_propagation( - ψψ, - mts_init; - contract_kwargs=(; alg="exact"), - target_precision=1e-10, - niters=100, - verbose=true, - ) - end -end - -main() diff --git a/examples/belief_propagation/sqrt_bp.jl b/examples/belief_propagation/sqrt_bp.jl deleted file mode 100644 index 207fb941..00000000 --- a/examples/belief_propagation/sqrt_bp.jl +++ /dev/null @@ -1,78 +0,0 @@ -using ITensors -using ITensorNetworks -using Random -using SplitApplyCombine - -using ITensorNetworks: - approx_network_region, - belief_propagation, - sqrt_belief_propagation, - ising_network_state, - message_tensors - -function main(; n, niters, network="ising", β=nothing, h=nothing, χ=nothing) - g_dims = (n, n) - @show g_dims - g = named_grid(g_dims) - s = siteinds("S=1/2", g) - - Random.seed!(5467) - - ψ = if network == "ising" - ising_network_state(s, β; h) - elseif network == "random" - randomITensorNetwork(s; link_space=χ) - else - error("Network type $network not supported.") - end - - ψψ = norm_network(ψ) - - # Site to take expectation value on - v = (n ÷ 2, n ÷ 2) - @show v - - #Now do Simple Belief Propagation to Measure Sz on Site v - mts = message_tensors( - ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ)))) - ) - - mts = @time belief_propagation(ψψ, mts; niters, contract_kwargs=(; alg="exact")) - numerator_network = approx_network_region( - ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork([apply(op("Sz", s[v]), ψ[v])]) - ) - denominator_network = approx_network_region(ψψ, mts, [(v, 1)]) - sz_bp = - contract( - numerator_network; sequence=contraction_sequence(numerator_network; alg="optimal") - )[] / contract( - denominator_network; sequence=contraction_sequence(denominator_network; alg="optimal") - )[] - - println( - "Simple Belief Propagation Gives Sz on Site " * string(v) * " as " * string(sz_bp) - ) - - mts_sqrt = message_tensors( - ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ)))) - ) - - mts_sqrt = @time sqrt_belief_propagation(ψ, mts_sqrt; niters) - - numerator_network = approx_network_region( - ψψ, mts_sqrt, [(v, 1)]; verts_tn=ITensorNetwork([apply(op("Sz", s[v]), ψ[v])]) - ) - denominator_network = approx_network_region(ψψ, mts_sqrt, [(v, 1)]) - sz_sqrt_bp = - contract( - numerator_network; sequence=contraction_sequence(numerator_network; alg="optimal") - )[] / contract( - denominator_network; sequence=contraction_sequence(denominator_network; alg="optimal") - )[] - - println( - "Sqrt Belief Propagation Gives Sz on Site " * string(v) * " as " * string(sz_sqrt_bp) - ) - - return (; sz_bp, sz_sqrt_bp) -end diff --git a/examples/boundary.jl b/examples/boundary.jl index c10aecbf..7c43712a 100644 --- a/examples/boundary.jl +++ b/examples/boundary.jl @@ -8,8 +8,9 @@ tn = ITensorNetwork(named_grid((6, 3)); link_space=4) @visualize tn -g = subgraph_vertices(tn; nvertices_per_partition=2) -sub_vs_1, sub_vs_2 = g[1], g[2] +ptn = PartitionedGraph(tn; nvertices_per_partition=2) +sub_vs_1 = vertices(ptn, PartitionVertex(1)) +sub_vs_2 = vertices(ptn, PartitionVertex(2)) @show (1, 1) ∈ sub_vs_1 @show (6, 3) ∈ sub_vs_2 diff --git a/examples/distances.jl b/examples/distances.jl index 38df7296..7767651f 100644 --- a/examples/distances.jl +++ b/examples/distances.jl @@ -14,4 +14,4 @@ t = dijkstra_tree(ψ, only(center(ψ))) @show a_star(ψ, (2, 1), (2, 5)) @show mincut_partitions(ψ) @show mincut_partitions(ψ, (1, 1), (3, 5)) -@show subgraph_vertices(ψ; npartitions=2) +@show partitioned_vertices(ψ; npartitions=2) diff --git a/examples/dynamics/heavy_hex_ising_real_tebd.jl b/examples/dynamics/heavy_hex_ising_real_tebd.jl deleted file mode 100644 index e8dcc863..00000000 --- a/examples/dynamics/heavy_hex_ising_real_tebd.jl +++ /dev/null @@ -1,134 +0,0 @@ -using ITensors -using ITensorNetworks -using NamedGraphs -using NamedGraphs: rem_edges!, add_edge!, decorate_graph_edges, hexagonal_lattice_graph -using Graphs - -using ITensorNetworks: - initialize_bond_tensors, - vidal_itn_canonicalness, - vidal_to_symmetric_gauge, - vidal_gauge, - approx_network_region, - optimal_contraction_sequence, - contract, - symmetric_to_vidal_gauge, - norm_network - -function heavy_hex_lattice_graph(n::Int64, m::Int64) - g = hexagonal_lattice_graph(n, m) - g = decorate_graph_edges(g) - return g -end - -function ibm_processor_graph(n::Int64, m::Int64) - g = heavy_hex_lattice_graph(n, m) - dims = maximum(vertices(hexagonal_lattice_graph(n, m))) - v1, v2 = (dims[1], 1), (1, dims[2]) - add_vertices!(g, [v1, v2]) - add_edge!(g, v1 => v1 .- (1, 0)) - add_edge!(g, v2 => v2 .+ (1, 0)) - - return g -end - -eagle_processor_graph() = ibm_processor_graph(3, 6) -hummingbird_processor_graph() = ibm_processor_graph(2, 4) -osprey_processor_graph() = ibm_processor_graph(6, 12) - -"""Take the expectation value of o on an ITN using belief propagation""" -function expect_state_SBP( - o::ITensor, ψ::AbstractITensorNetwork, ψψ::AbstractITensorNetwork, mts::DataGraph -) - Oψ = apply(o, ψ; cutoff=1e-16) - ψ = copy(ψ) - s = siteinds(ψ) - vs = vertices(s)[findall(i -> (length(commoninds(s[i], inds(o))) != 0), vertices(s))] - vs_braket = [(v, 1) for v in vs] - - numerator_network = approx_network_region( - ψψ, mts, vs_braket; verts_tn=ITensorNetwork(ITensor[Oψ[v] for v in vs]) - ) - denominator_network = approx_network_region(ψψ, mts, vs_braket) - num_seq = contraction_sequence(numerator_network; alg="optimal") - den_seq = contraction_sequence(numerator_network; alg="optimal") - return contract(numerator_network; sequence=num_seq)[] / - contract(denominator_network; sequence=den_seq)[] -end - -function main(θh::Float64, no_trotter_steps::Int64; apply_kwargs...) - - #Build the graph - g = eagle_processor_graph() - - #Do this if measuring a Z based expectation value (i.e. ignore ZZ_gates in final layer as they are irrelevant) - shortened_final_layer = true - - s = siteinds("S=1/2", g) - - #Gauging parameters - target_c = 1e-3 - gauge_freq = 1 - - #State initialisation - ψ = ITensorNetwork(s, v -> "↑") - bond_tensors = initialize_bond_tensors(ψ) - - #Build gates - HX, HZZ = ising(g; h=θh / 2, J1=0), ising(g; h=0, J1=-pi / 4) - RX, RZZ = exp(-im * HX; alg=Trotter{1}()), exp(-im * HZZ; alg=Trotter{1}()) - RX_gates, RZZ_gates = Vector{ITensor}(RX, s), Vector{ITensor}(RZZ, s) - - for i in 1:no_trotter_steps - println("On Trotter Step $i") - gates = if (shortened_final_layer && i == no_trotter_steps) - RX_gates - else - vcat(RX_gates, RZZ_gates) - end - - for gate in gates - ψ, bond_tensors = apply(gate, ψ, bond_tensors; normalize=true, apply_kwargs...) - end - - cur_C = vidal_itn_canonicalness(ψ, bond_tensors) - if ((i + 1) % gauge_freq) == 0 && (cur_C >= target_c) - println("Too far from the Vidal gauge. Regauging the state!") - ψ, _ = vidal_to_symmetric_gauge(ψ, bond_tensors) - ψ, bond_tensors = vidal_gauge( - ψ; target_canonicalness=target_c, niters=20, cutoff=1e-14 - ) - end - end - - #Calculate on-site magnetisations - ψ, mts = vidal_to_symmetric_gauge(ψ, bond_tensors) - ψψ = norm_network(ψ) - mag_dict = Dict( - zip( - [v for v in vertices(ψ)], - [expect_state_SBP(op("Z", s[v]), ψ, ψψ, mts) for v in vertices(ψ)], - ), - ) - - return mag_dict -end - -ibm_processor_graph(3, 6) - -θh = pi / 4 -no_trotter_steps = 5 -χ = 32 -apply_kwargs = (; maxdim=χ, cutoff=1e-12) - -println( - "Simulating $no_trotter_steps Steps of the Kicked Ising Model on the Heavy Hex Layout." -) -println("θh is $θh. Maxdim is $χ.") - -mags = main(θh, no_trotter_steps; apply_kwargs...) - -println( - "After $no_trotter_steps steps. Average magnetisation is " * - string(sum(values(mags)) / length(values(mags))), -) diff --git a/examples/gauging/gauging_itns.jl b/examples/gauging/gauging_itns.jl deleted file mode 100644 index 7463d54a..00000000 --- a/examples/gauging/gauging_itns.jl +++ /dev/null @@ -1,144 +0,0 @@ -using Compat -using ITensors -using Metis -using ITensorNetworks -using Random -using SplitApplyCombine - -using ITensorNetworks: - message_tensors, - nested_graph_leaf_vertices, - belief_propagation_iteration, - belief_propagation, - find_subgraph, - vidal_gauge, - symmetric_gauge, - vidal_itn_canonicalness, - vidal_to_symmetric_gauge, - initialize_bond_tensors, - vidal_itn_isometries, - norm_network, - edge_sequence - -using NamedGraphs -using NamedGraphs: add_edges!, rem_vertex!, hexagonal_lattice_graph -using Graphs - -"""Eager Gauging""" -function eager_gauging(ψ::ITensorNetwork, bond_tensors::DataGraph, mts::DataGraph) - isometries = vidal_itn_isometries(ψ, bond_tensors) - - ψ = copy(ψ) - mts = copy(mts) - for e in edges(ψ) - s1, s2 = find_subgraph((src(e), 1), mts), find_subgraph((dst(e), 1), mts) - normalize!(isometries[e]) - normalize!(isometries[reverse(e)]) - mts[s1 => s2], mts[s2 => s1] = ITensorNetwork(isometries[e]), - ITensorNetwork(isometries[reverse(e)]) - end - - ψ, bond_tensors = vidal_gauge(ψ, mts, bond_tensors) - - return ψ, bond_tensors, mts -end - -"""Bring an ITN into the Vidal gauge, various methods possible. Result is timed""" -function benchmark_state_gauging( - ψ::ITensorNetwork; - mode="belief_propagation", - no_iterations=50, - BP_update_order::String="sequential", -) - s = siteinds(ψ) - - C = zeros((no_iterations)) - times_iters = zeros((no_iterations)) - times_gauging = zeros((no_iterations)) - - ψψ = norm_network(ψ) - ψinit = copy(ψ) - vertex_groups = nested_graph_leaf_vertices(partition(ψψ, group(v -> v[1], vertices(ψψ)))) - mts = message_tensors(partition(ψψ; subgraph_vertices=vertex_groups)) - bond_tensors = initialize_bond_tensors(ψ) - for e in edges(mts) - mts[e] = ITensorNetwork(dense(delta(inds(ITensors.contract(ITensor(mts[e])))))) - end - - for i in 1:no_iterations - println("On Iteration " * string(i)) - - if mode == "belief_propagation" - if BP_update_order != "parallel" - times_iters[i] = @elapsed mts, _ = belief_propagation_iteration( - ψψ, mts; contract_kwargs=(; alg="exact") - ) - else - times_iters[i] = @elapsed mts, _ = belief_propagation_iteration( - ψψ, mts; contract_kwargs=(; alg="exact"), edges=edge_sequence(mts; alg="parallel") - ) - end - - times_gauging[i] = @elapsed ψ, bond_tensors = vidal_gauge(ψinit, mts) - elseif mode == "eager" - times_iters[i] = @elapsed ψ, bond_tensors, mts = eager_gauging(ψ, bond_tensors, mts) - else - times_iters[i] = @elapsed begin - for e in edges(ψ) - ψ, bond_tensors = apply(e, ψ, bond_tensors; normalize=true, cutoff=1e-16) - end - end - end - - C[i] = vidal_itn_canonicalness(ψ, bond_tensors) - end - @show times_iters, time - simulation_times = cumsum(times_iters) + times_gauging - - return simulation_times, C -end - -L, χ = 10, 10 -g = named_grid((L, L)) -s = siteinds("S=1/2", g) -ψ = randomITensorNetwork(s; link_space=χ) -no_iterations = 30 - -BPG_simulation_times, BPG_Cs = benchmark_state_gauging( - ψ; no_iterations, BP_update_order="parallel" -) -BPG_sequential_simulation_times, BPG_sequential_Cs = benchmark_state_gauging( - ψ; no_iterations -) -Eager_simulation_times, Eager_Cs = benchmark_state_gauging(ψ; mode="eager", no_iterations) -SU_simulation_times, SU_Cs = benchmark_state_gauging(ψ; mode="SU", no_iterations) - -epsilon = 1e-10 - -println( - "Time for BPG (with parallel updates) to reach C < epsilon was " * - string(BPG_simulation_times[findfirst(x -> x < 0, BPG_Cs .- epsilon)]) * - " seconds. No iters was " * - string(findfirst(x -> x < 0, BPG_Cs .- epsilon)), -) -println( - "Time for BPG (with sequential updates) to reach C < epsilon was " * - string( - BPG_sequential_simulation_times[findfirst(x -> x < 0, BPG_sequential_Cs .- epsilon)] - ) * - " seconds. No iters was " * - string(findfirst(x -> x < 0, BPG_sequential_Cs .- epsilon)), -) - -println( - "Time for Eager Gauging to reach C < epsilon was " * - string(Eager_simulation_times[findfirst(x -> x < 0, Eager_Cs .- epsilon)]) * - " seconds. No iters was " * - string(findfirst(x -> x < 0, Eager_Cs .- epsilon)), -) -println( - "Time for SU Gauging (with sequential updates) to reach C < epsilon was " * - string(SU_simulation_times[findfirst(x -> x < 0, SU_Cs .- epsilon)]) * - " seconds. No iters was " * - string(findfirst(x -> x < 0, SU_Cs .- epsilon)), -) diff --git a/examples/group_partition.jl b/examples/group_partition.jl deleted file mode 100644 index 1a3c12ae..00000000 --- a/examples/group_partition.jl +++ /dev/null @@ -1,17 +0,0 @@ -using ITensors -using Graphs -using NamedGraphs -using ITensorNetworks -using SplitApplyCombine -using Metis - -s = siteinds("S=1/2", named_grid(8)) -tn = ITensorNetwork(s; link_space=2) -Z = prime(tn; sites=[]) ⊗ tn -vertex_groups = group(v -> v[1], vertices(Z)) -# Create two layers of partitioning -Z_p = partition(partition(Z, vertex_groups); nvertices_per_partition=2) -# Flatten the partitioned partitions -Z_verts = [ - reduce(vcat, (vertices(Z_p[vp][v]) for v in vertices(Z_p[vp]))) for vp in vertices(Z_p) -] diff --git a/examples/partition/kahypar_vs_metis.jl b/examples/partition/kahypar_vs_metis.jl deleted file mode 100644 index 36301a5b..00000000 --- a/examples/partition/kahypar_vs_metis.jl +++ /dev/null @@ -1,44 +0,0 @@ -using Graphs -using KaHyPar -using Metis -using ITensorNetworks - -g = grid((16,)) -npartitions = 4 - -kahypar_partitions = subgraph_vertices(g; npartitions, backend="KaHyPar") -metis_partitions = subgraph_vertices(g; npartitions, backend="Metis") -@show kahypar_partitions, length(kahypar_partitions) -@show metis_partitions, length(metis_partitions) - -g_parts = partition(g; npartitions) -@show nv(g_parts) == 4 -@show nv(g_parts[1]) == 4 -@show nv(g_parts[2]) == 4 -@show nv(g_parts[3]) == 4 -@show nv(g_parts[4]) == 4 -@show issetequal(metis_partitions[2], vertices(g_parts[2])) - -using ITensorNetworks -tn = ITensorNetwork(named_grid((4, 2)); link_space=3); - -# subgraph_vertices -tn_sv = subgraph_vertices(tn; npartitions=2) # Same as `partition_vertices(tn; nvertices_per_partition=4)` - -# partition_vertices -tn_pv = partition_vertices(tn; npartitions=2); -typeof(tn_pv) -tn_pv[1] -edges(tn_pv) -tn_pv[1 => 2] - -# subgraphs -tn_sg = subgraphs(tn; npartitions=2); -typeof(tn_sg) -tn_sg[1] - -# partition -tn_pg = partition(tn; npartitions=2); -typeof(tn_pg) -tn_pg[1] -tn_pg[1 => 2][:edges] diff --git a/examples/partition/partitioning.jl b/examples/partition/partitioning.jl deleted file mode 100644 index 0667d02b..00000000 --- a/examples/partition/partitioning.jl +++ /dev/null @@ -1,31 +0,0 @@ -using Graphs -using ITensors -using ITensorNetworks -using ITensorUnicodePlots -using KaHyPar -using Suppressor - -system_dims = (4, 4) -g = hypercubic_lattice_graph(system_dims) - -s = siteinds("S=1/2", g) - -χ = 5 -ψ = ITensorNetwork(s; link_space=χ) - -@visualize ψ edge_labels = (; plevs=true) - -ψ′ = ψ' -@visualize ψ′ edge_labels = (; plevs=true) - -v = (2, 2) -neighbor_edges = [v => nv for nv in neighbors(ψ, v)] -@show siteinds(ψ, v) -@show [e => linkinds(ψ, e) for e in neighbor_edges] - -npartitions = 4 -partitions = partition(ψ; npartitions) - -@show partitions - -nothing diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 9d5f14d2..24fea311 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -70,7 +70,6 @@ include("observers.jl") include("visualize.jl") include("graphs.jl") include("itensors.jl") -include("partition.jl") include("lattices.jl") include("abstractindsnetwork.jl") include("indextags.jl") @@ -86,11 +85,12 @@ include("tebd.jl") include("itensornetwork.jl") include("mincut.jl") include("contract_deltas.jl") -include("binary_tree_partition.jl") include(joinpath("approx_itensornetwork", "utils.jl")) include(joinpath("approx_itensornetwork", "density_matrix.jl")) include(joinpath("approx_itensornetwork", "ttn_svd.jl")) include(joinpath("approx_itensornetwork", "approx_itensornetwork.jl")) +include(joinpath("approx_itensornetwork", "partition.jl")) +include(joinpath("approx_itensornetwork", "binary_tree_partition.jl")) include("contract.jl") include("utility.jl") include("specialitensornetworks.jl") @@ -98,7 +98,6 @@ include("renameitensornetwork.jl") include("boundarymps.jl") include(joinpath("beliefpropagation", "beliefpropagation.jl")) include(joinpath("beliefpropagation", "beliefpropagation_schedule.jl")) -include(joinpath("beliefpropagation", "sqrt_beliefpropagation.jl")) include("contraction_tree_to_graph.jl") include("gauging.jl") include("utils.jl") @@ -131,12 +130,6 @@ include(joinpath("treetensornetworks", "solvers", "tree_sweeping.jl")) include("exports.jl") function __init__() - @require KaHyPar = "2a6221f6-aa48-11e9-3542-2d9e0ef01880" include( - joinpath("requires", "kahypar.jl") - ) - @require Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" include( - joinpath("requires", "metis.jl") - ) @require OMEinsumContractionOrders = "6f22d1fd-8eed-4bb7-9776-e7d684900715" include( joinpath("requires", "omeinsumcontractionorders.jl") ) diff --git a/src/apply.jl b/src/apply.jl index 676f5867..c7745973 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -297,7 +297,7 @@ end #In the future we will try to unify this into apply() above but currently leave it mostly as a separate function """Apply() function for an ITN in the Vidal Gauge. Hence the bond tensors are required. Gate does not necessarily need to be passed. Can supply an edge to do an identity update instead. Uses Simple Update procedure assuming gate is two-site""" -function ITensors.apply( +function vidal_apply( o::Union{ITensor,NamedEdge}, ψ::AbstractITensorNetwork, bond_tensors::DataGraph; @@ -336,8 +336,8 @@ function ITensors.apply( ind_to_replace = commonind(V, S) ind_to_replace_with = commonind(U, S) - replaceind!(S, ind_to_replace, ind_to_replace_with') - replaceind!(V, ind_to_replace, ind_to_replace_with) + S = replaceind(S, ind_to_replace => ind_to_replace_with') + V = replaceind(V, ind_to_replace => ind_to_replace_with) ψv1, bond_tensors[e], ψv2 = U * Qᵥ₁, S, V * Qᵥ₂ diff --git a/src/approx_itensornetwork/approx_itensornetwork.jl b/src/approx_itensornetwork/approx_itensornetwork.jl index 5041a8a5..c1cbfb20 100644 --- a/src/approx_itensornetwork/approx_itensornetwork.jl +++ b/src/approx_itensornetwork/approx_itensornetwork.jl @@ -69,7 +69,7 @@ function approx_itensornetwork( contraction_sequence_alg="optimal", contraction_sequence_kwargs=(;), ) - par = partition(tn, inds_btree; alg="mincut_recursive_bisection") + par = _partition(tn, inds_btree; alg="mincut_recursive_bisection") output_tn, log_root_norm = approx_itensornetwork( alg, par; diff --git a/src/binary_tree_partition.jl b/src/approx_itensornetwork/binary_tree_partition.jl similarity index 96% rename from src/binary_tree_partition.jl rename to src/approx_itensornetwork/binary_tree_partition.jl index cff0c31e..4a7b5846 100644 --- a/src/binary_tree_partition.jl +++ b/src/approx_itensornetwork/binary_tree_partition.jl @@ -59,7 +59,7 @@ Note: for a given binary tree with n indices, the output partition will contain Note: name of vertices in the output partition are the same as the name of vertices in `inds_btree`. """ -function partition( +function _partition( ::Algorithm"mincut_recursive_bisection", tn::ITensorNetwork, inds_btree::DataGraph ) @assert _is_rooted_directed_binary_tree(inds_btree) @@ -115,7 +115,7 @@ function partition( end tn_deltas = ITensorNetwork(vcat(output_deltas_vector...)) tn_scalars = ITensorNetwork(vcat(scalars_vector...)) - par = partition(disjoint_union(out_tn, tn_deltas, tn_scalars), subgraph_vs) + par = _partition(disjoint_union(out_tn, tn_deltas, tn_scalars), subgraph_vs) @assert is_tree(par) name_map = Dict() for (i, v) in enumerate(pre_order_dfs_vertices(inds_btree, root)) @@ -124,6 +124,6 @@ function partition( return rename_vertices(par, name_map) end -function partition(tn::ITensorNetwork, inds_btree::DataGraph; alg::String) - return partition(Algorithm(alg), tn, inds_btree) +function _partition(tn::ITensorNetwork, inds_btree::DataGraph; alg::String) + return _partition(Algorithm(alg), tn, inds_btree) end diff --git a/src/approx_itensornetwork/partition.jl b/src/approx_itensornetwork/partition.jl new file mode 100644 index 00000000..9f89d063 --- /dev/null +++ b/src/approx_itensornetwork/partition.jl @@ -0,0 +1,99 @@ +function _partition(g::AbstractGraph, subgraph_vertices) + partitioned_graph = DataGraph( + NamedGraph(eachindex(subgraph_vertices)), + map(vs -> subgraph(g, vs), Dictionary(subgraph_vertices)), + ) + for e in edges(g) + s1 = findfirst_on_vertices(subgraph -> src(e) ∈ vertices(subgraph), partitioned_graph) + s2 = findfirst_on_vertices(subgraph -> dst(e) ∈ vertices(subgraph), partitioned_graph) + if (!has_edge(partitioned_graph, s1, s2) && s1 ≠ s2) + add_edge!(partitioned_graph, s1, s2) + partitioned_graph[s1 => s2] = Dictionary( + [:edges, :edge_data], + [Vector{edgetype(g)}(), Dictionary{edgetype(g),edge_data_type(g)}()], + ) + end + if has_edge(partitioned_graph, s1, s2) + push!(partitioned_graph[s1 => s2][:edges], e) + if isassigned(g, e) + set!(partitioned_graph[s1 => s2][:edge_data], e, g[e]) + end + end + end + return partitioned_graph +end + +""" +Find all vertices `v` such that `f(graph[v]) == true` +""" +function findall_on_vertices(f::Function, graph::AbstractDataGraph) + return findall(f, vertex_data(graph)) +end + +""" +Find the vertex `v` such that `f(graph[v]) == true` +""" +function findfirst_on_vertices(f::Function, graph::AbstractDataGraph) + return findfirst(f, vertex_data(graph)) +end + +""" +Find all edges `e` such that `f(graph[e]) == true` +""" +function findall_on_edges(f::Function, graph::AbstractDataGraph) + return findall(f, edge_data(graph)) +end + +""" +Find the edge `e` such that `f(graph[e]) == true` +""" +function findfirst_on_edges(f::Function, graph::AbstractDataGraph) + return findfirst(f, edge_data(graph)) +end + +# function subgraphs(g::AbstractSimpleGraph, subgraph_vertices) +# return subgraphs(NamedGraph(g), subgraph_vertices) +# end + +# """ +# subgraphs(g::AbstractGraph, subgraph_vertices) + +# Return a collection of subgraphs of `g` defined by the subgraph +# vertices `subgraph_vertices`. +# """ +# function subgraphs(g::AbstractGraph, subgraph_vertices) +# return map(vs -> subgraph(g, vs), subgraph_vertices) +# end + +# """ +# subgraphs(g::AbstractGraph; npartitions::Integer, kwargs...) + +# Given a graph `g`, partition `g` into `npartitions` partitions +# or into partitions with `nvertices_per_partition` vertices per partition, +# returning a list of subgraphs. +# Try to keep all subgraphs the same size and minimise edges cut between them. +# A graph partitioning backend such as Metis or KaHyPar needs to be installed for this function to work. +# """ +# function subgraphs( +# g::AbstractGraph; npartitions=nothing, nvertices_per_partition=nothing, kwargs... +# ) +# return subgraphs(g, subgraph_vertices(g; npartitions, nvertices_per_partition, kwargs...)) +# end + +""" + TODO: do we want to make it a public function? +""" +function _noncommoninds(partition::DataGraph) + networks = [Vector{ITensor}(partition[v]) for v in vertices(partition)] + network = vcat(networks...) + return noncommoninds(network...) +end + +# Util functions for partition +function _commoninds(partition::DataGraph) + networks = [Vector{ITensor}(partition[v]) for v in vertices(partition)] + network = vcat(networks...) + outinds = noncommoninds(network...) + allinds = mapreduce(t -> [i for i in inds(t)], vcat, network) + return Vector(setdiff(allinds, outinds)) +end diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl index 8a1f17b4..3b16e3b3 100644 --- a/src/beliefpropagation/beliefpropagation.jl +++ b/src/beliefpropagation/beliefpropagation.jl @@ -1,104 +1,61 @@ -function message_tensors( - tn::ITensorNetwork; - nvertices_per_partition=nothing, - npartitions=nothing, - subgraph_vertices=nothing, - kwargs..., -) - return message_tensors( - partition(tn; nvertices_per_partition, npartitions, subgraph_vertices); kwargs... - ) -end - -function message_tensors_skeleton(subgraphs::DataGraph) - mts = DataGraph{vertextype(subgraphs),vertex_data_type(subgraphs),ITensorNetwork}( - directed_graph(underlying_graph(subgraphs)) - ) - for v in vertices(mts) - mts[v] = subgraphs[v] - end - - return mts +default_mt_constructor(inds_e) = ITensor[denseblocks(delta(inds_e))] +default_bp_cache(ptn::PartitionedGraph) = Dict() +function default_contractor(contract_list::Vector{ITensor}; kwargs...) + return contract_exact(contract_list; kwargs...) end +default_contractor_kwargs() = (; normalize=true, contraction_sequence_alg="optimal") -function message_tensors( - subgraphs::DataGraph; - itensor_constructor=inds_e -> ITensor[dense(delta(i)) for i in inds_e], +function message_tensor( + ptn::PartitionedGraph, edge::PartitionEdge; mt_constructor=default_mt_constructor ) - mts = message_tensors_skeleton(subgraphs) - for e in edges(subgraphs) - inds_e = commoninds(subgraphs[src(e)], subgraphs[dst(e)]) - itensors = itensor_constructor(inds_e) - mts[e] = ITensorNetwork(itensors) - mts[reverse(e)] = dag(mts[e]) - end - return mts + src_e_itn = subgraph(ptn, src(edge)) + dst_e_itn = subgraph(ptn, dst(edge)) + inds_e = commoninds(src_e_itn, dst_e_itn) + return mt_constructor(inds_e) end """ -DO a single update of a message tensor using the current subgraph and the incoming mts +Do a single update of a message tensor using the current subgraph and the incoming mts """ function update_message_tensor( - tn::ITensorNetwork, - subgraph_vertices::Vector, - mts::Vector{ITensorNetwork}; - contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1), + ptn::PartitionedGraph, + edge::PartitionEdge, + mts; + contractor=default_contractor, + contractor_kwargs=default_contractor_kwargs(), + mt_constructor=default_mt_constructor, ) - mts_itensors = reduce(vcat, ITensor.(mts); init=ITensor[]) + pb_edges = partitionedges(ptn, boundary_edges(ptn, vertices(ptn, src(edge)); dir=:in)) - contract_list = ITensor[mts_itensors; ITensor[tn[v] for v in subgraph_vertices]] - tn = if isone(length(contract_list)) - copy(only(contract_list)) - else - ITensorNetwork(contract_list) - end + incoming_messages = [ + e_in ∈ keys(mts) ? mts[e_in] : message_tensor(ptn, e_in; mt_constructor) for + e_in in setdiff(pb_edges, [reverse(edge)]) + ] + incoming_messages = reduce(vcat, incoming_messages; init=ITensor[]) - if contract_kwargs.alg != "exact" - contract_output = contract(tn; contract_kwargs...) - else - contract_output = contract(tn; sequence=contraction_sequence(tn; alg="optimal")) - end + contract_list = ITensor[ + incoming_messages + Vector{ITensor}(subgraph(ptn, src(edge))) + ] - itn = if typeof(contract_output) == ITensor - ITensorNetwork(contract_output) - else - first(contract_output) - end - normalize!.(vertex_data(itn)) - - return itn -end - -function update_message_tensor( - tn::ITensorNetwork, subgraph::ITensorNetwork, mts::Vector{ITensorNetwork}; kwargs... -) - return update_message_tensor(tn, vertices(subgraph), mts; kwargs...) + return contractor(contract_list; contractor_kwargs...) end """ Do a sequential update of message tensors on `edges` for a given ITensornetwork and its partition into sub graphs """ function belief_propagation_iteration( - tn::ITensorNetwork, - mts::DataGraph, - edges::Vector{<:AbstractEdge}; - contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1), - compute_norm=false, + ptn::PartitionedGraph, mts, edges::Vector{<:PartitionEdge}; compute_norm=false, kwargs... ) new_mts = copy(mts) c = 0 for e in edges - environment_tensornetworks = ITensorNetwork[ - new_mts[e_in] for - e_in in setdiff(boundary_edges(new_mts, [src(e)]; dir=:in), [reverse(e)]) - ] - new_mts[src(e) => dst(e)] = update_message_tensor( - tn, new_mts[src(e)], environment_tensornetworks; contract_kwargs - ) + new_mts[e] = update_message_tensor(ptn, e, new_mts; kwargs...) if compute_norm - LHS, RHS = ITensors.contract(ITensor(mts[src(e) => dst(e)])), - ITensors.contract(ITensor(new_mts[src(e) => dst(e)])) + LHS = e ∈ keys(mts) ? contract(mts[e]) : contract(message_tensor(ptn, e)) + RHS = contract(new_mts[e]) + #This line only makes sense if the message tensors are rank 2??? Should fix this. LHS /= sum(diag(LHS)) RHS /= sum(diag(RHS)) c += 0.5 * norm(denseblocks(LHS) - denseblocks(RHS)) @@ -113,18 +70,12 @@ Currently we send the full message tensor data struct to belief_propagation_iter mts relevant to that subgraph. """ function belief_propagation_iteration( - tn::ITensorNetwork, - mts::DataGraph, - edge_groups::Vector{<:Vector{<:AbstractEdge}}; - contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1), - compute_norm=false, + ptn::PartitionedGraph, mts, edge_groups::Vector{<:Vector{<:PartitionEdge}}; kwargs... ) new_mts = copy(mts) c = 0 for edges in edge_groups - updated_mts, ct = belief_propagation_iteration( - tn, mts, edges; contract_kwargs, compute_norm - ) + updated_mts, ct = belief_propagation_iteration(ptn, mts, edges; kwargs...) for e in edges new_mts[e] = updated_mts[e] end @@ -134,30 +85,26 @@ function belief_propagation_iteration( end function belief_propagation_iteration( - tn::ITensorNetwork, - mts::DataGraph; - contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1), - compute_norm=false, - edges=edge_sequence(mts), + ptn::PartitionedGraph, mts; edges=default_edge_sequence(ptn), kwargs... ) - return belief_propagation_iteration(tn, mts, edges; contract_kwargs, compute_norm) + return belief_propagation_iteration(ptn, mts, edges; kwargs...) end function belief_propagation( - tn::ITensorNetwork, - mts::DataGraph; - contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1), - niters=default_bp_niters(mts), + ptn::PartitionedGraph, + mts; + niters=default_bp_niters(partitioned_graph(ptn)), target_precision=nothing, - edges=edge_sequence(mts), + edges=default_edge_sequence(ptn), verbose=false, + kwargs..., ) compute_norm = !isnothing(target_precision) if isnothing(niters) error("You need to specify a number of iterations for BP!") end for i in 1:niters - mts, c = belief_propagation_iteration(tn, mts, edges; contract_kwargs, compute_norm) + mts, c = belief_propagation_iteration(ptn, mts, edges; compute_norm, kwargs...) if compute_norm && c <= target_precision if verbose println("BP converged to desired precision after $i iterations.") @@ -168,49 +115,29 @@ function belief_propagation( return mts end -function belief_propagation( - tn::ITensorNetwork; - contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1), - nvertices_per_partition=nothing, - npartitions=nothing, - subgraph_vertices=nothing, - niters=default_bp_niters(mts), - target_precision=nothing, - verbose=false, -) - mts = message_tensors(tn; nvertices_per_partition, npartitions, subgraph_vertices) - return belief_propagation(tn, mts; contract_kwargs, niters, target_precision, verbose) +function belief_propagation(ptn::PartitionedGraph; bp_cache=default_bp_cache, kwargs...) + mts = bp_cache(ptn) + return belief_propagation(ptn, mts; kwargs...) end - """ -Given a subet of vertices of a given Tensor Network and the Message Tensors for that network, return a Dictionary with the involved subgraphs as keys and the vector of tensors associated with that subgraph as values -Specifically, the contraction of the environment tensors and tn[vertices] will be a scalar. +Given a subet of partitionvertices of a ptn get the incoming message tensors to that region """ -function get_environment(tn::ITensorNetwork, mts::DataGraph, verts::Vector; dir=:in) - subgraphs = unique([find_subgraph(v, mts) for v in verts]) - - if dir == :out - return get_environment(tn, mts, setdiff(vertices(tn), verts)) - end - - env_tns = ITensorNetwork[mts[e] for e in boundary_edges(mts, subgraphs; dir=:in)] - central_tn = ITensorNetwork( - ITensor[tn[v] for v in setdiff(flatten([vertices(mts[s]) for s in subgraphs]), verts)] - ) - return ITensorNetwork(vcat(env_tns, ITensorNetwork[central_tn])) +function environment_tensors(ptn::PartitionedGraph, mts, verts::Vector) + partition_verts = partitionvertices(ptn, verts) + central_verts = vertices(ptn, partition_verts) + + pedges = partitionedges(ptn, boundary_edges(ptn, central_verts; dir=:in)) + env_tensors = [mts[e] for e in pedges] + env_tensors = reduce(vcat, env_tensors; init=ITensor[]) + central_tensors = ITensor[ + (unpartitioned_graph(ptn))[v] for v in setdiff(central_verts, verts) + ] + + return vcat(env_tensors, central_tensors) end -""" -Calculate the contraction of a tensor network centred on the vertices verts. Using message tensors. -Defaults to using tn[verts] as the local network but can be overriden -""" -function approx_network_region( - tn::ITensorNetwork, - mts::DataGraph, - verts::Vector; - verts_tn=ITensorNetwork([tn[v] for v in verts]), +function environment_tensors( + ptn::PartitionedGraph, mts, partition_verts::Vector{<:PartitionVertex} ) - environment_tn = get_environment(tn, mts, verts) - - return environment_tn ⊗ verts_tn + return environment_tensors(ptn, mts, vertices(ptn, partition_verts)) end diff --git a/src/beliefpropagation/beliefpropagation_schedule.jl b/src/beliefpropagation/beliefpropagation_schedule.jl index cbc5f206..a56814cd 100644 --- a/src/beliefpropagation/beliefpropagation_schedule.jl +++ b/src/beliefpropagation/beliefpropagation_schedule.jl @@ -1,4 +1,7 @@ default_edge_sequence_alg() = "forest_cover" +function default_edge_sequence(pg::PartitionedGraph) + return PartitionEdge.(edge_sequence(partitioned_graph(pg))) +end @traitfn default_bp_niters(g::::(!IsDirected)) = is_tree(g) ? 1 : nothing @traitfn function default_bp_niters(g::::IsDirected) diff --git a/src/beliefpropagation/sqrt_beliefpropagation.jl b/src/beliefpropagation/sqrt_beliefpropagation.jl deleted file mode 100644 index 20ff3d8e..00000000 --- a/src/beliefpropagation/sqrt_beliefpropagation.jl +++ /dev/null @@ -1,149 +0,0 @@ -# using ITensors: scalartype -# using ITensorNetworks: find_subgraph, map_diag, sqrt_diag, boundary_edges - -function sqrt_belief_propagation_iteration( - tn::ITensorNetwork, sqrt_mts::DataGraph, edges::Vector{<:AbstractEdge} -) - new_sqrt_mts = copy(sqrt_mts) - c = 0.0 - for e in edges - environment_tensornetworks = ITensorNetwork[ - new_sqrt_mts[e_in] for - e_in in setdiff(boundary_edges(new_sqrt_mts, [src(e)]; dir=:in), [reverse(e)]) - ] - - new_sqrt_mts[src(e) => dst(e)] = update_sqrt_message_tensor( - tn, new_sqrt_mts[src(e)], environment_tensornetworks; - ) - - # if compute_norm - # LHS, RHS = ITensors.contract(ITensor(sqrt_mts[src(e) => dst(e)])), - # ITensors.contract(ITensor(new_sqrt_mts[src(e) => dst(e)])) - # LHS /= sum(diag(LHS)) - # RHS /= sum(diag(RHS)) - # c += 0.5 * norm(LHS - RHS) - # end - end - return new_sqrt_mts, c / (length(edges)) -end - -function sqrt_belief_propagation_iteration( - tn::ITensorNetwork, sqrt_mts::DataGraph, edges::Vector{<:Vector{<:AbstractEdge}} -) - new_sqrt_mts = copy(sqrt_mts) - c = 0.0 - for e_group in edges - updated_sqrt_mts, ct = sqrt_belief_propagation_iteration(tn, sqr_mts, e_group) - for e in e_group - new_sqrt_mts[e] = updated_sqrt_mts[e] - end - c += ct - end - return new_sqrt_mts, c / (length(edges)) -end - -function sqrt_belief_propagation_iteration( - tn::ITensorNetwork, sqrt_mts::DataGraph; edges=edge_sequence(mts) -) - return sqrt_belief_propagation_iteration(tn, sqrt_mts, edges) -end - -function sqrt_belief_propagation( - tn::ITensorNetwork, - mts::DataGraph; - niters=default_bp_niters(mts), - edges=edge_sequence(mts), - # target_precision::Union{Float64,Nothing}=nothing, -) - # compute_norm = target_precision == nothing ? false : true - sqrt_mts = sqrt_message_tensors(tn, mts) - if isnothing(niters) - error("You need to specify a number of iterations for BP!") - end - for i in 1:niters - sqrt_mts, c = sqrt_belief_propagation_iteration(tn, sqrt_mts, edges) #; compute_norm) - # if compute_norm && c <= target_precision - # println( - # "Belief Propagation finished. Reached a canonicalness of " * - # string(c) * - # " after $i iterations. ", - # ) - # break - # end - end - return sqr_message_tensors(sqrt_mts) -end - -function update_sqrt_message_tensor( - tn::ITensorNetwork, subgraph_vertices::Vector, sqrt_mts::Vector{ITensorNetwork}; -) - sqrt_mts_itensors = reduce(vcat, ITensor.(sqrt_mts); init=ITensor[]) - v = only(unique(first.(subgraph_vertices))) - site_itensor = tn[v] - contract_list = ITensor[sqrt_mts_itensors; site_itensor] - contract_output = contract( - contract_list; sequence=contraction_sequence(contract_list; alg="optimal") - ) - left_inds = [uniqueinds(contract_output, site_itensor); siteinds(tn, v)] - Q, R = qr(contract_output, left_inds) - normalize!(R) - return ITensorNetwork(R) -end - -function update_sqrt_message_tensor( - tn::ITensorNetwork, subgraph::ITensorNetwork, mts::Vector{ITensorNetwork}; kwargs... -) - return update_sqrt_message_tensor(tn, vertices(subgraph), mts; kwargs...) -end - -function sqrt_message_tensors( - ψ::ITensorNetwork, - mts::DataGraph; - eigen_message_tensor_cutoff=10 * eps(real(scalartype(ψ))), - regularization=10 * eps(real(scalartype(ψ))), -) - sqrt_mts = copy(mts) - for e in edges(ψ) - vsrc, vdst = src(e), dst(e) - ψvsrc, ψvdst = ψ[vsrc], ψ[vdst] - - s1, s2 = find_subgraph((vsrc, 1), mts), find_subgraph((vdst, 1), mts) - edge_ind = commoninds(ψ[vsrc], ψ[vdst]) - edge_ind_sim = sim(edge_ind) - - X_D, X_U = eigen( - contract(ITensor(mts[s1 => s2])); ishermitian=true, cutoff=eigen_message_tensor_cutoff - ) - Y_D, Y_U = eigen( - contract(ITensor(mts[s2 => s1])); ishermitian=true, cutoff=eigen_message_tensor_cutoff - ) - X_D, Y_D = map_diag(x -> x + regularization, X_D), - map_diag(x -> x + regularization, Y_D) - - rootX_D, rootY_D = sqrt_diag(X_D), sqrt_diag(Y_D) - rootX = X_U * rootX_D * prime(dag(X_U)) - rootY = Y_U * rootY_D * prime(dag(Y_U)) - - sqrt_mts[s1 => s2] = ITensorNetwork(rootX) - sqrt_mts[s2 => s1] = ITensorNetwork(rootY) - end - return sqrt_mts -end - -function sqr_message_tensors(sqrt_mts::DataGraph) - mts = copy(sqrt_mts) - for e in edges(sqrt_mts) - sqrt_mt_tn = sqrt_mts[e] - sqrt_mt = sqrt_mt_tn[only(vertices(sqrt_mt_tn))] - sqrt_mt_rev_tn = sqrt_mts[reverse(e)] - sqrt_mt_rev = sqrt_mt_rev_tn[only(vertices(sqrt_mt_rev_tn))] - l = commoninds(sqrt_mt, sqrt_mt_rev) - mt = dag(prime(sqrt_mt, l)) * sqrt_mt - normalize!(mt) - mt_rev = dag(prime(sqrt_mt_rev, l)) * sqrt_mt_rev - normalize!(mt_rev) - mts[e] = ITensorNetwork(mt) - mts[reverse(e)] = ITensorNetwork(mt_rev) - end - return mts -end diff --git a/src/boundarymps.jl b/src/boundarymps.jl index 449442d6..8e27868e 100644 --- a/src/boundarymps.jl +++ b/src/boundarymps.jl @@ -1,4 +1,4 @@ -#GIVEN AN ITENSOR NETWORK CORRESPONDING to an Lx*Ly grid with sites indexed as (i,j) then perform contraction using a sequence of mps-mpo contractions +#Given an ITensorNetwork on an Lx*Ly grid with sites indexed as (i,j) then perform contraction using a sequence of mps-mpo contractions function contract_boundary_mps(tn::ITensorNetwork; kwargs...) dims = maximum(vertices(tn)) d1, d2 = dims diff --git a/src/contract.jl b/src/contract.jl index 97265d77..77c4ee1a 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -17,3 +17,30 @@ function contract( ) return approx_itensornetwork(alg, tn, output_structure; kwargs...) end + +function contract_density_matrix( + contract_list::Vector{ITensor}; normalize=true, contractor_kwargs... +) + tn, _ = contract( + ITensorNetwork(contract_list); alg="density_matrix", contractor_kwargs... + ) + out = Vector{ITensor}(tn) + if normalize + out .= normalize!.(copy.(out)) + end + return out +end + +function contract_exact( + contract_list::Vector{ITensor}; + contraction_sequence_alg="optimal", + normalize=true, + contractor_kwargs..., +) + seq = contraction_sequence(contract_list; alg=contraction_sequence_alg) + out = ITensors.contract(contract_list; sequence=seq, contractor_kwargs...) + if normalize + normalize!(out) + end + return ITensor[out] +end diff --git a/src/exports.jl b/src/exports.jl index f298994d..97194466 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -51,6 +51,7 @@ export AbstractITensorNetwork, itensors, reverse_bfs_edges, data_graph, + default_bp_cache, flatten_networks, inner_network, norm_network, diff --git a/src/gauging.jl b/src/gauging.jl index 00347046..c6d3c3f6 100644 --- a/src/gauging.jl +++ b/src/gauging.jl @@ -1,10 +1,10 @@ """initialize bond tensors of an ITN to identity matrices""" function initialize_bond_tensors(ψ::ITensorNetwork; index_map=prime) - bond_tensors = DataGraph{vertextype(ψ),ITensor,ITensor}(underlying_graph(ψ)) + bond_tensors = DataGraph{vertextype(ψ),Nothing,ITensor}(underlying_graph(ψ)) for e in edges(ψ) index = commoninds(ψ[src(e)], ψ[dst(e)]) - bond_tensors[e] = dense(delta(index, index_map(index))) + bond_tensors[e] = denseblocks(delta(index, index_map(index))) end return bond_tensors @@ -13,7 +13,8 @@ end """Use an ITensorNetwork ψ, its bond tensors and gauging mts to put ψ into the vidal gauge, return the bond tensors and ψ_vidal.""" function vidal_gauge( ψ::ITensorNetwork, - mts::DataGraph, + pψψ::PartitionedGraph, + mts, bond_tensors::DataGraph; eigen_message_tensor_cutoff=10 * eps(real(scalartype(ψ))), regularization=10 * eps(real(scalartype(ψ))), @@ -26,15 +27,13 @@ function vidal_gauge( vsrc, vdst = src(e), dst(e) ψvsrc, ψvdst = ψ_vidal[vsrc], ψ_vidal[vdst] - s1, s2 = find_subgraph((vsrc, 1), mts), find_subgraph((vdst, 1), mts) - edge_ind = commoninds(ψ_vidal[vsrc], ψ_vidal[vdst]) + pe = partitionedge(pψψ, (vsrc, 1) => (vdst, 1)) + edge_ind = commoninds(ψvsrc, ψvdst) edge_ind_sim = sim(edge_ind) - X_D, X_U = eigen( - ITensor(mts[s1 => s2]); ishermitian=true, cutoff=eigen_message_tensor_cutoff - ) + X_D, X_U = eigen(only(mts[pe]); ishermitian=true, cutoff=eigen_message_tensor_cutoff) Y_D, Y_U = eigen( - ITensor(mts[s2 => s1]); ishermitian=true, cutoff=eigen_message_tensor_cutoff + only(mts[reverse(pe)]); ishermitian=true, cutoff=eigen_message_tensor_cutoff ) X_D, Y_D = map_diag(x -> x + regularization, X_D), map_diag(x -> x + regularization, Y_D) @@ -77,7 +76,8 @@ end """Use an ITensorNetwork ψ in the symmetric gauge and its mts to put ψ into the vidal gauge. Return the bond tensors and ψ_vidal.""" function vidal_gauge( ψ::ITensorNetwork, - mts::DataGraph; + pψψ::PartitionedGraph, + mts; eigen_message_tensor_cutoff=10 * eps(real(scalartype(ψ))), regularization=10 * eps(real(scalartype(ψ))), edges=NamedGraphs.edges(ψ), @@ -85,7 +85,14 @@ function vidal_gauge( ) bond_tensors = initialize_bond_tensors(ψ) return vidal_gauge( - ψ, mts, bond_tensors; eigen_message_tensor_cutoff, regularization, edges, svd_kwargs... + ψ, + pψψ, + mts, + bond_tensors; + eigen_message_tensor_cutoff, + regularization, + edges, + svd_kwargs..., ) end @@ -100,41 +107,36 @@ function vidal_gauge( svd_kwargs..., ) ψψ = norm_network(ψ) - Z = partition(ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ))))) - mts = message_tensors(Z) + pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) - mts = belief_propagation( - ψψ, - mts; - contract_kwargs=(; alg="exact"), - niters, - target_precision=target_canonicalness, - verbose, + mts = belief_propagation(pψψ; niters, target_precision=target_canonicalness, verbose) + return vidal_gauge( + ψ, pψψ, mts; eigen_message_tensor_cutoff, regularization, svd_kwargs... ) - return vidal_gauge(ψ, mts; eigen_message_tensor_cutoff, regularization, svd_kwargs...) end -"""Transform from an ITensor in the Vidal Gauge (bond tensors) to the Symmetric Gauge (message tensors)""" +"""Transform from an ITensor in the Vidal Gauge (bond tensors) to the Symmetric Gauge (partitionedgraph, message tensors)""" function vidal_to_symmetric_gauge(ψ::ITensorNetwork, bond_tensors::DataGraph) ψsymm = copy(ψ) ψψsymm = norm_network(ψsymm) - Z = partition( - ψψsymm; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψsymm)))) - ) - ψsymm_mts = message_tensors_skeleton(Z) + pψψsymm = PartitionedGraph(ψψsymm, group(v -> v[1], vertices(ψψsymm))) + ψsymm_mts = default_bp_cache(pψψsymm) for e in edges(ψsymm) vsrc, vdst = src(e), dst(e) - s1, s2 = find_subgraph((vsrc, 1), ψsymm_mts), find_subgraph((vdst, 1), ψsymm_mts) + pe = partitionedge(pψψsymm, NamedEdge((vsrc, 1) => (vdst, 1))) root_S = sqrt_diag(bond_tensors[e]) setindex_preserve_graph!(ψsymm, noprime(root_S * ψsymm[vsrc]), vsrc) setindex_preserve_graph!(ψsymm, noprime(root_S * ψsymm[vdst]), vdst) - ψsymm_mts[s1 => s2], ψsymm_mts[s2 => s1] = ITensorNetwork(bond_tensors[e]), - ITensorNetwork(bond_tensors[e]) + ψsymm_mts[pe], ψsymm_mts[reverse(pe)] = copy(ITensor[dense(bond_tensors[e])]), + copy(ITensor[dense(bond_tensors[e])]) end - return ψsymm, ψsymm_mts + ψψsymm = norm_network(ψsymm) + pψψsymm = PartitionedGraph(ψψsymm, group(v -> v[1], vertices(ψψsymm))) + + return ψsymm, pψψsymm, ψsymm_mts end """Put an ITensorNetwork into the symmetric gauge and also return the message tensors (which are the diagonal bond matrices from the Vidal Gauge)""" @@ -146,7 +148,7 @@ function symmetric_gauge( target_canonicalness::Union{Nothing,Float64}=nothing, svd_kwargs..., ) - ψsymm, bond_tensors = vidal_gauge( + ψ_vidal, bond_tensors = vidal_gauge( ψ; eigen_message_tensor_cutoff, regularization, @@ -155,21 +157,24 @@ function symmetric_gauge( svd_kwargs..., ) - return vidal_to_symmetric_gauge(ψsymm, bond_tensors) + return vidal_to_symmetric_gauge(ψ_vidal, bond_tensors) end """Transform from the Symmetric Gauge (message tensors) to the Vidal Gauge (bond tensors)""" function symmetric_to_vidal_gauge( - ψ::ITensorNetwork, mts::DataGraph; regularization=10 * eps(real(scalartype(ψ))) + ψ::ITensorNetwork, + pψψ::PartitionedGraph, + mts; + regularization=10 * eps(real(scalartype(ψ))), ) - bond_tensors = DataGraph{vertextype(ψ),ITensor,ITensor}(underlying_graph(ψ)) + bond_tensors = DataGraph{vertextype(ψ),Nothing,ITensor}(underlying_graph(ψ)) ψ_vidal = copy(ψ) for e in edges(ψ) vsrc, vdst = src(e), dst(e) - s1, s2 = find_subgraph((vsrc, 1), mts), find_subgraph((vdst, 1), mts) - bond_tensors[e] = ITensor(mts[s1 => s2]) + pe = partitionedge(pψψ, NamedEdge((vsrc, 1) => (vdst, 1))) + bond_tensors[e], bond_tensors[reverse(e)] = only(mts[pe]), only(mts[pe]) invroot_S = invsqrt_diag(map_diag(x -> x + regularization, bond_tensors[e])) setindex_preserve_graph!(ψ_vidal, noprime(invroot_S * ψ_vidal[vsrc]), vsrc) setindex_preserve_graph!(ψ_vidal, noprime(invroot_S * ψ_vidal[vdst]), vdst) @@ -184,7 +189,7 @@ function vidal_itn_isometries( bond_tensors::DataGraph; edges=vcat(NamedGraphs.edges(ψ), reverse.(NamedGraphs.edges(ψ))), ) - isometries = DataGraph{vertextype(ψ),ITensor,ITensor}(directed_graph(underlying_graph(ψ))) + isometries = Dict() for e in edges vsrc, vdst = src(e), dst(e) @@ -207,19 +212,19 @@ function vidal_itn_canonicalness(ψ::ITensorNetwork, bond_tensors::DataGraph) isometries = vidal_itn_isometries(ψ, bond_tensors) - for e in edges(isometries) + for e in keys(isometries) LHS = isometries[e] / sum(diag(isometries[e])) id = dense(delta(inds(LHS))) id /= sum(diag(id)) f += 0.5 * norm(id - LHS) end - return f / (length(edges(isometries))) + return f / (length(keys(isometries))) end """Function to measure the 'canonicalness' of a state in the Symmetric Gauge""" -function symmetric_itn_canonicalness(ψ::ITensorNetwork, mts::DataGraph) - ψ_vidal, bond_tensors = symmetric_to_vidal_gauge(ψ, mts) +function symmetric_itn_canonicalness(ψ::ITensorNetwork, pψψ::PartitionedGraph, mts) + ψ_vidal, bond_tensors = symmetric_to_vidal_gauge(ψ, pψψ, mts) return vidal_itn_canonicalness(ψ_vidal, bond_tensors) end diff --git a/src/itensornetwork.jl b/src/itensornetwork.jl index 0b8178c3..6c3e346d 100644 --- a/src/itensornetwork.jl +++ b/src/itensornetwork.jl @@ -264,11 +264,6 @@ end ITensorNetwork(itns::Vector{ITensorNetwork}) = reduce(⊗, itns) -function ITensors.ITensor(ψ::ITensorNetwork) - vs = vertices(ψ) - if length(vs) == 1 - return ψ[first(vs)] - else - return ITensor[ψ[v] for v in vertices(ψ)] - end +function Vector{ITensor}(ψ::ITensorNetwork) + return ITensor[ψ[v] for v in vertices(ψ)] end diff --git a/src/partition.jl b/src/partition.jl deleted file mode 100644 index 0709b77d..00000000 --- a/src/partition.jl +++ /dev/null @@ -1,341 +0,0 @@ -""" -Graph partitioning backend -""" -struct Backend{T} end - -Backend(s::Symbol) = Backend{s}() -Backend(s::String) = Backend(Symbol(s)) -Backend(backend::Backend) = backend - -macro Backend_str(s) - return :(Backend{$(Expr(:quote, Symbol(s)))}) -end - -""" -Current default graph partitioning backend -""" -const CURRENT_PARTITIONING_BACKEND = Ref{Union{Missing,Backend}}(missing) - -""" -Get the graph partitioning backend -""" -current_partitioning_backend() = CURRENT_PARTITIONING_BACKEND[] - -""" -Set the graph partitioning backend -""" -function set_partitioning_backend!(backend::Union{Missing,Backend,String}) - CURRENT_PARTITIONING_BACKEND[] = Backend(backend) - return nothing -end - -# KaHyPar configuration options -# -# configurations = readdir(joinpath(pkgdir(KaHyPar), "src", "config")) -# "cut_kKaHyPar_sea20.ini" -# "cut_rKaHyPar_sea20.ini" -# "km1_kKaHyPar-E_sea20.ini" -# "km1_kKaHyPar_eco_sea20.ini" -# "km1_kKaHyPar_sea20.ini" -# "km1_rKaHyPar_sea20.ini" -# -const kahypar_configurations = Dict([ - (objective="edge_cut", alg="kway") => "cut_kKaHyPar_sea20.ini", - (objective="edge_cut", alg="recursive") => "cut_rKaHyPar_sea20.ini", - (objective="connectivity", alg="kway") => "km1_kKaHyPar_sea20.ini", - (objective="connectivity", alg="recursive") => "km1_rKaHyPar_sea20.ini", -]) - -# Metis configuration options -const metis_algs = Dict(["kway" => :KWAY, "recursive" => :RECURSIVE]) - -function _npartitions( - g::AbstractGraph, npartitions::Integer, nvertices_per_partition::Nothing -) - return npartitions -end - -function _npartitions( - g::AbstractGraph, npartitions::Nothing, nvertices_per_partition::Integer -) - return nv(g) ÷ nvertices_per_partition -end - -function _npartitions(g::AbstractGraph, npartitions::Int, nvertices_per_partition::Int) - return error("Can't specify both `npartitions` and `nvertices_per_partition`") -end - -function _npartitions( - g::AbstractGraph, npartitions::Nothing, nvertices_per_partition::Nothing -) - return error("Must specify either `npartitions` or `nvertices_per_partition`") -end - -function subgraph_vertices( - g::Graph; - npartitions=nothing, - nvertices_per_partition=nothing, - backend=current_partitioning_backend(), - kwargs..., -) - #Metis cannot handle the edge case npartitions = 1, so we will fix it here for now - if (_npartitions(g, npartitions, nvertices_per_partition) == 1) - return group(v -> 1, collect(vertices(g))) - end - - return subgraph_vertices( - Backend(backend), g, _npartitions(g, npartitions, nvertices_per_partition); kwargs... - ) -end - -function subgraph_vertices( - g::NamedGraph; npartitions=nothing, nvertices_per_partition=nothing, kwargs... -) - vertex_partitions = subgraph_vertices( - parent_graph(g); npartitions, nvertices_per_partition, kwargs... - ) - #[inv(vertex_to_parent_vertex(g))[v] for v in partitions] - # TODO: output the reverse of this dictionary (a Vector of Vector - # of the vertices in each partition). - # return Dictionary(vertices(g), partitions) - return [ - parent_vertices_to_vertices(g, vertex_partition) for - vertex_partition in vertex_partitions - ] -end - -function subgraph_vertices( - g::AbstractDataGraph; npartitions=nothing, nvertices_per_partition=nothing, kwargs... -) - return subgraph_vertices( - underlying_graph(g); npartitions, nvertices_per_partition, kwargs... - ) -end - -""" - partition_vertices(g::AbstractGraph, subgraph_vertices::Vector) - -Given a graph (`g`) and groups of vertices defining subgraphs of that -graph (`subgraph_vertices`), return a DataGraph storing the subgraph -vertices on the vertices of the graph and with edges denoting -which subgraphs of the original graph have edges connecting them, along with -edge data storing the original edges that were connecting the subgraphs. -""" -function partition_vertices(g::AbstractGraph, subgraph_vertices) - partitioned_vertices = DataGraph( - NamedGraph(eachindex(subgraph_vertices)), Dictionary(subgraph_vertices) - ) - for e in edges(g) - s1 = findfirst_on_vertices( - subgraph_vertices -> src(e) ∈ subgraph_vertices, partitioned_vertices - ) - s2 = findfirst_on_vertices( - subgraph_vertices -> dst(e) ∈ subgraph_vertices, partitioned_vertices - ) - if (!has_edge(partitioned_vertices, s1, s2) && s1 ≠ s2) - add_edge!(partitioned_vertices, s1, s2) - partitioned_vertices[s1 => s2] = Vector{edgetype(g)}() - end - if has_edge(partitioned_vertices, s1, s2) - push!(partitioned_vertices[s1 => s2], e) - end - end - return partitioned_vertices -end - -""" - partition_vertices(g::AbstractGraph; npartitions, nvertices_per_partition, kwargs...) - -Given a graph `g`, partition the vertices of `g` into 'npartitions' partitions -or into partitions with `nvertices_per_partition` vertices per partition. -Try to keep all subgraphs the same size and minimise edges cut between them -Returns a datagraph where each vertex contains the list of vertices involved in that subgraph. The edges state which subgraphs are connected. -A graph partitioning backend such as Metis or KaHyPar needs to be installed for this function to work. -""" -function partition_vertices( - g::AbstractGraph; npartitions=nothing, nvertices_per_partition=nothing, kwargs... -) - return partition_vertices( - g, subgraph_vertices(g; npartitions, nvertices_per_partition, kwargs...) - ) -end - -""" -Find all vertices `v` such that `f(graph[v]) == true` -""" -function findall_on_vertices(f::Function, graph::AbstractDataGraph) - return findall(f, vertex_data(graph)) -end - -""" -Find the vertex `v` such that `f(graph[v]) == true` -""" -function findfirst_on_vertices(f::Function, graph::AbstractDataGraph) - return findfirst(f, vertex_data(graph)) -end - -""" -Find all edges `e` such that `f(graph[e]) == true` -""" -function findall_on_edges(f::Function, graph::AbstractDataGraph) - return findall(f, edge_data(graph)) -end - -""" -Find the edge `e` such that `f(graph[e]) == true` -""" -function findfirst_on_edges(f::Function, graph::AbstractDataGraph) - return findfirst(f, edge_data(graph)) -end - -""" -Find the subgraph which contains the specified vertex. - -TODO: Rename something more general, like: - -findfirst_in_vertex_data(item, graph::AbstractDataGraph) -""" -function find_subgraph(vertex, subgraphs::DataGraph) - return findfirst_on_vertices(subgraph -> vertex ∈ vertices(subgraph), subgraphs) -end - -function subgraphs(g::AbstractSimpleGraph, subgraph_vertices) - return subgraphs(NamedGraph(g), subgraph_vertices) -end - -""" - subgraphs(g::AbstractGraph, subgraph_vertices) - -Return a collection of subgraphs of `g` defined by the subgraph -vertices `subgraph_vertices`. -""" -function subgraphs(g::AbstractGraph, subgraph_vertices) - return map(vs -> subgraph(g, vs), subgraph_vertices) -end - -""" - subgraphs(g::AbstractGraph; npartitions::Integer, kwargs...) - -Given a graph `g`, partition `g` into `npartitions` partitions -or into partitions with `nvertices_per_partition` vertices per partition, -returning a list of subgraphs. -Try to keep all subgraphs the same size and minimise edges cut between them. -A graph partitioning backend such as Metis or KaHyPar needs to be installed for this function to work. -""" -function subgraphs( - g::AbstractGraph; npartitions=nothing, nvertices_per_partition=nothing, kwargs... -) - return subgraphs(g, subgraph_vertices(g; npartitions, nvertices_per_partition, kwargs...)) -end - -function partition(g::AbstractSimpleGraph, subgraph_vertices) - return partition(NamedGraph(g), subgraph_vertices) -end - -function partition(g::AbstractGraph, subgraph_vertices) - partitioned_graph = DataGraph( - NamedGraph(eachindex(subgraph_vertices)), subgraphs(g, Dictionary(subgraph_vertices)) - ) - for e in edges(g) - s1 = findfirst_on_vertices(subgraph -> src(e) ∈ vertices(subgraph), partitioned_graph) - s2 = findfirst_on_vertices(subgraph -> dst(e) ∈ vertices(subgraph), partitioned_graph) - if (!has_edge(partitioned_graph, s1, s2) && s1 ≠ s2) - add_edge!(partitioned_graph, s1, s2) - partitioned_graph[s1 => s2] = Dictionary( - [:edges, :edge_data], - [Vector{edgetype(g)}(), Dictionary{edgetype(g),edge_data_type(g)}()], - ) - end - if has_edge(partitioned_graph, s1, s2) - push!(partitioned_graph[s1 => s2][:edges], e) - if isassigned(g, e) - set!(partitioned_graph[s1 => s2][:edge_data], e, g[e]) - end - end - end - return partitioned_graph -end - -""" - partition(g::AbstractGraph; npartitions::Integer, kwargs...) - partition(g::AbstractGraph, subgraph_vertices) - -Given a graph `g`, partition `g` into `npartitions` partitions -or into partitions with `nvertices_per_partition` vertices per partition. -The partitioning tries to keep all subgraphs the same size and minimize -edges cut between them. - -Alternatively, specify a desired partitioning with a collection of sugraph -vertices. - -Returns a data graph where each vertex contains the corresponding subgraph as vertex data. -The edges indicates which subgraphs are connected, and the edge data stores a dictionary -with two fields. The field `:edges` stores a list of the edges of the original graph -that were connecting the two subgraphs, and `:edge_data` stores a dictionary -mapping edges of the original graph to the data living on the edges of the original -graph, if it existed. - -Therefore, one should be able to extract that data and recreate the original -graph from the results of `partition`. - -A graph partitioning backend such as Metis or KaHyPar needs to be installed for this function to work -if the subgraph vertices aren't specified explicitly. -""" -function partition( - g::AbstractGraph; - npartitions=nothing, - nvertices_per_partition=nothing, - subgraph_vertices=nothing, - kwargs..., -) - if count(isnothing, (npartitions, nvertices_per_partition, subgraph_vertices)) != 2 - error( - "Error: Cannot give multiple/ no partitioning options. Please specify exactly one." - ) - end - - if isnothing(subgraph_vertices) - subgraph_vertices = ITensorNetworks.subgraph_vertices( - g; npartitions, nvertices_per_partition, kwargs... - ) - end - - return partition(g, subgraph_vertices) -end - -"""Given a nested graph fetch all the vertices down to the lowest levels and return the grouping -at the highest level. Keyword argument is used to state whether we are at the top""" -function nested_graph_leaf_vertices(g; toplevel=true) - vertex_groups = [] - for v in vertices(g) - if hasmethod(vertices, Tuple{typeof(g[v])}) - if !toplevel - push!(vertex_groups, nested_graph_leaf_vertices(g[v]; toplevel=false)...) - else - push!(vertex_groups, nested_graph_leaf_vertices(g[v]; toplevel=false)) - end - else - push!(vertex_groups, [v]...) - end - end - - return vertex_groups -end - -""" - TODO: do we want to make it a public function? -""" -function _noncommoninds(partition::DataGraph) - networks = [Vector{ITensor}(partition[v]) for v in vertices(partition)] - network = vcat(networks...) - return noncommoninds(network...) -end - -# Util functions for partition -function _commoninds(partition::DataGraph) - networks = [Vector{ITensor}(partition[v]) for v in vertices(partition)] - network = vcat(networks...) - outinds = noncommoninds(network...) - allinds = mapreduce(t -> [i for i in inds(t)], vcat, network) - return Vector(setdiff(allinds, outinds)) -end diff --git a/src/requires/kahypar.jl b/src/requires/kahypar.jl deleted file mode 100644 index fd4f7b38..00000000 --- a/src/requires/kahypar.jl +++ /dev/null @@ -1,33 +0,0 @@ -set_partitioning_backend!(Backend"KaHyPar"()) - -# https://github.com/kahypar/KaHyPar.jl/issues/20 -KaHyPar.HyperGraph(g::SimpleGraph) = incidence_matrix(g) - -""" - subgraph_vertices(::Backend"KaHyPar", g::Graph, npartiations::Integer; objective="edge_cut", alg="kway", kwargs...) - -- default_configuration => "cut_kKaHyPar_sea20.ini" -- :edge_cut => "cut_kKaHyPar_sea20.ini" -- :connectivity => "km1_kKaHyPar_sea20.ini" -- imbalance::Number=0.03 -""" -function subgraph_vertices( - ::Backend"KaHyPar", - g::SimpleGraph, - npartitions::Integer; - objective="edge_cut", - alg="kway", - configuration=nothing, - kwargs..., -) - if isnothing(configuration) - configuration = joinpath( - pkgdir(KaHyPar), - "src", - "config", - kahypar_configurations[(objective=objective, alg=alg)], - ) - end - partitions = @suppress KaHyPar.partition(g, npartitions; configuration, kwargs...) - return groupfind(partitions .+ 1) -end diff --git a/src/requires/metis.jl b/src/requires/metis.jl deleted file mode 100644 index 19985484..00000000 --- a/src/requires/metis.jl +++ /dev/null @@ -1,33 +0,0 @@ -set_partitioning_backend!(Backend"Metis"()) - -""" - subgraph_vertices(::Backend"Metis", g::AbstractGraph, npartitions::Integer; alg="recursive") - -Partition the graph `G` in `n` parts. -The partition algorithm is defined by the `alg` keyword: - - :KWAY: multilevel k-way partitioning - - :RECURSIVE: multilevel recursive bisection -""" -function subgraph_vertices( - ::Backend"Metis", g::SimpleGraph, npartitions::Integer; alg="recursive", kwargs... -) - metis_alg = metis_algs[alg] - partitions = Metis.partition(g, npartitions; alg=metis_alg, kwargs...) - return groupfind(Int.(partitions)) -end - -## #= -## Metis.partition(G, n; alg = :KWAY) -## -## Partition the graph `G` in `n` parts. -## The partition algorithm is defined by the `alg` keyword: -## - :KWAY: multilevel k-way partitioning -## - :RECURSIVE: multilevel recursive bisection -## =# -## function partition(g::Metis.Graph, npartitions::Integer) -## return Metis.partition(g, npartitions; alg=:KWAY) -## end -## -## function partition(g::Graph, npartitions::Integer) -## return partition(Metis.graph(adjacency_matrix(g)), npartitions) -## end diff --git a/test/test_apply.jl b/test/test_apply.jl index 64909adf..2df2f491 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,11 +1,10 @@ using ITensorNetworks using ITensorNetworks: belief_propagation, - get_environment, + environment_tensors, contract_inner, - message_tensors, - nested_graph_leaf_vertices, vidal_gauge, + vidal_apply, vidal_to_symmetric_gauge, norm_network using Test @@ -30,24 +29,16 @@ using SplitApplyCombine ψψ = norm_network(ψ) #Simple Belief Propagation Grouping - vertex_groupsSBP = nested_graph_leaf_vertices( - partition(partition(ψψ, group(v -> v[1], vertices(ψψ))); nvertices_per_partition=1) - ) - Z = partition(ψψ; subgraph_vertices=vertex_groupsSBP) - mtsSBP = message_tensors(Z) - mtsSBP = belief_propagation(ψψ, mtsSBP; contract_kwargs=(; alg="exact"), niters=50) - envsSBP = get_environment(ψψ, mtsSBP, [(v1, 1), (v1, 2), (v2, 1), (v2, 2)]) + pψψ_SBP = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) + mtsSBP = belief_propagation(pψψ_SBP; niters=20) + envsSBP = environment_tensors(pψψ_SBP, mtsSBP, PartitionVertex.([v1, v2])) - ψ_vidal, bond_tensors = vidal_gauge(ψ, mtsSBP) + ψ_vidal, bond_tensors = vidal_gauge(ψ, pψψ_SBP, mtsSBP) #This grouping will correspond to calculating the environments exactly (each column of the grid is a partition) - vertex_groupsGBP = nested_graph_leaf_vertices( - partition(partition(ψψ, group(v -> v[1][1], vertices(ψψ))); nvertices_per_partition=1) - ) - Z = partition(ψψ; subgraph_vertices=vertex_groupsSBP) - mtsGBP = message_tensors(Z) - mtsGBP = belief_propagation(ψψ, mtsGBP; contract_kwargs=(; alg="exact"), niters=50) - envsGBP = get_environment(ψψ, mtsGBP, [(v1, 1), (v1, 2), (v2, 1), (v2, 2)]) + pψψ_GBP = PartitionedGraph(ψψ, group(v -> v[1][1], vertices(ψψ))) + mtsGBP = belief_propagation(pψψ_GBP; niters=20) + envsGBP = environment_tensors(pψψ_GBP, mtsGBP, [(v1, 1), (v1, 2), (v2, 1), (v2, 2)]) ngates = 5 @@ -64,7 +55,9 @@ using SplitApplyCombine print_fidelity_loss=true, envisposdef=true, ) - ψOVidal, bond_tensors_t = apply(o, ψ_vidal, bond_tensors; maxdim=χ, normalize=true) + ψOVidal, bond_tensors_t = vidal_apply( + o, ψ_vidal, bond_tensors; maxdim=χ, normalize=true + ) ψOVidal_symm, _ = vidal_to_symmetric_gauge(ψOVidal, bond_tensors_t) ψOGBP = apply( o, diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index d3747598..cab105b9 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -2,21 +2,18 @@ using ITensorNetworks using ITensorNetworks: ising_network, belief_propagation, - nested_graph_leaf_vertices, - approx_network_region, split_index, contract_inner, contract_boundary_mps, - message_tensors + environment_tensors using Test using Compat using ITensors -#ASSUME ONE CAN INSTALL THIS, MIGHT FAIL FOR WINDOWS -using Metis using LinearAlgebra using NamedGraphs using SplitApplyCombine using Random +using Metis ITensors.disable_warn_order() @@ -38,20 +35,16 @@ ITensors.disable_warn_order() Oψ[v] = apply(op("Sz", s[v]), ψ[v]) exact_sz = contract_inner(Oψ, ψ) / contract_inner(ψ, ψ) - Z = partition(ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ))))) - mts = message_tensors(Z) - mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact")) - - numerator_network = approx_network_region( - ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork(ITensor[apply(op("Sz", s[v]), ψ[v])]) - ) - denominator_network = approx_network_region(ψψ, mts, [(v, 1)]) - bp_sz = ITensors.contract(numerator_network)[] / ITensors.contract(denominator_network)[] + pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) + mts = belief_propagation(pψψ) + env_tensors = environment_tensors(pψψ, mts, [PartitionVertex(v)]) + numerator = contract(vcat(env_tensors, ITensor[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))]))[] + denominator = contract(vcat(env_tensors, ITensor[ψ[v], op("I", s[v]), dag(prime(ψ[v]))]))[] - @test abs.(bp_sz - exact_sz) <= 1e-14 + @test abs.((numerator / denominator) - exact_sz) <= 1e-14 #Now test on a tree, should also be exact - g = named_comb_tree((6, 6)) + g = named_comb_tree((4, 4)) s = siteinds("S=1/2", g) χ = 2 Random.seed!(1564) @@ -65,19 +58,15 @@ ITensors.disable_warn_order() Oψ[v] = apply(op("Sz", s[v]), ψ[v]) exact_sz = contract_inner(Oψ, ψ) / contract_inner(ψ, ψ) - Z = partition(ψψ; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ))))) - mts = message_tensors(Z) - mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact")) - - numerator_network = approx_network_region( - ψψ, mts, [(v, 1)]; verts_tn=ITensorNetwork(ITensor[apply(op("Sz", s[v]), ψ[v])]) - ) - denominator_network = approx_network_region(ψψ, mts, [(v, 1)]) - bp_sz = ITensors.contract(numerator_network)[] / ITensors.contract(denominator_network)[] + pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) + mts = belief_propagation(pψψ) + env_tensors = environment_tensors(pψψ, mts, [PartitionVertex(v)]) + numerator = contract(vcat(env_tensors, ITensor[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))]))[] + denominator = contract(vcat(env_tensors, ITensor[ψ[v], op("I", s[v]), dag(prime(ψ[v]))]))[] - @test abs.(bp_sz - exact_sz) <= 1e-14 + @test abs.((numerator / denominator) - exact_sz) <= 1e-14 - #Now test two-site expec taking on the partition function of the Ising model. Not exact, but close + # #Now test two-site expec taking on the partition function of the Ising model. Not exact, but close g_dims = (3, 4) g = named_grid(g_dims) s = IndsNetwork(g; link_space=2) @@ -91,61 +80,46 @@ ITensors.disable_warn_order() ITensors.contract(ψOψ; sequence=contract_seq)[] / ITensors.contract(ψψ; sequence=contract_seq)[] - nsites = 2 - Z = partition(ψψ; nvertices_per_partition=nsites) - mts = message_tensors(Z) - mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"), niters=20) - numerator_network = approx_network_region( - ψψ, mts, vs; verts_tn=ITensorNetwork(ITensor[ψOψ[v] for v in vs]) - ) + pψψ = PartitionedGraph(ψψ; nvertices_per_partition=2, backend="Metis") + mts = belief_propagation(pψψ; niters=20) - denominator_network = approx_network_region(ψψ, mts, vs) - bp_szsz = - ITensors.contract(numerator_network)[] / ITensors.contract(denominator_network)[] + env_tensors = environment_tensors(pψψ, mts, vs) + numerator = contract(vcat(env_tensors, ITensor[ψOψ[v] for v in vs]))[] + denominator = contract(vcat(env_tensors, ITensor[ψψ[v] for v in vs]))[] - @test abs.(bp_szsz - actual_szsz) <= 0.05 + @test abs.((numerator / denominator) - actual_szsz) <= 0.05 - #Test forming a two-site RDM. Check it has the correct size, trace 1 and is PSD - g_dims = (4, 4) + # #Test forming a two-site RDM. Check it has the correct size, trace 1 and is PSD + g_dims = (3, 3) g = named_grid(g_dims) s = siteinds("S=1/2", g) - vs = [(2, 2), (2, 3), (2, 4)] + vs = [(2, 2), (2, 3)] χ = 3 ψ = randomITensorNetwork(s; link_space=χ) ψψ = ψ ⊗ prime(dag(ψ); sites=[]) - nsites = 2 - Zp = partition( - partition(ψψ, group(v -> v[1], vertices(ψψ))); nvertices_per_partition=nsites - ) - Zpp = partition(ψψ; subgraph_vertices=nested_graph_leaf_vertices(Zp)) - mts = message_tensors(Zpp) - mts = belief_propagation(ψψ, mts; contract_kwargs=(; alg="exact"), niters=20) + pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) + mts = belief_propagation(pψψ; niters=20) ψψsplit = split_index(ψψ, NamedEdge.([(v, 1) => (v, 2) for v in vs])) + env_tensors = environment_tensors(pψψ, mts, [(v, 2) for v in vs]) rdm = ITensors.contract( - approx_network_region( - ψψ, - mts, - [(v, 2) for v in vs]; - verts_tn=ITensorNetwork(ITensor[ψψsplit[vp] for vp in [(v, 2) for v in vs]]), - ), + vcat(env_tensors, ITensor[ψψsplit[vp] for vp in [(v, 2) for v in vs]]) ) rdm = array((rdm * combiner(inds(rdm; plev=0)...)) * combiner(inds(rdm; plev=1)...)) rdm /= tr(rdm) eigs = eigvals(rdm) - @test all(>=(0), real(eigs)) && all(==(0), imag(eigs)) @test size(rdm) == (2^length(vs), 2^length(vs)) + @test all(>=(0), real(eigs)) && all(==(0), imag(eigs)) - #Test more advanced block BP with MPS message tensors on a grid - g_dims = (5, 4) + # #Test more advanced block BP with MPS message tensors on a grid + g_dims = (4, 3) g = named_grid(g_dims) s = siteinds("S=1/2", g) χ = 2 ψ = randomITensorNetwork(s; link_space=χ) - maxdim = 16 v = (2, 2) ψψ = flatten_networks(ψ, dag(ψ); combine_linkinds=false, map_bra_linkinds=prime) @@ -156,24 +130,19 @@ ITensors.disable_warn_order() combiners = linkinds_combiners(ψψ) ψψ = combine_linkinds(ψψ, combiners) ψOψ = combine_linkinds(ψOψ, combiners) - - Z = partition(ψψ, group(v -> v[1], vertices(ψψ))) - mts = message_tensors(Z) + pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) mts = belief_propagation( - ψψ, - mts; - contract_kwargs=(; - alg="density_matrix", output_structure=path_graph_structure, cutoff=1e-16, maxdim - ), + pψψ; + contractor=ITensorNetworks.contract_density_matrix, + contractor_kwargs=(; cutoff=1e-6, maxdim=4), ) - numerator_network = approx_network_region(ψψ, mts, [v]; verts_tn=ITensorNetwork(ψOψ[v])) - - denominator_network = approx_network_region(ψψ, mts, [v]) - bp_sz = ITensors.contract(numerator_network)[] / ITensors.contract(denominator_network)[] + env_tensors = environment_tensors(pψψ, mts, [v]) + numerator = contract(vcat(env_tensors, ITensor[ψOψ[v]]))[] + denominator = contract(vcat(env_tensors, ITensor[ψψ[v]]))[] exact_sz = contract_boundary_mps(ψOψ; cutoff=1e-16) / contract_boundary_mps(ψψ; cutoff=1e-16) - @test abs.(bp_sz - exact_sz) <= 1e-5 + @test abs.((numerator / denominator) - exact_sz) <= 1e-5 end diff --git a/test/test_binary_tree_partition.jl b/test/test_binary_tree_partition.jl index 32942b4a..4fe55250 100644 --- a/test/test_binary_tree_partition.jl +++ b/test/test_binary_tree_partition.jl @@ -1,5 +1,6 @@ using ITensors, OMEinsumContractionOrders using Graphs, NamedGraphs +using ITensorNetworks using ITensors: contract using ITensorNetworks: _root, @@ -8,7 +9,9 @@ using ITensorNetworks: _is_rooted_directed_binary_tree, _contract_deltas_ignore_leaf_partitions, _rem_vertex!, - _DensityMartrixAlgGraph + _DensityMartrixAlgGraph, + _partition +using Test @testset "test mincut functions on top of MPS" begin i = Index(2, "i") @@ -61,7 +64,7 @@ end @testset "test partition with mincut_recursive_bisection alg of disconnected tn" begin inds = [Index(2, "$i") for i in 1:5] tn = ITensorNetwork([randomITensor(i) for i in inds]) - par = partition(tn, binary_tree_structure(tn); alg="mincut_recursive_bisection") + par = _partition(tn, binary_tree_structure(tn); alg="mincut_recursive_bisection") networks = [Vector{ITensor}(par[v]) for v in vertices(par)] network = vcat(networks...) @test isapprox(contract(Vector{ITensor}(tn)), contract(network...)) @@ -80,7 +83,7 @@ end out1 = contract(network...) tn = ITensorNetwork(network) inds_btree = binary_tree_structure(tn) - par = partition(tn, inds_btree; alg="mincut_recursive_bisection") + par = _partition(tn, inds_btree; alg="mincut_recursive_bisection") par = _contract_deltas_ignore_leaf_partitions(par; root=_root(inds_btree)) networks = [Vector{ITensor}(par[v]) for v in vertices(par)] network2 = vcat(networks...) @@ -113,7 +116,7 @@ end M = MPS(T, (i, j, k, l, m); cutoff=1e-5, maxdim=5) tn = ITensorNetwork(M[:]) out_tree = path_graph_structure(tn) - input_partition = partition(tn, out_tree; alg="mincut_recursive_bisection") + input_partition = _partition(tn, out_tree; alg="mincut_recursive_bisection") underlying_tree = underlying_graph(input_partition) # Change type of each partition[v] since they will be updated # with potential data type chage. diff --git a/test/test_contract_deltas.jl b/test/test_contract_deltas.jl index d48eb2d5..151aa6c3 100644 --- a/test/test_contract_deltas.jl +++ b/test/test_contract_deltas.jl @@ -32,7 +32,7 @@ end end tn = ITensorNetwork(vec(tn[:, :, 1])) for inds_tree in [binary_tree_structure(tn), path_graph_structure(tn)] - par = partition(tn, inds_tree; alg="mincut_recursive_bisection") + par = _partition(tn, inds_tree; alg="mincut_recursive_bisection") root = _root(inds_tree) par_contract_deltas = _contract_deltas_ignore_leaf_partitions(par; root=root) @test Set(_noncommoninds(par)) == Set(_noncommoninds(par_contract_deltas)) diff --git a/test/test_examples/test_apply_bp.jl b/test/test_examples/test_apply_bp.jl deleted file mode 100644 index 6558fbe5..00000000 --- a/test/test_examples/test_apply_bp.jl +++ /dev/null @@ -1,28 +0,0 @@ -using ITensorNetworks -using Suppressor -using Test - -include(joinpath(pkgdir(ITensorNetworks), "examples", "apply", "apply_bp", "apply_bp.jl")) - -@testset "Test apply_bp example" begin - opnames = ["Id", "RandomUnitary"] - graphs = (named_comb_tree, named_grid) - dims = (6, 6) - @testset "$opname, $graph" for opname in opnames, graph in graphs - ψ_bp, mts_bp, ψ_vidal, mts_vidal = @suppress main(; - seed=1234, - opname, - graph, - dims, - χ=2, - nlayers=2, - variational_optimization_only=false, - regauge=false, - reduced=true, - ) - v = dims .÷ 2 - sz_bp = expect_bp("Sz", v, ψ_bp, mts_bp) - sz_vidal = expect_bp("Sz", v, ψ_vidal, mts_vidal) - @test sz_bp ≈ sz_vidal rtol = 1e-5 - end -end diff --git a/test/test_examples/test_examples.jl b/test/test_examples/test_examples.jl index 8afc3be4..b805db50 100644 --- a/test/test_examples/test_examples.jl +++ b/test/test_examples/test_examples.jl @@ -8,15 +8,11 @@ using Test "boundary.jl", "distances.jl", "examples.jl", - "group_partition.jl", "mincut.jl", "mps.jl", "peps.jl", "steiner_tree.jl", - joinpath("belief_propagation", "bpexample.jl"), - joinpath("belief_propagation", "bpsequences.jl"), joinpath("dynamics", "2d_ising_imag_tebd.jl"), - joinpath("dynamics", "heavy_hex_ising_real_tebd.jl"), joinpath("treetensornetworks", "comb_tree.jl"), joinpath("treetensornetworks", "spanning_tree.jl"), joinpath("treetensornetworks", "ttn_basics.jl"), @@ -26,11 +22,7 @@ using Test @suppress include(joinpath(pkgdir(ITensorNetworks), "examples", example_file)) end if !Sys.iswindows() - example_files = [ - joinpath("contraction_sequence", "contraction_sequence.jl"), - joinpath("partition", "kahypar_vs_metis.jl"), - joinpath("partition", "partitioning.jl"), - ] + example_files = [joinpath("contraction_sequence", "contraction_sequence.jl")] @testset "Test $example_file (using KaHyPar, so no Windows support)" for example_file in example_files @suppress include(joinpath(pkgdir(ITensorNetworks), "examples", example_file)) diff --git a/test/test_examples/test_sqrt_bp.jl b/test/test_examples/test_sqrt_bp.jl deleted file mode 100644 index 19d90aa5..00000000 --- a/test/test_examples/test_sqrt_bp.jl +++ /dev/null @@ -1,10 +0,0 @@ -using ITensorNetworks -using Suppressor -using Test - -include(joinpath(pkgdir(ITensorNetworks), "examples", "belief_propagation", "sqrt_bp.jl")) - -@testset "Test sqrt_bp example" begin - (; sz_bp, sz_sqrt_bp) = @suppress main(; n=8, niters=10, β=0.28, h=0.5) - @test sz_bp ≈ sz_sqrt_bp -end diff --git a/test/test_gauging.jl b/test/test_gauging.jl index 5d5e6c77..6441b766 100644 --- a/test/test_gauging.jl +++ b/test/test_gauging.jl @@ -5,11 +5,9 @@ using ITensorNetworks: contract_inner, symmetric_gauge, symmetric_to_vidal_gauge, - message_tensors, vidal_itn_canonicalness, vidal_gauge, - symmetric_itn_canonicalness, - belief_propagation_iteration + symmetric_itn_canonicalness using NamedGraphs using Test using Compat @@ -25,32 +23,26 @@ using SplitApplyCombine Random.seed!(5467) ψ = randomITensorNetwork(s; link_space=χ) - ψ_symm, ψ_symm_mts = symmetric_gauge(ψ; niters=50) + ψ_symm, pψψ_symm, ψ_symm_mts = symmetric_gauge(ψ; niters=50) + + @test symmetric_itn_canonicalness(ψ_symm, pψψ_symm, ψ_symm_mts) < 1e-5 #Test we just did a gauge transform and didn't change the overall network @test contract_inner(ψ_symm, ψ) / sqrt(contract_inner(ψ_symm, ψ_symm) * contract_inner(ψ, ψ)) ≈ 1.0 - ψψ_symm = ψ_symm ⊗ prime(dag(ψ_symm); sites=[]) - Z = partition( - ψψ_symm; subgraph_vertices=collect(values(group(v -> v[1], vertices(ψψ_symm)))) - ) - ψ_symm_mts_V2 = message_tensors(Z) - ψ_symm_mts_V2 = belief_propagation( - ψψ_symm, ψ_symm_mts_V2; contract_kwargs=(; alg="exact"), niters=50 - ) + ψψ_symm_V2 = ψ_symm ⊗ prime(dag(ψ_symm); sites=[]) + pψψ_symm_V2 = PartitionedGraph(ψψ_symm_V2, group(v -> v[1], vertices(ψψ_symm_V2))) + ψ_symm_mts_V2 = belief_propagation(pψψ_symm_V2; niters=50) - for e in edges(ψ_symm_mts_V2) + for m_e in values(ψ_symm_mts_V2) #Test all message tensors are approximately diagonal - m_e = ψ_symm_mts_V2[e][first(vertices(ψ_symm_mts_V2[e]))] - @test diagITensor(vector(diag(m_e)), inds(m_e)) ≈ m_e atol = 1e-8 + @test diagITensor(vector(diag(only(m_e))), inds(only(m_e))) ≈ only(m_e) atol = 1e-8 end - @test symmetric_itn_canonicalness(ψ_symm, ψ_symm_mts) < 1e-5 - ψ_vidal, bond_tensors = vidal_gauge(ψ; target_canonicalness=1e-6) @test vidal_itn_canonicalness(ψ_vidal, bond_tensors) < 1e-5 - ψ_vidal, bond_tensors = symmetric_to_vidal_gauge(ψ_symm, ψ_symm_mts) + ψ_vidal, bond_tensors = symmetric_to_vidal_gauge(ψ_symm, pψψ_symm, ψ_symm_mts) @test vidal_itn_canonicalness(ψ_vidal, bond_tensors) < 1e-5 end