diff --git a/examples/apply/apply_bp/apply_bp.jl b/examples/apply/apply_bp/apply_bp.jl deleted file mode 100644 index 4fd6159c..00000000 --- a/examples/apply/apply_bp/apply_bp.jl +++ /dev/null @@ -1,169 +0,0 @@ -using ITensorNetworks -using ITensorNetworks: - approx_network_region, - belief_propagation, - get_environment, - contract_inner, - message_tensors, - neighbor_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, ψ, pψψ, mts) - s = siteinds(ψ) - numerator_tensors = approx_network_region( - pψψ, mts, [(v, 1)]; verts_tensors=ITensor[apply(op(opname, s[v]), ψ[v])] - ) - denominator_tensors = approx_network_region(pψψ, mts, [(v, 1)]) - return contract(numerator_tensors)[] / contract(denominator_tensors)[] -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(ψ) - pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) - mts = belief_propagation( - pψψ; contract_kwargs=(; alg="exact"), niters=50, target_precision=1e-5 - ) - edges = PartitionEdge.(NamedGraphs.edges(partitioned_graph(pψψ))) - for layer in eachindex(os) - @show layer - o⃗ = os[layer] - for o in o⃗ - v⃗ = neighbor_vertices(ψ, o) - for e in edges - @assert order(only(mts[e])) == 2 - @assert order(only(mts[reverse(e)])) == 2 - end - - @assert length(v⃗) == 2 - v1, v2 = v⃗ - - pe = partitionedge(pψψ, NamedEdge((v1, 1) => (v2, 1))) - envs = get_environment(pψψ, 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(ψ) - pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) - mts[pe] = dense.(obs.singular_values) - mts[reverse(pe)] = dense.(obs.singular_values) - end - if regauge - println("regauge") - mts = belief_propagation( - pψψ, mts; contract_kwargs=(; alg="exact"), niters=50, target_precision=1e-5 - ) - end - end - return ψ, pψψ, mts -end - -function simple_update_vidal(os, ψ::ITensorNetwork; maxdim, regauge=false) - println("Simple update, Vidal gauge") - ψψ = norm_network(ψ) - pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) - mts = belief_propagation( - pψψ; contract_kwargs=(; alg="exact"), niters=50, target_precision=1e-5 - ) - ψ, bond_tensors = vidal_gauge(ψ, pψψ, 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, pψψ_symmetric, mts = vidal_to_symmetric_gauge(ψ, bond_tensors) - mts = belief_propagation( - pψψ_symmetric, - mts; - contract_kwargs=(; alg="exact"), - niters=50, - target_precision=1e-5, - ) - ψ, bond_tensors = vidal_gauge(ψ_symmetric, pψψ_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, pψψ_bp, mts_bp = 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, pψψ_vidal, mts_vidal = vidal_to_symmetric_gauge(ψ_vidal, bond_tensors) - - return ψ_bp, pψψ_bp, mts_bp, ψ_vidal, pψψ_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 7c6ff8df..00000000 --- a/examples/apply/apply_bp/apply_bp_run.jl +++ /dev/null @@ -1,74 +0,0 @@ -include("apply_bp.jl") - -# opname = "Id" -opname = "RandomUnitary" - -@show opname - -# graph = named_comb_tree -graph = named_grid - -dims = (6, 6) - -ψ_bp, pψψ_bp, mts_bp, ψ_vidal, pψψ_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, pψψ_bp, mts_bp) -sz_vidal = @show expect_bp("Sz", v, ψ_vidal, pψψ_vidal, mts_vidal) -@show abs(sz_bp - sz_vidal) / abs(sz_vidal) - -# Run BP again -mts_bp = belief_propagation( - pψψ_bp, - mts_bp; - contract_kwargs=(; alg="exact"), - niters=50, - target_precision=1e-7, - verbose=true, -) -mts_vidal = belief_propagation( - pψψ_vidal, - mts_vidal; - contract_kwargs=(; alg="exact"), - niters=50, - target_precision=1e-7, - verbose=true, -) - -sz_bp = @show expect_bp("Sz", v, ψ_bp, pψψ_bp, mts_bp) -sz_vidal = @show expect_bp("Sz", v, ψ_vidal, pψψ_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 62eba8fa..00000000 --- a/examples/belief_propagation/bpexample.jl +++ /dev/null @@ -1,97 +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 - -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 - pψψ = PartitionedGraph(ψψ, collect(values(group(v -> v[1], vertices(ψψ))))) - mts = belief_propagation( - pψψ; contract_kwargs=(; alg="exact"), verbose=true, niters=10, target_precision=1e-3 - ) - numerator_tensors = approx_network_region( - pψψ, mts, [(v, 1)]; verts_tensors=[apply(op("Sz", s[v]), ψ[v])] - ) - denominator_tensors = approx_network_region(pψψ, mts, [(v, 1)]) - sz_bp = contract(numerator_tensors)[] / contract(denominator_tensors)[] - - println( - "Simple Belief Propagation Gives Sz on Site " * string(v) * " as " * string(sz_bp) - ) - - #Now do Column-wise General Belief Propagation to Measure Sz on Site v - pψψ = PartitionedGraph(ψψ, collect(values(group(v -> v[1][1], vertices(ψψ))))) - mts = belief_propagation( - pψψ; contract_kwargs=(; alg="exact"), verbose=true, niters=10, target_precision=1e-3 - ) - numerator_tensors = approx_network_region( - pψψ, mts, [(v, 1)]; verts_tensors=[apply(op("Sz", s[v]), ψ[v])] - ) - denominator_tensors = approx_network_region(pψψ, mts, [(v, 1)]) - sz_gen_bp = contract(numerator_tensors)[] / contract(denominator_tensors)[] - - println( - "General Belief Propagation Gives Sz on Site " * string(v) * " as " * string(sz_gen_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) - - pψψ = PartitionedGraph(ψψ, collect(values(group(v -> v[1], vertices(ψψ))))) - mts = belief_propagation( - pψψ; - itensor_constructor=inds_e -> ITensor[dense(delta(i)) for i in inds_e], - contract_kwargs=(; - alg="density_matrix", - output_structure=path_graph_structure, - maxdim=8, - contraction_sequence_alg="optimal", - ), - ) - numerator_tensors = approx_network_region(pψψ, mts, [v]; verts_tensors=[ψOψ[v]]) - denominator_tensors = approx_network_region(pψψ, mts, [v]) - sz_MPS_bp = contract(numerator_tensors)[] / contract(denominator_tensors)[] - - println( - "Column-Wise MPS Belief Propagation Gives Sz on Site " * - string(v) * - " as " * - string(sz_gen_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 8f3b1438..00000000 --- a/examples/belief_propagation/bpsequences.jl +++ /dev/null @@ -1,80 +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, 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 - pψψ = PartitionedGraph(ψψ, collect(values(group(v -> v[1], vertices(ψψ))))) - mts_init = message_tensors(pψψ) - - 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( - pψψ, - mts_init; - contract_kwargs=(; alg="exact"), - target_precision=1e-10, - niters=100, - edges=[ - PartitionEdge.(e) for e in edge_sequence(partitioned_graph(pψψ); alg="parallel") - ], - verbose=true, - ) - print("Sequential updates (sequence is default edge list of the message tensors): ") - belief_propagation( - pψψ, - mts_init; - contract_kwargs=(; alg="exact"), - target_precision=1e-10, - niters=100, - edges=PartitionEdge.([ - e for - e in vcat(edges(partitioned_graph(pψψ)), reverse.(edges(partitioned_graph(pψψ)))) - ]), - verbose=true, - ) - print("Sequential updates (sequence is our custom sequence finder): ") - belief_propagation( - pψψ, - 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 512e5981..00000000 --- a/examples/belief_propagation/sqrt_bp.jl +++ /dev/null @@ -1,71 +0,0 @@ -using ITensors -using ITensorNetworks -using Random -using SplitApplyCombine -using NamedGraphs - -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) - 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) - - #Now do Simple Belief Propagation to Measure Sz on Site v - pψψ = PartitionedGraph(ψψ, collect(values(group(v -> v[1], vertices(ψψ))))) - - mts = @time belief_propagation(pψψ; niters, contract_kwargs=(; alg="exact")) - numerator_network = approx_network_region( - pψψ, mts, [(v, 1)]; verts_tensors=ITensor[apply(op("Sz", s[v]), ψ[v])] - ) - denominator_network = approx_network_region(pψψ, mts, [(v, 1)]) - sz_bp = - ITensors.contract( - numerator_network; sequence=contraction_sequence(numerator_network; alg="optimal") - )[] / ITensors.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 = @time sqrt_belief_propagation(pψψ; niters) - - numerator_network = approx_network_region( - pψψ, mts_sqrt, [(v, 1)]; verts_tensors=ITensor[apply(op("Sz", s[v]), ψ[v])] - ) - denominator_network = approx_network_region(pψψ, 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 71071119..6305744d 100644 --- a/examples/boundary.jl +++ b/examples/boundary.jl @@ -8,10 +8,10 @@ tn = ITensorNetwork(named_grid((6, 3)); link_space=4) @visualize tn -ptn = PartitionedGraph( - tn, partitioned_vertices(underlying_graph(tn); nvertices_per_partition=2) -) -sub_vs_1, sub_vs_2 = vertices(ptn, PartitionVertex(1)), vertices(ptn, PartitionVertex(2)) +pvertices = partitioned_vertices(underlying_graph(tn); nvertices_per_partition=2) +ptn = PartitionedGraph(tn, pvertices) +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/dynamics/heavy_hex_ising_real_tebd.jl b/examples/dynamics/heavy_hex_ising_real_tebd.jl deleted file mode 100644 index d17084a9..00000000 --- a/examples/dynamics/heavy_hex_ising_real_tebd.jl +++ /dev/null @@ -1,131 +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, pψψ::PartitionedGraph, mts) - 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_tensors = approx_network_region( - pψψ, mts, vs_braket; verts_tensors=ITensor[Oψ[v] for v in vs] - ) - denominator_tensors = approx_network_region(pψψ, mts, vs_braket) - num_seq = contraction_sequence(numerator_tensors; alg="optimal") - den_seq = contraction_sequence(denominator_tensors; alg="optimal") - return contract(numerator_tensors; sequence=num_seq)[] / - contract(denominator_tensors; 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 - ψ, pψψ, mts = vidal_to_symmetric_gauge(ψ, bond_tensors) - mag_dict = Dict( - zip( - [v for v in vertices(ψ)], - [expect_state_SBP(op("Z", s[v]), ψ, pψψ, 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 6c988753..00000000 --- a/examples/gauging/gauging_itns.jl +++ /dev/null @@ -1,148 +0,0 @@ -using Compat -using ITensors -using Metis -using ITensorNetworks -using Random -using SplitApplyCombine - -using ITensorNetworks: - message_tensors, - belief_propagation, - vidal_gauge, - symmetric_gauge, - vidal_itn_canonicalness, - initialize_bond_tensors, - vidal_itn_isometries, - norm_network, - edge_sequence, - belief_propagation_iteration - -using NamedGraphs -using NamedGraphs: add_edges!, rem_vertex!, hexagonal_lattice_graph -using Graphs - -"""Eager Gauging""" -function eager_gauging( - ψ::ITensorNetwork, pψψ::PartitionedGraph, bond_tensors::DataGraph, mts -) - isometries = vidal_itn_isometries(ψ, bond_tensors) - - ψ = copy(ψ) - mts = copy(mts) - for e in edges(ψ) - vsrc, vdst = src(e), dst(e) - pe = which_partitionedge(pψψ, NamedEdge((vsrc, 1) => (vdst, 1))) - normalize!(isometries[e]) - normalize!(isometries[reverse(e)]) - mts[pe], mts[reverse(pe)] = ITensorNetwork(isometries[e]), - ITensorNetwork(isometries[reverse(e)]) - end - - ψ, bond_tensors = vidal_gauge(ψ, pψψ, 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(ψ) - pψψ = PartitionedGraph(ψψ, collect(values(group(v -> v[1], vertices(ψψ))))) - mts = message_tensors(pψψ) - bond_tensors = initialize_bond_tensors(ψ) - - 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( - pψψ, mts; contract_kwargs=(; alg="exact") - ) - else - times_iters[i] = @elapsed mts, _ = belief_propagation_iteration( - pψψ, - mts; - contract_kwargs=(; alg="exact"), - edges=[ - PartitionEdge.(e) for e in edge_sequence(partitioned_graph(pψψ); alg="parallel") - ], - ) - end - - times_gauging[i] = @elapsed ψ, bond_tensors = vidal_gauge(ψinit, pψψ, mts) - elseif mode == "eager" - times_iters[i] = @elapsed ψ, bond_tensors, mts = eager_gauging( - ψ, pψψ, 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/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 47f14b4e..b457864f 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -99,7 +99,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") diff --git a/src/apply.jl b/src/apply.jl index 676f5867..0923229a 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -297,10 +297,10 @@ 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; + bond_tensors; normalize=false, apply_kwargs..., ) @@ -314,13 +314,13 @@ function ITensors.apply( for vn in neighbors(ψ, src(e)) if (vn != dst(e)) - ψv1 = noprime(ψv1 * bond_tensors[vn => src(e)]) + ψv1 = noprime(ψv1 * bond_tensors[NamedEdge(vn => src(e))]) end end for vn in neighbors(ψ, dst(e)) if (vn != src(e)) - ψv2 = noprime(ψv2 * bond_tensors[vn => dst(e)]) + ψv2 = noprime(ψv2 * bond_tensors[NamedEdge(vn => dst(e))]) end end @@ -339,17 +339,17 @@ function ITensors.apply( replaceind!(S, ind_to_replace, ind_to_replace_with') replaceind!(V, ind_to_replace, ind_to_replace_with) - ψv1, bond_tensors[e], ψv2 = U * Qᵥ₁, S, V * Qᵥ₂ + ψv1, bond_tensors[e], bond_tensors[reverse(e)], ψv2 = U * Qᵥ₁, S, S, V * Qᵥ₂ for vn in neighbors(ψ, src(e)) if (vn != dst(e)) - ψv1 = noprime(ψv1 * inv_diag(bond_tensors[vn => src(e)])) + ψv1 = noprime(ψv1 * inv_diag(bond_tensors[NamedEdge(vn => src(e))])) end end for vn in neighbors(ψ, dst(e)) if (vn != src(e)) - ψv2 = noprime(ψv2 * inv_diag(bond_tensors[vn => dst(e)])) + ψv2 = noprime(ψv2 * inv_diag(bond_tensors[NamedEdge(vn => dst(e))])) end end diff --git a/src/beliefpropagation/beliefpropagation.jl b/src/beliefpropagation/beliefpropagation.jl index 5a0b7b79..cfbafdf6 100644 --- a/src/beliefpropagation/beliefpropagation.jl +++ b/src/beliefpropagation/beliefpropagation.jl @@ -3,8 +3,8 @@ function message_tensors( ) mts = Dict() for e in partitionedges(ptn) - src_e_itn = unpartitioned_graph(subgraph(ptn, [src(e)])) - dst_e_itn = unpartitioned_graph(subgraph(ptn, [dst(e)])) + src_e_itn = subgraph(ptn, src(e)) + dst_e_itn = subgraph(ptn, dst(e)) inds_e = commoninds(src_e_itn, dst_e_itn) mts[e] = itensor_constructor(inds_e) mts[reverse(e)] = dag.(mts[e]) @@ -16,10 +16,7 @@ end Do a single update of a message tensor using the current subgraph and the incoming mts """ function update_message_tensor( - ptn::PartitionedGraph, - edge::PartitionEdge, - mts; - contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1), + ptn::PartitionedGraph, edge::PartitionEdge, mts; contract_kwargs=(; alg="exact") ) pedges = setdiff( partitionedges(ptn, boundary_edges(ptn, vertices(ptn, src(edge)); dir=:in)), @@ -30,7 +27,7 @@ function update_message_tensor( contract_list = ITensor[ incoming_messages - ITensor(unpartitioned_graph(subgraph(ptn, [src(edge)]))) + ITensor(subgraph(ptn, src(edge))) ] if contract_kwargs.alg != "exact" @@ -54,7 +51,7 @@ function belief_propagation_iteration( ptn::PartitionedGraph, mts, edges::Vector{<:PartitionEdge}; - contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1), + contract_kwargs=(; alg="exact"), compute_norm=false, ) new_mts = copy(mts) @@ -82,7 +79,7 @@ function belief_propagation_iteration( ptn::PartitionedGraph, mts, edge_groups::Vector{<:Vector{<:PartitionEdge}}; - contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1), + contract_kwargs=(; alg="exact"), compute_norm=false, ) new_mts = copy(mts) @@ -102,7 +99,7 @@ end function belief_propagation_iteration( ptn::PartitionedGraph, mts; - contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1), + contract_kwargs=(; alg="exact"), compute_norm=false, edges=PartitionEdge.(edge_sequence(partitioned_graph(ptn))), ) @@ -112,7 +109,7 @@ end function belief_propagation( ptn::PartitionedGraph, mts; - contract_kwargs=(; alg="density_matrix", output_structure=path_graph_structure, maxdim=1), + contract_kwargs=(; alg="exact"), niters=default_bp_niters(partitioned_graph(ptn)), target_precision=nothing, edges=PartitionEdge.(edge_sequence(partitioned_graph(ptn))), @@ -143,17 +140,12 @@ function belief_propagation( 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(ptn::PartitionedGraph, mts, verts::Vector; dir=:in) +function get_environment_tensors(ptn::PartitionedGraph, mts, verts::Vector) partition_verts = partitionvertices(ptn, verts) central_verts = vertices(ptn, partition_verts) - if dir == :out - return get_environment(ptn, mts, setdiff(vertices(ptn), verts)) - end - 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[]) @@ -164,6 +156,12 @@ function get_environment(ptn::PartitionedGraph, mts, verts::Vector; dir=:in) return vcat(env_tensors, central_tensors) end +function get_environment_tensors( + ptn::PartitionedGraph, mts, partition_verts::Vector{<:PartitionVertex} +) + return get_environment_tensors(ptn, mts, vertices(ptn, partition_verts)) +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 @@ -174,7 +172,18 @@ function approx_network_region( verts::Vector; verts_tensors=ITensor[(unpartitioned_graph(ptn))[v] for v in verts], ) - environment_tensors = get_environment(ptn, mts, verts) + environment_tensors = get_environment_tensors(ptn, mts, verts) return vcat(environment_tensors, verts_tensors) end + +function approx_network_region( + ptn::PartitionedGraph, + mts, + partition_verts::Vector{<:PartitionVertex}; + verts_tensors=ITensor[ + (unpartitioned_graph(ptn))[v] for v in vertices(ptn, partition_verts) + ], +) + return approx_network_region(ptn, mts, vertices(ptn, partition_verts); verts_tensors) +end diff --git a/src/beliefpropagation/sqrt_beliefpropagation.jl b/src/beliefpropagation/sqrt_beliefpropagation.jl deleted file mode 100644 index dc13713e..00000000 --- a/src/beliefpropagation/sqrt_beliefpropagation.jl +++ /dev/null @@ -1,141 +0,0 @@ -# using ITensors: scalartype -# using ITensorNetworks: find_subgraph, map_diag, sqrt_diag, boundary_edges - -function sqrt_belief_propagation_iteration( - pψψ::PartitionedGraph, sqrt_mts, edges::Vector{<:PartitionEdge} -) - new_sqrt_mts = copy(sqrt_mts) - c = 0.0 - for e in edges - new_sqrt_mts[e] = ITensor[update_sqrt_message_tensor(pψψ, e, new_sqrt_mts;)] - - # 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( - pψψ::PartitionedGraph, sqrt_mts, edges::Vector{<:Vector{<:PartitionEdge}} -) - new_sqrt_mts = copy(sqrt_mts) - c = 0.0 - for e_group in edges - updated_sqrt_mts, ct = sqrt_belief_propagation_iteration(pψψ, 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( - pψψ::PartitionedGraph, sqrt_mts; edges=edge_sequence(partitioned_graph(pψψ)) -) - return sqrt_belief_propagation_iteration(pψψ, sqrt_mts, edges) -end - -function sqrt_belief_propagation( - pψψ::PartitionedGraph, - mts; - niters=default_bp_niters(partitioned_graph(pψψ)), - edges=PartitionEdge.(edge_sequence(partitioned_graph(pψψ))), - # target_precision::Union{Float64,Nothing}=nothing, -) - # compute_norm = target_precision == nothing ? false : true - sqrt_mts = sqrt_message_tensors(pψψ, 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(pψψ, 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 sqrt_belief_propagation(ptn::PartitionedGraph; kwargs...) - mts = message_tensors(ptn) - return sqrt_belief_propagation(ptn, mts; kwargs...) -end - -function update_sqrt_message_tensor(pψψ::PartitionedGraph, edge::PartitionEdge, sqrt_mts;) - v = only(filter(v -> v[2] == 1, vertices(pψψ, src(edge)))) - site_itensor = unpartitioned_graph(pψψ)[v] - p_edges = setdiff(partitionedges(pψψ, boundary_edges(pψψ, [v]; dir=:in)), [reverse(edge)]) - incoming_messages = [sqrt_mts[e_in] for e_in in p_edges] - incoming_messages = reduce(vcat, incoming_messages; init=ITensor[]) - contract_list = ITensor[incoming_messages; site_itensor] - contract_output = contract( - contract_list; sequence=contraction_sequence(contract_list; alg="optimal") - ) - left_inds = [ - uniqueinds(contract_output, site_itensor) - siteinds(unpartitioned_graph(pψψ), v) - ] - Q, R = qr(contract_output, left_inds) - normalize!(R) - return R -end - -function sqrt_message_tensors( - pψψ::PartitionedGraph, - mts; - eigen_message_tensor_cutoff=10 * eps(real(scalartype(unpartitioned_graph(pψψ)))), - regularization=10 * eps(real(scalartype(unpartitioned_graph(pψψ)))), -) - sqrt_mts = copy(mts) - for e in partitionedges(pψψ) - vsrc, vdst = filter(v -> v[2] == 1, vertices(pψψ, src(e))), - filter(v -> v[2] == 1, vertices(pψψ, dst(e))) - ψvsrc, ψvdst = unpartitioned_graph(pψψ)[only(vsrc)], - unpartitioned_graph(pψψ)[only(vdst)] - - edge_ind = commoninds(ψvsrc, ψvdst) - edge_ind_sim = sim(edge_ind) - - X_D, X_U = eigen(only(mts[e]); ishermitian=true, cutoff=eigen_message_tensor_cutoff) - Y_D, Y_U = eigen( - only(mts[reverse(e)]); 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[e] = ITensor[rootX] - sqrt_mts[reverse(e)] = ITensor[rootY] - end - return sqrt_mts -end - -function sqr_message_tensors(sqrt_mts) - mts = copy(sqrt_mts) - for e in keys(sqrt_mts) - sqrt_mt = only(sqrt_mts[e]) - sqrt_mt_rev = only(sqrt_mts[reverse(e)]) - 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] = ITensor[mt] - mts[reverse(e)] = ITensor[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/exports.jl b/src/exports.jl index 24623f1c..f298994d 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -73,7 +73,6 @@ export AbstractITensorNetwork, TreeTensorNetwork, TTN, random_ttn, - PartitionedITensorNetwork, ProjTTN, ProjTTNSum, ProjTTNApply, diff --git a/src/gauging.jl b/src/gauging.jl index 25fce5d0..4b36badf 100644 --- a/src/gauging.jl +++ b/src/gauging.jl @@ -1,10 +1,11 @@ """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 = Dict() for e in edges(ψ) index = commoninds(ψ[src(e)], ψ[dst(e)]) bond_tensors[e] = dense(delta(index, index_map(index))) + bond_tensors[reverse(e)] = bond_tensors[e] end return bond_tensors @@ -15,7 +16,7 @@ function vidal_gauge( ψ::ITensorNetwork, pψψ::PartitionedGraph, mts, - bond_tensors::DataGraph; + bond_tensors; eigen_message_tensor_cutoff=10 * eps(real(scalartype(ψ))), regularization=10 * eps(real(scalartype(ψ))), edges=NamedGraphs.edges(ψ), @@ -67,7 +68,7 @@ function vidal_gauge( [commoninds(S, U)..., commoninds(S, V)...] => [new_edge_ind..., prime(new_edge_ind)...], ) - bond_tensors[e] = S + bond_tensors[e], bond_tensors[reverse(e)] = S, S end return ψ_vidal, bond_tensors @@ -107,7 +108,7 @@ function vidal_gauge( svd_kwargs..., ) ψψ = norm_network(ψ) - pψψ = PartitionedGraph(ψψ, collect(values(group(v -> v[1], vertices(ψψ))))) + pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) mts = message_tensors(pψψ) mts = belief_propagation( @@ -124,10 +125,10 @@ function vidal_gauge( end """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) +function vidal_to_symmetric_gauge(ψ::ITensorNetwork, bond_tensors) ψsymm = copy(ψ) ψψsymm = norm_network(ψsymm) - pψψsymm = PartitionedGraph(ψψsymm, collect(values(group(v -> v[1], vertices(ψψsymm))))) + pψψsymm = PartitionedGraph(ψψsymm, group(v -> v[1], vertices(ψψsymm))) ψsymm_mts = message_tensors(pψψsymm) for e in edges(ψsymm) @@ -142,7 +143,7 @@ function vidal_to_symmetric_gauge(ψ::ITensorNetwork, bond_tensors::DataGraph) end ψψsymm = norm_network(ψsymm) - pψψsymm = PartitionedGraph(ψψsymm, collect(values(group(v -> v[1], vertices(ψψsymm))))) + pψψsymm = PartitionedGraph(ψψsymm, group(v -> v[1], vertices(ψψsymm))) return ψsymm, pψψsymm, ψsymm_mts end @@ -175,14 +176,14 @@ function symmetric_to_vidal_gauge( mts; regularization=10 * eps(real(scalartype(ψ))), ) - bond_tensors = DataGraph{vertextype(ψ),ITensor,ITensor}(underlying_graph(ψ)) + bond_tensors = Dict() ψ_vidal = copy(ψ) for e in edges(ψ) vsrc, vdst = src(e), dst(e) pe = partitionedge(pψψ, NamedEdge((vsrc, 1) => (vdst, 1))) - bond_tensors[e] = only(mts[pe]) + 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) @@ -194,16 +195,16 @@ end """Function to measure the 'isometries' of a state in the Vidal Gauge""" function vidal_itn_isometries( ψ::ITensorNetwork, - bond_tensors::DataGraph; + bond_tensors; 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) ψv = copy(ψ[vsrc]) for vn in setdiff(neighbors(ψ, vsrc), [vdst]) - ψv = noprime(ψv * bond_tensors[vn => vsrc]) + ψv = noprime(ψv * bond_tensors[NamedEdge(vn => vsrc)]) end ψvdag = dag(ψv) @@ -215,19 +216,19 @@ function vidal_itn_isometries( end """Function to measure the 'canonicalness' of a state in the Vidal Gauge""" -function vidal_itn_canonicalness(ψ::ITensorNetwork, bond_tensors::DataGraph) +function vidal_itn_canonicalness(ψ::ITensorNetwork, bond_tensors) f = 0 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""" diff --git a/test/test_apply.jl b/test/test_apply.jl index 4cdc6d8b..55bca9aa 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,10 +1,11 @@ using ITensorNetworks using ITensorNetworks: belief_propagation, - get_environment, + get_environment_tensors, contract_inner, message_tensors, vidal_gauge, + vidal_apply, vidal_to_symmetric_gauge, norm_network using Test @@ -29,16 +30,16 @@ using SplitApplyCombine ψψ = norm_network(ψ) #Simple Belief Propagation Grouping - pψψ_SBP = PartitionedGraph(ψψ, collect(values(group(v -> v[1], vertices(ψψ))))) + pψψ_SBP = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) mtsSBP = belief_propagation(pψψ_SBP; contract_kwargs=(; alg="exact"), niters=50) - envsSBP = get_environment(pψψ_SBP, mtsSBP, [(v1, 1), (v1, 2), (v2, 1), (v2, 2)]) + envsSBP = get_environment_tensors(pψψ_SBP, mtsSBP, PartitionVertex.([v1, v2])) ψ_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) - pψψ_GBP = PartitionedGraph(ψψ, collect(values(group(v -> v[1][1], vertices(ψψ))))) + pψψ_GBP = PartitionedGraph(ψψ, group(v -> v[1][1], vertices(ψψ))) mtsGBP = belief_propagation(pψψ_GBP; contract_kwargs=(; alg="exact"), niters=50) - envsGBP = get_environment(pψψ_GBP, mtsGBP, [(v1, 1), (v1, 2), (v2, 1), (v2, 2)]) + envsGBP = get_environment_tensors(pψψ_GBP, mtsGBP, [(v1, 1), (v1, 2), (v2, 1), (v2, 2)]) ngates = 5 @@ -55,7 +56,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 f2e25f4b..ffcd2430 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -10,12 +10,11 @@ using ITensorNetworks: 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() @@ -37,20 +36,18 @@ ITensors.disable_warn_order() Oψ[v] = apply(op("Sz", s[v]), ψ[v]) exact_sz = contract_inner(Oψ, ψ) / contract_inner(ψ, ψ) - pψψ = PartitionedGraph(ψψ, collect(values(group(v -> v[1], vertices(ψψ))))) - mts = belief_propagation( - pψψ; contract_kwargs=(; alg="exact"), verbose=true, niters=10, target_precision=1e-3 - ) + pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) + mts = belief_propagation(pψψ) numerator_tensors = approx_network_region( - pψψ, mts, [(v, 1)]; verts_tensors=[apply(op("Sz", s[v]), ψ[v])] + pψψ, mts, [PartitionVertex(v)]; verts_tensors=[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))] ) - denominator_tensors = approx_network_region(pψψ, mts, [(v, 1)]) + denominator_tensors = approx_network_region(pψψ, mts, [PartitionVertex(v)]) bp_sz = contract(numerator_tensors)[] / contract(denominator_tensors)[] @test abs.(bp_sz - 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) @@ -64,14 +61,12 @@ ITensors.disable_warn_order() Oψ[v] = apply(op("Sz", s[v]), ψ[v]) exact_sz = contract_inner(Oψ, ψ) / contract_inner(ψ, ψ) - pψψ = PartitionedGraph(ψψ, collect(values(group(v -> v[1], vertices(ψψ))))) - mts = belief_propagation( - pψψ; contract_kwargs=(; alg="exact"), verbose=true, niters=10, target_precision=1e-3 - ) + pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) + mts = belief_propagation(pψψ) numerator_tensors = approx_network_region( - pψψ, mts, [(v, 1)]; verts_tensors=[apply(op("Sz", s[v]), ψ[v])] + pψψ, mts, [PartitionVertex(v)]; verts_tensors=[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))] ) - denominator_tensors = approx_network_region(pψψ, mts, [(v, 1)]) + denominator_tensors = approx_network_region(pψψ, mts, [PartitionVertex(v)]) bp_sz = contract(numerator_tensors)[] / contract(denominator_tensors)[] @test abs.(bp_sz - exact_sz) <= 1e-14 @@ -90,10 +85,8 @@ ITensors.disable_warn_order() ITensors.contract(ψOψ; sequence=contract_seq)[] / ITensors.contract(ψψ; sequence=contract_seq)[] - nsites = 2 - p_vertices = partitioned_vertices(underlying_graph(ψψ); nvertices_per_partition=nsites) - pψψ = PartitionedGraph(ψψ, p_vertices) - mts = belief_propagation(pψψ; contract_kwargs=(; alg="exact"), niters=20) + pψψ = PartitionedGraph(ψψ; nvertices_per_partition=2, backend="Metis") + mts = belief_propagation(pψψ; niters=20) numerator_network = approx_network_region( pψψ, mts, vs; verts_tensors=ITensor[ψOψ[v] for v in vs] ) @@ -105,16 +98,16 @@ ITensors.disable_warn_order() @test abs.(bp_szsz - 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) + 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=[]) - pψψ = PartitionedGraph(ψψ, collect(values(group(v -> v[1], vertices(ψψ))))) - mts = belief_propagation(pψψ; 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])) rdm = ITensors.contract( @@ -134,12 +127,12 @@ ITensors.disable_warn_order() @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) + g_dims = (4, 3) g = named_grid(g_dims) s = siteinds("S=1/2", g) χ = 2 ψ = randomITensorNetwork(s; link_space=χ) - maxdim = 16 + maxdim = 8 v = (2, 2) ψψ = flatten_networks(ψ, dag(ψ); combine_linkinds=false, map_bra_linkinds=prime) @@ -150,12 +143,11 @@ ITensors.disable_warn_order() combiners = linkinds_combiners(ψψ) ψψ = combine_linkinds(ψψ, combiners) ψOψ = combine_linkinds(ψOψ, combiners) - - pψψ = PartitionedGraph(ψψ, collect(values(group(v -> v[1], vertices(ψψ))))) + pψψ = PartitionedGraph(ψψ, group(v -> v[1], vertices(ψψ))) mts = belief_propagation( pψψ; contract_kwargs=(; - alg="density_matrix", output_structure=path_graph_structure, cutoff=1e-16, maxdim + alg="density_matrix", output_structure=path_graph_structure, cutoff=1e-12, maxdim ), ) diff --git a/test/test_examples/test_apply_bp.jl b/test/test_examples/test_apply_bp.jl deleted file mode 100644 index 5e9cabb8..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, pψψ_bp, mts_bp, ψ_vidal, pψψ_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, pψψ_bp, mts_bp) - sz_vidal = expect_bp("Sz", v, ψ_vidal, pψψ_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 124bffc6..b805db50 100644 --- a/test/test_examples/test_examples.jl +++ b/test/test_examples/test_examples.jl @@ -12,10 +12,7 @@ using Test "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"), diff --git a/test/test_examples/test_sqrt_bp.jl b/test/test_examples/test_sqrt_bp.jl deleted file mode 100644 index 127eda83..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) = 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 04067b43..fa003896 100644 --- a/test/test_gauging.jl +++ b/test/test_gauging.jl @@ -34,9 +34,7 @@ using SplitApplyCombine sqrt(contract_inner(ψ_symm, ψ_symm) * contract_inner(ψ, ψ)) ≈ 1.0 ψψ_symm_V2 = ψ_symm ⊗ prime(dag(ψ_symm); sites=[]) - pψψ_symm_V2 = PartitionedGraph( - ψψ_symm_V2, collect(values(group(v -> v[1], vertices(ψψ_symm_V2)))) - ) + pψψ_symm_V2 = PartitionedGraph(ψψ_symm_V2, group(v -> v[1], vertices(ψψ_symm_V2))) ψ_symm_mts_V2 = belief_propagation( pψψ_symm_V2; contract_kwargs=(; alg="exact"), niters=50 ) diff --git a/test/test_indsnetwork.jl b/test/test_indsnetwork.jl index 29a5d858..1a9bc27f 100644 --- a/test/test_indsnetwork.jl +++ b/test/test_indsnetwork.jl @@ -1,7 +1,6 @@ using Dictionaries using ITensors using ITensorNetworks -using ITensorNetworks: dim using Random using Test