diff --git a/Project.toml b/Project.toml index d8915558..8b6b3c90 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TensorNetworkQuantumSimulator" uuid = "4de3b72a-362e-43dd-83ff-3f381eda9f9c" license = "MIT" -version = "0.3.9" +version = "0.3.10" authors = ["JoeyT1994 ", "MSRudolph ", "and contributors"] description = "A Julia package for quantum simulation with tensor networks of near-arbitrary topology." diff --git a/docs/src/api.md b/docs/src/api.md index 1f3ba7bc..ee8e9ca6 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -83,6 +83,12 @@ update update_iteration! ``` +## Loop Cluster Expansion + +```@docs +connected_edgeinduced_subgraphs_no_leaves +``` + ## Utilities ```@docs diff --git a/examples/loopcorrections.jl b/examples/loopcorrections.jl index 8b599a94..945a2d41 100644 --- a/examples/loopcorrections.jl +++ b/examples/loopcorrections.jl @@ -10,7 +10,7 @@ using Random Random.seed!(1634) function main() - nx, ny = 4, 4 + nx, ny = 4,4 χ = 3 ITensors.disable_warn_order() gs = [ @@ -19,18 +19,32 @@ function main() (named_grid((nx, ny)), "square", 4), ] for (g, g_str, smallest_loop_size) in gs + max_configuration_size = 2*smallest_loop_size -1 + println("\n") + println("-----------------------") + obs = ("Z", first(center(g))) println("Testing for $g_str lattice with $(nv(g)) vertices") - ψ = random_tensornetworkstate(ComplexF32, g, "S=1/2"; bond_dimension = χ) + ψ = random_tensornetworkstate(ComplexF64, g, "S=1/2"; bond_dimension = χ) ψ = normalize(ψ; alg = "bp") norm_bp = norm(ψ; alg = "bp") - norm_loopcorrected = norm(ψ; alg = "loopcorrections", max_configuration_size = 2 * (smallest_loop_size) - 1) + norm_loopcorrected = norm(ψ; alg = "loopcorrections", max_configuration_size) norm_exact = norm(ψ; alg = "exact") println("Bp Value for norm is $norm_bp") println("1st Order Loop Corrected Value for norm is $norm_loopcorrected") println("Exact Value for norm is $norm_exact") + + sz_bp = expect(ψ, obs; alg = "bp") + sz_loopcorrected = expect(ψ, obs; alg = "loopcorrections",max_configuration_size) + sz_exact = expect(ψ, obs; alg = "exact") + + println("\n") + + println("Bp Value for sz is $sz_bp") + println("1st Order Loop Corrected Value for sz is $sz_loopcorrected") + println("Exact Value for sz is $sz_exact") end return end diff --git a/src/Apply/apply_gates.jl b/src/Apply/apply_gates.jl index e9d7beeb..a7f0270f 100644 --- a/src/Apply/apply_gates.jl +++ b/src/Apply/apply_gates.jl @@ -127,4 +127,9 @@ function apply_gate!( return ψ_bpc, err end +function apply_gate(gate::ITensor, ψ_bpc::BeliefPropagationCache; kwargs...) + ψ_bpc = copy(ψ_bpc) + return apply_gate!(gate, ψ_bpc; kwargs...) +end + const apply_circuit = apply_gates diff --git a/src/MessagePassing/LoopExpansion/graph_clusters.jl b/src/MessagePassing/LoopExpansion/graph_clusters.jl new file mode 100644 index 00000000..036e07d7 --- /dev/null +++ b/src/MessagePassing/LoopExpansion/graph_clusters.jl @@ -0,0 +1,177 @@ +# Edge-induced no-leaf subgraph enumeration for the loop / linked-cluster expansion. +# +# This file is PURE GRAPH CODE: it never touches a tensor and operates only on a +# NamedGraph. A "generalized loop" is a connected edge set in which EVERY vertex has +# induced degree >= 2 (no leaves). These are exactly the polymers of the hard-core loop +# gas whose log is the free-energy correction `F − ln Z_BP`, and the configurations whose +# antiprojector weights sum to `Z / Z_BP − 1` (see `weight`/`weights` and +# `loopcorrected_free_energy` in `MessagePassing/loopcorrection.jl`). +# +# NamedGraphs' `edgeinduced_subgraphs_no_leaves` builds candidates as unions of simple +# cycles, so it is connected-only AND silently drops any no-leaf subgraph that needs a +# non-cycle (bridge) edge to be no-leaf — e.g. a "dumbbell" of two triangles joined by a +# single bridge edge (the bridge lies on no cycle, so it is never added; the two bare +# triangles fail the connectivity filter). Those diagrams have all degrees >= 2 and a +# nonzero antiprojector weight, so omitting them under-counts the loop series. The +# enumerator here grows connected edge sets directly, so it is connected *by construction* +# and catches the bridge diagrams. +# +# Enumeration uses ESU (Wernicke, "A faster algorithm for detecting network motifs", +# WABI 2006) on the line graph: connected edge sets of `g` are connected vertex sets of the +# line graph, and ESU visits each connected vertex set EXACTLY ONCE via the "smallest-index +# vertex is the root, only grow with strictly-larger exclusive neighbours" rule. That makes +# the de-duplication structural — no `seen` hash set, and no subgraph is ever generated more +# than once — and the per-node bookkeeping (induced degrees, leaf count) is maintained +# incrementally on integer-indexed arrays instead of rebuilt from a dictionary each step. + +# Apply edge `i` to the running induced-degree state, updating the leaf counter +# (`nleaf` = number of vertices whose induced degree is exactly 1). +@inline function _add_edge!(i, esrc, edst, vdeg, nleaf) + @inbounds for v in (esrc[i], edst[i]) + d = vdeg[v] + vdeg[v] = d + 1 + d == 0 ? (nleaf[] += 1) : d == 1 ? (nleaf[] -= 1) : nothing + end + return nothing +end + +# Inverse of `_add_edge!`: remove edge `i` from the running state. +@inline function _del_edge!(i, esrc, edst, vdeg, nleaf) + @inbounds for v in (esrc[i], edst[i]) + d = vdeg[v] + vdeg[v] = d - 1 + d == 1 ? (nleaf[] -= 1) : d == 2 ? (nleaf[] += 1) : nothing + end + return nothing +end + +# ESU recursion over the line graph. `Vsub` is the current edge-index set (a stack), `ext` +# its extension list (candidate edges, all with index > `root`), and `marked[u]` records +# whether edge `u` is already in `Vsub` or has entered the extension somewhere on the path +# from the root (so it is never re-queued — this is what enforces exact-once visiting). +function _esu_extend!( + out, Vsub, ext, root, eadj, esrc, edst, vdeg, nleaf, marked, max_edges, + ) + # `nleaf == 0` with a non-empty edge set <=> every touched vertex has degree >= 2. + (nleaf[] == 0 && !isempty(Vsub)) && push!(out, copy(Vsub)) + length(Vsub) >= max_edges && return + # Leaf-prune: each added edge clears at most 2 degree-1 vertices, so a set with more + # leaves than 2*(remaining budget) can never become leaf-free within budget. + nleaf[] > 2 * (max_edges - length(Vsub)) && return + added = Int[] + while !isempty(ext) + w = pop!(ext) + newext = copy(ext) # siblings still to be processed (ESU: V_ext \ {w}) + empty!(added) + @inbounds for u in eadj[w] # plus the exclusive neighbours of w, index > root + if u > root && !marked[u] + marked[u] = true + push!(newext, u) + push!(added, u) + end + end + _add_edge!(w, esrc, edst, vdeg, nleaf) + push!(Vsub, w) + _esu_extend!(out, Vsub, newext, root, eadj, esrc, edst, vdeg, nleaf, marked, max_edges) + pop!(Vsub) + _del_edge!(w, esrc, edst, vdeg, nleaf) + @inbounds for u in added # unmark only what THIS w introduced; w stays + marked[u] = false # marked for its siblings (parent unmarks it) + end + end + return nothing +end + +""" + connected_edgeinduced_subgraphs_no_leaves(g, max_edges; anchor = nothing) -> Vector + +All connected edge-induced subgraphs of `g` with at most `max_edges` edges in which every +vertex has induced degree `>= 2` (no leaves) — the "generalized loops" of the loop/linked- +cluster series. Unlike `edgeinduced_subgraphs_no_leaves` (which unions simple cycles and so +misses bridge-joined no-leaf subgraphs such as two triangles linked by a single edge), this +grows connected edge sets directly and is therefore both connected-by-construction and +bridge-complete. + +Enumeration is ESU on the line graph: every connected edge set is visited exactly once (no +deduplication hashing), and induced degrees / the leaf count are tracked incrementally. With +`anchor` set, only subgraphs containing that vertex are produced: the anchor-incident edges +are given the smallest indices and used as the ESU roots, so a subgraph is generated from +its minimum edge iff that edge touches the anchor — i.e. iff the subgraph touches the anchor. +With `anchor = nothing` the whole graph is enumerated. +""" +function connected_edgeinduced_subgraphs_no_leaves( + g::AbstractGraph, max_edges::Integer; anchor = nothing, + ) + E0 = collect(edges(g)) + m = length(E0) + (m == 0 || max_edges < 3) && return typeof(g)[] # smallest no-leaf subgraph is a 3-cycle + + # Integer-index the vertices for O(1) degree bookkeeping. + V = collect(vertices(g)) + vidx = Dict{eltype(V), Int}() + for (k, v) in enumerate(V) + vidx[v] = k + end + + # Order the edge list so that, when anchored, the anchor-incident edges occupy the + # smallest indices `1:nseed`; ESU rooted at `1:nseed` then yields exactly the subgraphs + # whose minimum edge touches the anchor == the subgraphs that touch the anchor. + if anchor === nothing + E = E0 + nseed = m + else + ai = vidx[anchor] + anchored_first = sort(eachindex(E0); by = i -> (vidx[src(E0[i])] != ai && vidx[dst(E0[i])] != ai)) + E = E0[anchored_first] + nseed = count(e -> vidx[src(e)] == ai || vidx[dst(e)] == ai, E) + end + + esrc = Vector{Int}(undef, m) + edst = Vector{Int}(undef, m) + inc = [Int[] for _ in 1:length(V)] # incident edge indices per vertex + for (i, e) in enumerate(E) + a = vidx[src(e)] + b = vidx[dst(e)] + esrc[i] = a + edst[i] = b + push!(inc[a], i) + push!(inc[b], i) + end + # Line-graph adjacency: two edges are adjacent iff they share a vertex. + eadj = [Int[] for _ in 1:m] + for ev in inc, i in ev, j in ev + i != j && push!(eadj[i], j) + end + for i in 1:m + unique!(eadj[i]) + end + + vdeg = zeros(Int, length(V)) + nleaf = Ref(0) + marked = falses(m) + Vsub = Int[] + out = Vector{Int}[] + added = Int[] + for r in 1:nseed + _add_edge!(r, esrc, edst, vdeg, nleaf) + push!(Vsub, r) + marked[r] = true + ext = Int[] + empty!(added) + @inbounds for u in eadj[r] + if u > r && !marked[u] + marked[u] = true + push!(ext, u) + push!(added, u) + end + end + _esu_extend!(out, Vsub, ext, r, eadj, esrc, edst, vdeg, nleaf, marked, max_edges) + pop!(Vsub) + _del_edge!(r, esrc, edst, vdeg, nleaf) + marked[r] = false + @inbounds for u in added + marked[u] = false + end + end + return [edge_subgraph(g, E[S]) for S in out] +end diff --git a/src/MessagePassing/LoopExpansion/loopcorrection.jl b/src/MessagePassing/LoopExpansion/loopcorrection.jl new file mode 100644 index 00000000..3969b86b --- /dev/null +++ b/src/MessagePassing/LoopExpansion/loopcorrection.jl @@ -0,0 +1,157 @@ +using NamedGraphs.GraphsExtensions: boundary_edges + +function loopcorrected_partitionfunction( + bp_cache::BeliefPropagationCache, + max_configuration_size::Integer, + ) + zbp = partitionfunction(bp_cache) + bp_cache = rescale(bp_cache) + egs = + connected_edgeinduced_subgraphs_no_leaves(graph(bp_cache), max_configuration_size) + isempty(egs) && return zbp + ws = weights(bp_cache, egs) + return zbp * (1 + sum(ws)) +end + +# Linked-cluster (free-energy) form of the loop correction. +# +# F = ln Z = ln Z_BP + ln(Z / Z_BP), Z / Z_BP = Σ_{polymer configs} ∏ w_C , +# +# where the polymers are the *connected* no-leaf edge subgraphs (generalized loops) and the +# hard-core exclusion is vertex-sharing. The linked-cluster theorem makes F extensive: its +# expansion is a sum over CONNECTED clusters only — the vertex-disjoint products that the +# partition-function series `loopcorrected_partitionfunction` would have to enumerate beyond +# total size `2·girth − 1` are resummed by the exponential and never appear here. To leading +# cumulant order this is `ln Z_BP + Σ_C w_C` over the connected no-leaf clusters (overlapping- +# cluster corrections, `−½ Σ_{C∼C'} w_C w_{C'} + …`, are higher order). Both this and +# `loopcorrected_partitionfunction` enumerate the same generalized loops via +# `connected_edgeinduced_subgraphs_no_leaves` (so both capture the bridge "dumbbell" +# diagrams); they differ only in how the cluster weights are resummed. +# +# Note `loopcorrected_free_energy` and `log(loopcorrected_partitionfunction)` differ at +# O(w²): they agree exactly only when no loops fit (Σw = 0), where both reduce to `ln Z_BP`. +function loopcorrected_free_energy( + bp_cache::BeliefPropagationCache, + max_configuration_size::Integer, + ) + zbp = partitionfunction(bp_cache) + F = log(complex(zbp)) + bp_cache = rescale(bp_cache) + g = graph(bp_cache) + egs = connected_edgeinduced_subgraphs_no_leaves(g, max_configuration_size) + isempty(egs) && return F + return F + sum(weights(bp_cache, egs)) +end + +# Free-energy (generating-function) estimate of a single-site observable, exposed through +# `expect(...; alg = "loopcorrections")`. +# +# ⟨Ô⟩ = ∂_ε ln⟨ψ|e^{ε Ô}|ψ⟩|_{ε=0} ≈ [F(ε) − F(−ε)] / (2ε), +# F(t) = ln⟨ψ|e^{t Ô}|ψ⟩ = ln‖e^{t Ô/2}|ψ⟩‖² (Hermitian Ô), +# +# Each F is the loop-corrected free energy of a *genuine norm network*: the gate +# e^{±ε Ô/2} is absorbed (un-normalized) into the ket at the observable site and BP is +# re-solved, so the loop expansion prunes ALL leaves — there is no protected operator +# vertex, hence no "anomalous" leaf-containing cluster, which is what tends to make this +# estimator converge smoothly. BP is re-solved on each shifted network so the loop +# corrections also pick up the linear response of the messages. + +# F(α-shifted) = loop-corrected free energy of e^{α Ô}|ψ⟩, built by absorbing the +# (un-normalized) one-site gate into the ket and re-solving BP for the perturbed network. +# `F = ln Z_BP + Σ_C w_C` is the additive linked-cluster free energy (`loopcorrected_free_energy`), +# the genuine extensive log-partition-function whose smooth ε-dependence makes this estimator +# converge well — NOT `log(loopcorrected_partitionfunction) = ln Z_BP + ln(1 + Σ_C w_C)`, which +# only agrees with it to O(w) and resums the same clusters multiplicatively instead. +function gated_lc_free_energy( + ψ_bpc::BeliefPropagationCache, op_string::String, v, α, max_configuration_size::Integer; + cache_update_kwargs, + ) + ψ = network(ψ_bpc) + s = only(siteinds(ψ)[v]) + G = ITensors.exp(α * ITensors.op(op_string, s); ishermitian = true) + # `normalize_tensors = false`: the un-normalized gated tensor is exactly what makes the + # squared norm equal the partition function ⟨ψ|e^{2α Ô}|ψ⟩ we want to differentiate. + ψ_bpc, _ = apply_gate(G, ψ_bpc; v⃗ = [v], apply_kwargs = (; normalize_tensors = false)) + ψ_bpc = update(ψ_bpc; cache_update_kwargs...) + return loopcorrected_free_energy(ψ_bpc, max_configuration_size) +end + +#Transform the indices in the given subgraph of the tensornetwork so that antiprojectors can be inserted without duplicate indices appearing +function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) + bpc = copy(bpc) + vs = collect(vertices(eg)) + es = + unique(collect(Iterators.flatten(boundary_edges(bpc, [v]; dir = :out) for v in vs))) + updated_es = NamedEdge[] + antiprojectors = ITensor[] + for e in es + if reverse(e) ∉ updated_es + mer = message(bpc, reverse(e)) + linds = filter(i -> plev(i) == 0, inds(mer)) + linds_sim = sim.(linds) + mer = replaceinds(mer, linds, linds_sim) + if network(bpc) isa TensorNetworkState + mer = replaceinds(mer, dag.(prime.(linds)), dag.(prime.(linds_sim))) + end + ms = messages(bpc) + set!(ms, reverse(e), mer) + t = network(bpc)[src(e)] + t_inds = filter(i -> i ∈ linds, inds(t)) + if !isempty(t_inds) + t_ind = only(t_inds) + t_ind_pos = findfirst(x -> x == t_ind, linds) + t = replaceind(t, t_ind, linds_sim[t_ind_pos]) + setindex_preserve!(bpc, t, src(e)) + end + push!(updated_es, e) + + if e ∈ edges(eg) || reverse(e) ∈ edges(eg) + row_inds, col_inds = linds, linds_sim + if network(bpc) isa TensorNetworkState + row_inds = vcat(row_inds, dag.(prime.(row_inds))) + col_inds = vcat(col_inds, dag.(prime.(col_inds))) + end + row_combiner, col_combiner = combiner(row_inds), combiner(col_inds) + ap = + adapt_like(message(bpc, e), denseblocks(delta(combinedind(col_combiner), dag(combinedind(row_combiner))))) + ap = ap * row_combiner * dag(col_combiner) + ap = ap - message(bpc, e) * mer + push!(antiprojectors, ap) + end + end + end + return bpc, antiprojectors +end + +#Get the all edges incident to the region specified by the vector of edges passed +function NamedGraphs.GraphsExtensions.boundary_edges( + bpc::BeliefPropagationCache, + es::Vector{<:NamedEdge}, + ) + vs = unique(vcat(src.(es), dst.(es))) + bpes = NamedEdge[] + for v in vs + incoming_es = NamedGraphs.GraphsExtensions.boundary_edges(bpc, [v]; dir = :in) + incoming_es = filter(e -> e ∉ es && reverse(e) ∉ es, incoming_es) + append!(bpes, incoming_es) + end + return bpes +end + +#Compute the contraction of the bp configuration specified by the edge induced subgraph eg +function weight(bpc::BeliefPropagationCache, eg) + vs = collect(vertices(eg)) + es = collect(edges(eg)) + bpc, antiprojectors = sim_edgeinduced_subgraph(bpc, eg) + incoming_ms = + ITensor[message(bpc, e) for e in boundary_edges(bpc, es)] + local_tensors = collect(Iterators.flatten(bp_factors(bpc, v) for v in vs)) + ts = [incoming_ms; local_tensors; antiprojectors] + seq = any(hasqns.(ts)) ? contraction_sequence(ts; alg = "optimal") : contraction_sequence(ts; alg = "einexpr", optimizer = Greedy()) + return scalar(contract(ts; sequence = seq)) +end + +#Vectorized version of weight +function weights(bpc::BeliefPropagationCache, egs) + return [weight(bpc, eg) for eg in egs] +end diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl deleted file mode 100644 index e3b82188..00000000 --- a/src/MessagePassing/loopcorrection.jl +++ /dev/null @@ -1,96 +0,0 @@ -using NamedGraphs.GraphsExtensions: boundary_edges - -function loopcorrected_partitionfunction( - bp_cache::BeliefPropagationCache, - max_configuration_size::Integer, - ) - zbp = partitionfunction(bp_cache) - bp_cache = rescale(bp_cache) - #TODO: Fix edgeinduced_subgraphs_no_leaves for PartitionedGraphView type - #Count the cycles using NamedGraphs - egs = - edgeinduced_subgraphs_no_leaves(graph(bp_cache), max_configuration_size) - isempty(egs) && return zbp - ws = weights(bp_cache, egs) - return zbp * (1 + sum(ws)) -end - -#Transform the indices in the given subgraph of the tensornetwork so that antiprojectors can be inserted without duplicate indices appearing -function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) - bpc = copy(bpc) - vs = collect(vertices(eg)) - es = - unique(collect(Iterators.flatten(boundary_edges(bpc, [v]; dir = :out) for v in vs))) - updated_es = NamedEdge[] - antiprojectors = ITensor[] - for e in es - if reverse(e) ∉ updated_es - mer = message(bpc, reverse(e)) - linds = filter(i -> plev(i) == 0, inds(mer)) - linds_sim = sim.(linds) - mer = replaceinds(mer, linds, linds_sim) - if network(bpc) isa TensorNetworkState - mer = replaceinds(mer, dag.(prime.(linds)), dag.(prime.(linds_sim))) - end - ms = messages(bpc) - set!(ms, reverse(e), mer) - t = network(bpc)[src(e)] - t_inds = filter(i -> i ∈ linds, inds(t)) - if !isempty(t_inds) - t_ind = only(t_inds) - t_ind_pos = findfirst(x -> x == t_ind, linds) - t = replaceind(t, t_ind, linds_sim[t_ind_pos]) - setindex_preserve!(bpc, t, src(e)) - end - push!(updated_es, e) - - if e ∈ edges(eg) || reverse(e) ∈ edges(eg) - row_inds, col_inds = linds, linds_sim - if network(bpc) isa TensorNetworkState - row_inds = vcat(row_inds, dag.(prime.(row_inds))) - col_inds = vcat(col_inds, dag.(prime.(col_inds))) - end - row_combiner, col_combiner = combiner(row_inds), combiner(col_inds) - ap = - adapt_like(message(bpc, e), denseblocks(delta(combinedind(col_combiner), dag(combinedind(row_combiner))))) - ap = ap * row_combiner * dag(col_combiner) - ap = ap - message(bpc, e) * mer - push!(antiprojectors, ap) - end - end - end - return bpc, antiprojectors -end - -#Get the all edges incident to the region specified by the vector of edges passed -function NamedGraphs.GraphsExtensions.boundary_edges( - bpc::BeliefPropagationCache, - es::Vector{<:NamedEdge}, - ) - vs = unique(vcat(src.(es), dst.(es))) - bpes = NamedEdge[] - for v in vs - incoming_es = NamedGraphs.GraphsExtensions.boundary_edges(bpc, [v]; dir = :in) - incoming_es = filter(e -> e ∉ es && reverse(e) ∉ es, incoming_es) - append!(bpes, incoming_es) - end - return bpes -end - -#Compute the contraction of the bp configuration specified by the edge induced subgraph eg -function weight(bpc::BeliefPropagationCache, eg) - vs = collect(vertices(eg)) - es = collect(edges(eg)) - bpc, antiprojectors = sim_edgeinduced_subgraph(bpc, eg) - incoming_ms = - ITensor[message(bpc, e) for e in boundary_edges(bpc, es)] - local_tensors = collect(Iterators.flatten(bp_factors(bpc, v) for v in vs)) - ts = [incoming_ms; local_tensors; antiprojectors] - seq = any(hasqns.(ts)) ? contraction_sequence(ts; alg = "optimal") : contraction_sequence(ts; alg = "einexpr", optimizer = Greedy()) - return scalar(contract(ts; sequence = seq)) -end - -#Vectorized version of weight -function weights(bpc::BeliefPropagationCache, egs) - return [weight(bpc, eg) for eg in egs] -end diff --git a/src/TensorNetworkQuantumSimulator.jl b/src/TensorNetworkQuantumSimulator.jl index 6d0c752b..1fc05706 100644 --- a/src/TensorNetworkQuantumSimulator.jl +++ b/src/TensorNetworkQuantumSimulator.jl @@ -16,7 +16,8 @@ include("Forms/quadraticform.jl") include("MessagePassing/abstractbeliefpropagationcache.jl") include("MessagePassing/beliefpropagationcache.jl") include("MessagePassing/boundarympscache.jl") -include("MessagePassing/loopcorrection.jl") +include("MessagePassing/LoopExpansion/graph_clusters.jl") +include("MessagePassing/LoopExpansion/loopcorrection.jl") include("graph_ops.jl") include("utils.jl") diff --git a/src/expect.jl b/src/expect.jl index c3cb0da0..b39893c2 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -155,6 +155,59 @@ function expect( return expect(alg, ψ_bmps, observable; bmps_messages_up_to_date = true, kwargs...) end +""" + expect(ψ, obs; alg = "loopcorrections", max_configuration_size, ε = 1e-4, cache_update_kwargs...) + +Free-energy / generating-function estimate of a **single-site Hermitian** observable +`obs = (op, vertex[, coeff])`, + + ⟨Ô⟩ = ∂_ε ln⟨ψ|e^{ε Ô}|ψ⟩|_0 ≈ [F(ε) − F(−ε)] / (2ε), + +with `F(t) = ln‖e^{t Ô/2}|ψ⟩‖²` the loop-corrected free energy of that norm network +(`loopcorrected_free_energy`, the additive linked-cluster form `ln Z_BP + Σ_C w_C`; +`max_configuration_size` counts EDGES). BP is re-solved for each shifted +network so the loop corrections include the linear response of the messages. + +`ε` is the central finite-difference step (default `1e-4`). `Ô` must be Hermitian (the +generating-function identity `⟨ψ|e^{εÔ}|ψ⟩ = ‖e^{εÔ/2}|ψ⟩‖²` relies on `Ô† = Ô`). Accepts a +`TensorNetworkState` or a `BeliefPropagationCache`, and a single observable or a vector of +them; multi-site observables are not supported. +""" +function expect( + alg::Algorithm"loopcorrections", + ψ_bpc::BeliefPropagationCache, + obs::Tuple; + max_configuration_size::Integer, ε::Real = 1e-4, + cache_update_kwargs = default_bp_update_kwargs(ψ_bpc), + ) + op_strings, obs_vs, coeff = collectobservable(obs, graph(ψ_bpc)) + iszero(coeff) && return zero(coeff) + length(obs_vs) == 1 || + error("expect(...; alg = \"loopcorrections\") supports single-site observables only; got $(length(obs_vs)) sites.") + v = only(obs_vs) + op_string = only(op_strings) + + Fp = gated_lc_free_energy(ψ_bpc, op_string, v, +ε / 2, max_configuration_size; cache_update_kwargs) + Fm = gated_lc_free_energy(ψ_bpc, op_string, v, -ε / 2, max_configuration_size; cache_update_kwargs) + return coeff * (Fp - Fm) / (2ε) +end + +function expect(alg::Algorithm"loopcorrections", + ψ::TensorNetworkState, obs::Tuple; cache_update_kwargs= default_bp_update_kwargs(ψ), kwargs...) + ψ_bpc = update(BeliefPropagationCache(ψ); cache_update_kwargs...) + return expect(alg, ψ_bpc, obs; kwargs...) +end + +function expect( + alg::Algorithm"loopcorrections", + ψ::Union{TensorNetworkState, BeliefPropagationCache}, + observables::Vector{<:Tuple}; + kwargs..., + ) + return map(obs -> expect(alg, ψ, obs; kwargs...), observables) +end + + #Process an observable into more readable form function collectobservable(obs::Tuple, g::NamedGraph) diff --git a/src/imports.jl b/src/imports.jl index 910796db..044c1ead 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -31,6 +31,7 @@ using NamedGraphs.GraphsExtensions: a_star, add_edge!, edgetype, + edge_subgraph, leaf_vertices, post_order_dfs_edges, decorate_graph_edges, diff --git a/src/utils.jl b/src/utils.jl index 58844c5e..0327c1a0 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -56,7 +56,7 @@ function algorithm_check(tns::Union{AbstractBeliefPropagationCache, TensorNetwor return error("Expected BeliefPropagationCache or TensorNetworkState for 'loop correction' algorithm, got $(typeof(tns))") end - if f ∈ ["normalize", "expect", "sample", "truncate", "rdm"] + if f ∈ ["normalize", "sample", "truncate", "rdm"] return error("Loop correction-based contraction not supported for this functionality yet") end elseif alg == "boundarymps"