diff --git a/src/beliefpropagationdmrg/graphsextensions.jl b/src/beliefpropagationdmrg/graphsextensions.jl deleted file mode 100644 index d605f58a..00000000 --- a/src/beliefpropagationdmrg/graphsextensions.jl +++ /dev/null @@ -1,132 +0,0 @@ -using Graphs: AbstractGraph -using NamedGraphs.GraphsExtensions: a_star, add_edge!, dfs_tree, post_order_dfs_edges -using NamedGraphs.GraphsExtensions: - default_root_vertex, degree, degrees, neighbors, edgetype, rem_edge!, vertextype - -#Given a graph, traverse it from start vertex to end vertex, covering each edge exactly once. -#Complexity is O(length(edges(g))) -function eulerian_path(g::AbstractGraph, start_vertex, end_vertex) - #Conditions on g for the required path to exist - if start_vertex != end_vertex - @assert isodd(degree(g, start_vertex) % 2) - @assert isodd(degree(g, end_vertex) % 2) - @assert all( - x -> iseven(x), degrees(g, setdiff(vertices(g), [start_vertex, end_vertex])) - ) - else - @assert all(x -> iseven(x), degrees(g)) - end - - path = vertextype(g)[] - stack = vertextype(g)[] - current_vertex = end_vertex - g_modified = copy(g) - while !isempty(stack) || !iszero(degree(g_modified, current_vertex)) - if iszero(degree(g_modified, current_vertex)) - push!(path, current_vertex) - last_vertex = pop!(stack) - current_vertex = last_vertex - else - push!(stack, current_vertex) - vn = first(neighbors(g_modified, current_vertex)) - rem_edge!(g_modified, edgetype(g_modified)(current_vertex, vn)) - current_vertex = vn - end - end - - push!(path, current_vertex) - - return edgetype(g_modified)[ - edgetype(g_modified)(path[i], path[i + 1]) for i in 1:(length(path) - 1) - ] -end - -function eulerian_cycle(g::AbstractGraph, start_vertex) - return eulerian_path(g, start_vertex, start_vertex) -end - -function make_all_degrees_even(g::AbstractGraph) - g_modified = copy(g) - vertices_odd_degree = collect(filter(v -> isodd(degree(g, v)), vertices(g_modified))) - while !isempty(vertices_odd_degree) - vertex_pairs = [ - (vertices_odd_degree[i], vertices_odd_degree[j]) for - i in 1:length(vertices_odd_degree) for j in (i + 1):length(vertices_odd_degree) - ] - vertex_pair = first(sort(vertex_pairs; by=vp -> length(a_star(g, vp...)))) - add_edge!(g_modified, edgetype(g_modified)(vertex_pair...)) - vertices_odd_degree = filter( - v -> v != first(vertex_pair) && v != last(vertex_pair), vertices_odd_degree - ) - end - return g_modified -end - -#Given a graph, traverse it in a cycle from start_vertex and try to minimise the number of edges traversed more than once -function _eulerian_cycle(g::AbstractGraph, start_vertex; add_additional_traversal=false) - g_modified = make_all_degrees_even(g::AbstractGraph) - path = eulerian_cycle(g_modified, start_vertex) - - !add_additional_traversal && return path - - modified_path = edgetype(g_modified)[] - - for e in path - if src(e) ∉ neighbors(g, dst(e)) - inner_path = a_star(g, src(e), dst(e)) - append!(modified_path, inner_path) - else - push!(modified_path, e) - end - end - - return modified_path -end - -function _bp_region_plan( - g::AbstractGraph, - start_vertex=default_root_vertex(g); - nsites::Int=1, - add_additional_traversal=false, -) - path = _eulerian_cycle(g, start_vertex; add_additional_traversal) - if nsites == 1 - regions = [[v] for v in vcat(src.(path))] - @assert all(v -> only(v) ∈ vertices(g), regions) - return vcat(regions, reverse(regions)) - else - regions = [e for e in path] - regions = filter(e -> e ∈ edges(g) || reverse(e) ∈ edges(g), regions) - regions = [[src(e), dst(e)] for e in regions] - return vcat(regions, reverse(regions)) - end -end - -function path_to_path(path) - verts = [] - for e in path - if isempty(verts) || src(e) ≠ last(verts) - push!(verts, src(e)) - push!(verts, dst(e)) - end - end - return [[v] for v in verts] -end - -function bp_region_plan( - g::AbstractGraph, - start_vertex=default_root_vertex(g); - nsites::Int=1, - add_additional_traversal=false, -) - path = post_order_dfs_edges(g, start_vertex) - if nsites == 1 - regions = path_to_path(path) - @assert all(v -> only(v) ∈ vertices(g), regions) - return vcat(regions, reverse(regions)) - else - regions = filter(e -> e ∈ edges(g) || reverse(e) ∈ edges(g), path) - @assert all(e -> e ∈ regions || reverse(e) ∈ regions, edges(g)) - return regions - end -end diff --git a/src/beliefpropagationdmrg/utils.jl b/src/beliefpropagationdmrg/utils.jl deleted file mode 100644 index eb10c096..00000000 --- a/src/beliefpropagationdmrg/utils.jl +++ /dev/null @@ -1,202 +0,0 @@ -using ITensors: siteinds, Op, prime, OpSum, apply -using Graphs: AbstractGraph, SimpleGraph, edges, vertices, is_tree, connected_components -using NamedGraphs: NamedGraph, NamedEdge, NamedGraphs, rename_vertices -using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph -using NamedGraphs.GraphsExtensions: - decorate_graph_edges, - forest_cover, - add_edges, - rem_edges, - add_vertices, - rem_vertices, - disjoint_union, - subgraph, - src, - dst -using NamedGraphs.PartitionedGraphs: PartitionVertex, partitionedge, unpartitioned_graph -using ITensorNetworks: - BeliefPropagationCache, - AbstractITensorNetwork, - AbstractFormNetwork, - IndsNetwork, - ITensorNetwork, - insert_linkinds, - ttn, - union_all_inds, - neighbor_vertices, - environment, - messages, - update_factor, - message, - partitioned_tensornetwork, - bra_vertex, - ket_vertex, - operator_vertex, - default_cache_update_kwargs, - dual_index_map, - region_scalar, - renormalize_messages, - scalar_factors_quotient, - norm_sqr_network -using DataGraphs: underlying_graph -using ITensorNetworks.ModelHamiltonians: heisenberg -using ITensors: - ITensor, - noprime, - dag, - noncommonind, - commonind, - replaceind, - dim, - noncommoninds, - delta, - replaceinds, - Trotter, - apply -using ITensors.NDTensors: denseblocks -using SplitApplyCombine: group - -using ITensors: siteinds, Op, prime, OpSum, apply -using Graphs: AbstractGraph, SimpleGraph, edges, vertices, is_tree, connected_components -using NamedGraphs: NamedGraph, NamedEdge, NamedGraphs, rename_vertices -using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph -using NamedGraphs.GraphsExtensions: - decorate_graph_edges, - forest_cover, - add_edges, - rem_edges, - add_vertices, - rem_vertices, - disjoint_union, - subgraph, - src, - dst -using NamedGraphs.PartitionedGraphs: PartitionVertex, partitionedge, unpartitioned_graph -using ITensorNetworks: - BeliefPropagationCache, - AbstractITensorNetwork, - AbstractFormNetwork, - IndsNetwork, - ITensorNetwork, - insert_linkinds, - ttn, - union_all_inds, - neighbor_vertices, - environment, - messages, - update_factor, - message, - partitioned_tensornetwork, - bra_vertex, - ket_vertex, - operator_vertex, - default_cache_update_kwargs, - dual_index_map -using DataGraphs: underlying_graph -using ITensorNetworks.ModelHamiltonians: heisenberg -using ITensors: - ITensor, - noprime, - dag, - noncommonind, - commonind, - replaceind, - dim, - noncommoninds, - delta, - replaceinds -using ITensors.NDTensors: denseblocks -using Dictionaries: set! - -function BP_apply( - o::ITensor, ψ::AbstractITensorNetwork, bpc::BeliefPropagationCache; apply_kwargs... -) - bpc = copy(bpc) - ψ = copy(ψ) - vs = neighbor_vertices(ψ, o) - envs = environment(bpc, PartitionVertex.(vs)) - singular_values! = Ref(ITensor()) - ψ = noprime(apply(o, ψ; envs, singular_values!, normalize=true, apply_kwargs...)) - ψdag = prime(dag(ψ); sites=[]) - if length(vs) == 2 - v1, v2 = vs - pe = partitionedge(bpc, (v1, "bra") => (v2, "bra")) - mts = messages(bpc) - ind1, ind2 = noncommonind(singular_values![], ψ[v1]), - commonind(singular_values![], ψ[v1]) - singular_values![] = denseblocks(replaceind(singular_values![], ind1, ind2')) - set!(mts, pe, ITensor[singular_values![]]) - set!(mts, reverse(pe), ITensor[singular_values![]]) - end - for v in vs - bpc = update_factor(bpc, (v, "ket"), ψ[v]) - bpc = update_factor(bpc, (v, "bra"), ψdag[v]) - end - return ψ, bpc -end - -function get_local_term(bpc::BeliefPropagationCache, v) - qf = tensornetwork(bpc) - return qf[(v, "ket")] * qf[(v, "operator")] * qf[(v, "bra")] -end - -function exact_energy(g::AbstractGraph, bpc::BeliefPropagationCache) - tn = ITensorNetwork(g) - for v in vertices(g) - tn[v] = get_local_term(bpc, v) - end - degree_two_sites = filter(v -> degree(tn, v) == 2, vertices(tn)) - while !isempty(degree_two_sites) - v = first(degree_two_sites) - vn = first(neighbors(g, v)) - tn = contract(tn, NamedEdge(v => vn); merged_vertex=vn) - degree_two_sites = filter(v -> degree(tn, v) == 2, vertices(tn)) - end - return ITensors.contract(ITensor[tn[v] for v in vertices(tn)]; sequence="automatic")[] -end - -function renamer(g) - vertex_rename = Dictionary() - for (i, v) in enumerate(vertices(g)) - set!(vertex_rename, v, (i,)) - end - return rename_vertices(v -> vertex_rename[v], g) -end - -function imaginary_time_evo( - s::IndsNetwork, - ψ::ITensorNetwork, - model::Function, - dbetas::Vector{<:Tuple}; - model_params, - bp_update_kwargs=(; maxiter=10, tol=1e-10), - apply_kwargs=(; cutoff=1e-12, maxdim=10), -) - ψ = copy(ψ) - g = underlying_graph(ψ) - - ℋ = model(g; model_params...) - ψψ = norm_sqr_network(ψ) - bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) - bpc = update(bpc; bp_update_kwargs...) - L = length(vertices(ψ)) - println("Starting Imaginary Time Evolution") - β = 0 - for (i, period) in enumerate(dbetas) - nbetas, dβ = first(period), last(period) - println("Entering evolution period $i , β = $β, dβ = $dβ") - U = exp(-dβ * ℋ; alg=Trotter{1}()) - gates = Vector{ITensor}(U, s) - for i in 1:nbetas - for gate in gates - ψ, bpc = BP_apply(gate, ψ, bpc; apply_kwargs...) - end - β += dβ - bpc = update(bpc; bp_update_kwargs...) - end - e = sum(expect(ψ, ℋ; alg="bp")) - println("Energy is $(e / L)") - end - - return ψ -end