From 9ad09bf3468fdc61afe585c1a55e8ac5b6ffc698 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Tue, 28 Oct 2025 16:40:33 -0400 Subject: [PATCH 01/36] local version PEPO --- .gitignore | 5 +- PEPO_2d_ising.jl | 130 +++++++++++++++++++++++++++++++++++++++++++ PEPO_2d_ising_cpu.jl | 119 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 253 insertions(+), 1 deletion(-) create mode 100644 PEPO_2d_ising.jl create mode 100644 PEPO_2d_ising_cpu.jl diff --git a/.gitignore b/.gitignore index 3a2a9d2..175be1d 100644 --- a/.gitignore +++ b/.gitignore @@ -6,4 +6,7 @@ Manifest.toml .vscode/ -.DS_Store \ No newline at end of file +.DS_Store +*.ipynb +*.ipynb_checkpoints +*.~ \ No newline at end of file diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl new file mode 100644 index 0000000..0146b09 --- /dev/null +++ b/PEPO_2d_ising.jl @@ -0,0 +1,130 @@ +using TensorNetworkQuantumSimulator +const TN = TensorNetworkQuantumSimulator + +using ITensorNetworks +const ITN = ITensorNetworks +using ITensors +using ITensors: @OpName_str, @SiteType_str, Algorithm, datatype + +using ITensorNetworks: AbstractBeliefPropagationCache, IndsNetwork, setindex_preserve_graph! +using NamedGraphs +using NamedGraphs: edges +using Graphs +const NG = NamedGraphs +const G = Graphs +using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph +using NamedGraphs.GraphsExtensions: add_edges, add_vertices + +using Random +using TOML + +using Base.Threads +using MKL +using LinearAlgebra +using NPZ + +using CUDA + +using Adapt +using Dictionaries + +BLAS.set_num_threads(min(6, Sys.CPU_THREADS)) +println("Julia is using "*string(nthreads())) +println("BLAS is using "*string(BLAS.get_num_threads())) + +#Gate : rho -> rho .X +function ITensors.op( + ::OpName"X", ::SiteType"Pauli" + ) + mat = zeros(ComplexF64, 4,4) + mat[1,2] = 1 + mat[2,1] = 1 + mat[3,4] = im + mat[4,3] = -im + return mat +end + +function main() + + n = 6 + g = named_grid((n,n)) + s = siteinds(g, "Pauli") + ρ = identitytensornetworkstate(ComplexF64, g, s) + ITensors.disable_warn_order() + + use_gpu = true + + δβ = 0.01 + hx = -3.1 + J = -1 + + # #Do a custom 4-way edge coloring then Trotterise the Hamiltonian into commuting groups + ec1 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 1:2:(n-1)] for i in 1:n]) + ec2 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 2:2:(n-1)] for i in 1:n]) + ec3 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 1:2:(n-1)] for i in 1:n]) + ec4 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 2:2:(n-1)] for i in 1:n]) + ec = [ec1, ec2, ec3, ec4] + + @assert length(reduce(vcat, ec)) == length(edges(g)) + nsteps = 50 + apply_kwargs = (; maxdim = 4, cutoff = 1e-12) + + MPS_message_rank = 10 + + β = 0 + for i in 1:nsteps + #Apply the singsite rotations half way + for v in vertices(g) + gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) + gate = adapt(datatype(ρ), gate) + setindex_preserve_graph!(ρ, normalize(apply(gate, ρ[v])), v) + end + + #Apply the two site rotations, use a boundary MPS cache to apply them (need to run column or row wise depending on the gates) + for (k, colored_edges) in enumerate(ec) + + #Only if you want to use GPU to do boundary MPS + if use_gpu + ρ_gpu =CUDA.cu(ρ) + ρρ = TN.BoundaryMPSCache(ρ_gpu, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) + else + ρρ = TN.BoundaryMPSCache(ρ, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) + end + ρρ = TN.update(ρρ) + TN.update_partitions!(ρρ, collect(TN.partitionvertices(TN.supergraph(ρρ)))) + + for pair in colored_edges + gate = TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s) + gate = adapt(datatype(ρ), gate) + envs = TN.incoming_messages(ρρ, [src(pair), dst(pair)]) + envs = adapt(datatype(ρ)).(envs) + ρv1, ρv2 = ITensorNetworks.full_update_bp(gate, TN.tensornetwork(ρ), [src(pair), dst(pair)]; envs, apply_kwargs...) + ρ[src(pair)], ρ[dst(pair)] = normalize(ρv1), normalize(ρv2) + end + end + + + for v in vertices(g) + gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) + gate = adapt(datatype(ρ), gate) + setindex_preserve_graph!(ρ, normalize(apply(gate, ρ[v])), v) + end + + β += δβ + + if use_gpu + sz_doubled = TN.expect(ρ, ("X", [(2,2)]); alg = "boundarymps", mps_bond_dimension = MPS_message_rank) + else + sz_doubled = TN.expect(CUDA.cu(ρ), ("X", [(2,2)]); alg = "boundarymps", mps_bond_dimension = MPS_message_rank) + end + + println("Inverse Temperature is $β") + println("Bond dimension of PEPO $(ITensorNetworks.maxlinkdim(ρ))") + + println("Expectation value at beta = $(2*β) is $(sz_doubled)") + end + + +end + +main() \ No newline at end of file diff --git a/PEPO_2d_ising_cpu.jl b/PEPO_2d_ising_cpu.jl new file mode 100644 index 0000000..1e9f288 --- /dev/null +++ b/PEPO_2d_ising_cpu.jl @@ -0,0 +1,119 @@ +using TensorNetworkQuantumSimulator +const TN = TensorNetworkQuantumSimulator + +using ITensorNetworks +const ITN = ITensorNetworks +using ITensors +using ITensors: @OpName_str, @SiteType_str, Algorithm, datatype + +using ITensorNetworks: AbstractBeliefPropagationCache, IndsNetwork, setindex_preserve_graph! +using NamedGraphs +using NamedGraphs: edges +using Graphs +const NG = NamedGraphs +const G = Graphs +using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph +using NamedGraphs.GraphsExtensions: add_edges, add_vertices + +using Random +using TOML + +using Base.Threads +using MKL +using LinearAlgebra +using NPZ + +using Adapt +using Dictionaries + +BLAS.set_num_threads(min(6, Sys.CPU_THREADS)) +println("Julia is using "*string(nthreads())) +println("BLAS is using "*string(BLAS.get_num_threads())) + +#Gate : rho -> rho .X +function ITensors.op( + ::OpName"X", ::SiteType"Pauli" + ) + mat = zeros(ComplexF64, 4,4) + mat[1,2] = 1 + mat[2,1] = 1 + mat[3,4] = im + mat[4,3] = -im + return mat +end + +function main() + + n = 6 + g = named_grid((n,n)) + s = siteinds(g, "Pauli") + ρ = identitytensornetworkstate(ComplexF64, g, s) + ITensors.disable_warn_order() + + use_gpu = true + + δβ = 0.01 + hx = -3.1 + J = -1 + + # #Do a custom 4-way edge coloring then Trotterise the Hamiltonian into commuting groups + ec1 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 1:2:(n-1)] for i in 1:n]) + ec2 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 2:2:(n-1)] for i in 1:n]) + ec3 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 1:2:(n-1)] for i in 1:n]) + ec4 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 2:2:(n-1)] for i in 1:n]) + ec = [ec1, ec2, ec3, ec4] + + @assert length(reduce(vcat, ec)) == length(edges(g)) + nsteps = 50 + apply_kwargs = (; maxdim = 4, cutoff = 1e-12) + + MPS_message_rank = 10 + + β = 0 + for i in 1:nsteps + #Apply the singsite rotations half way + for v in vertices(g) + gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) + gate = adapt(datatype(ρ), gate) + setindex_preserve_graph!(ρ, normalize(apply(gate, ρ[v])), v) + end + + #Apply the two site rotations, use a boundary MPS cache to apply them (need to run column or row wise depending on the gates) + for (k, colored_edges) in enumerate(ec) + + ρρ = TN.BoundaryMPSCache(ρ, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) + + ρρ = TN.update(ρρ) + TN.update_partitions!(ρρ, collect(TN.partitionvertices(TN.supergraph(ρρ)))) + + for pair in colored_edges + gate = TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s) + gate = adapt(datatype(ρ), gate) + envs = TN.incoming_messages(ρρ, [src(pair), dst(pair)]) + envs = adapt(datatype(ρ)).(envs) + ρv1, ρv2 = ITensorNetworks.full_update_bp(gate, TN.tensornetwork(ρ), [src(pair), dst(pair)]; envs, apply_kwargs...) + ρ[src(pair)], ρ[dst(pair)] = normalize(ρv1), normalize(ρv2) + end + end + + + for v in vertices(g) + gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) + gate = adapt(datatype(ρ), gate) + setindex_preserve_graph!(ρ, normalize(apply(gate, ρ[v])), v) + end + + β += δβ + + sz_doubled = TN.expect(ρ, ("X", [(2,2)]); alg = "boundarymps", mps_bond_dimension = MPS_message_rank) + + println("Inverse Temperature is $β") + println("Bond dimension of PEPO $(ITensorNetworks.maxlinkdim(ρ))") + + println("Expectation value at beta = $(2*β) is $(sz_doubled)") + end + + +end + +main() \ No newline at end of file From 27229facde6c760bdeceed85941b5f8d8ccf4bcc Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Tue, 28 Oct 2025 17:26:51 -0400 Subject: [PATCH 02/36] minor changes test --- PEPO_2d_ising.jl | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index 0146b09..ffa44c6 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -44,18 +44,14 @@ function ITensors.op( return mat end -function main() +function main(n, hx, nsteps; use_gpu::Bool = true, χ::Int=4, MPS_message_rank::Int = 10) - n = 6 g = named_grid((n,n)) s = siteinds(g, "Pauli") ρ = identitytensornetworkstate(ComplexF64, g, s) ITensors.disable_warn_order() - use_gpu = true - δβ = 0.01 - hx = -3.1 J = -1 # #Do a custom 4-way edge coloring then Trotterise the Hamiltonian into commuting groups @@ -66,10 +62,7 @@ function main() ec = [ec1, ec2, ec3, ec4] @assert length(reduce(vcat, ec)) == length(edges(g)) - nsteps = 50 apply_kwargs = (; maxdim = 4, cutoff = 1e-12) - - MPS_message_rank = 10 β = 0 for i in 1:nsteps From dd56cb2f0df1f2efd9e4e15c81442bf2c241bba1 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Fri, 31 Oct 2025 15:13:37 -0400 Subject: [PATCH 03/36] split main() into many functions, added bp --- PEPO_2d_ising.jl | 171 +++++++++++++++++++++++++++++-------------- PEPO_2d_ising_cpu.jl | 4 +- Project.toml | 25 +++++-- 3 files changed, 133 insertions(+), 67 deletions(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index ffa44c6..c9749c5 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -27,11 +27,14 @@ using CUDA using Adapt using Dictionaries +using JLD2 BLAS.set_num_threads(min(6, Sys.CPU_THREADS)) println("Julia is using "*string(nthreads())) println("BLAS is using "*string(BLAS.get_num_threads())) +DATA_DIR = "/mnt/ceph/users/gsommers/data/" + #Gate : rho -> rho .X function ITensors.op( ::OpName"X", ::SiteType"Pauli" @@ -44,80 +47,136 @@ function ITensors.op( return mat end -function main(n, hx, nsteps; use_gpu::Bool = true, χ::Int=4, MPS_message_rank::Int = 10) +function prep_edges(n::Int, g::AbstractNamedGraph) + # #Do a custom 4-way edge coloring then Trotterize the Hamiltonian into commuting groups + ec1 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 1:2:(n-1)] for i in 1:n]) + ec2 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 2:2:(n-1)] for i in 1:n]) + ec3 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 1:2:(n-1)] for i in 1:n]) + ec4 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 2:2:(n-1)] for i in 1:n]) + ec = [ec1, ec2, ec3, ec4] + + @assert length(reduce(vcat, ec)) == length(edges(g)) + ec +end + +# apply layer of single qubit gates +function apply_single_qubit_layer!(ρ::TensorNetworkState, gates::Dict) + for v=keys(gates) + setindex_preserve_graph!(ρ, normalize(apply(gates[v], ρ[v])), v) + end +end + +#Apply the two site rotations, use a boundary MPS cache to apply them (need to run column or row wise depending on the gates) +function apply_two_qubit_layer!(ρ::TensorNetworkState, ec::Array, gates::Dict; MPS_message_rank::Int, use_gpu::Bool=true, apply_kwargs...) + for (k, colored_edges) in enumerate(ec) + + #Only if you want to use GPU to do boundary MPS + if use_gpu + ρ_gpu =CUDA.cu(ρ) + ρρ = TN.BoundaryMPSCache(ρ_gpu, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) + else + ρρ = TN.BoundaryMPSCache(ρ, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) + end + @time ρρ = TN.update(ρρ) + TN.update_partitions!(ρρ, collect(TN.partitionvertices(TN.supergraph(ρρ)))) + + for pair in colored_edges + apply_two_qubit_gate!(ρ,ρρ, gates[pair], pair; apply_kwargs...) + end + end +end + +function apply_two_qubit_gate!(ρ::TensorNetworkState,ρρ::TN.BoundaryMPSCache, gate::ITensor, pair::NamedEdge; apply_kwargs...) + envs = TN.incoming_messages(ρρ, [src(pair), dst(pair)]) + envs = adapt(datatype(ρ)).(envs) + ρv1, ρv2 = ITensorNetworks.full_update_bp(gate, TN.tensornetwork(ρ), [src(pair), dst(pair)]; envs, apply_kwargs...) + ρ[src(pair)], ρ[dst(pair)] = normalize(ρv1), normalize(ρv2) +end + +function intermediate_save(sqrtρ, β; δβ::Float64, χ::Int, n::Int, MPS_message_rank, save_tag = "", hx = -3.04438) + dat = Dict("L"=>n, "δβ"=>δβ, "β"=>β, "χ"=>χ, "sqrtρ"=>sqrtρ, "mps_rank"=>MPS_message_rank, "hx"=>hx) + save(DATA_DIR * "$(save_tag)L$(n)_χ$(χ)_D$(MPS_message_rank)_step$(round(δβ,digits=3))_$(round(β,digits=3)).jld2", dat) +end + +function intermediate_save_bp(ρ, errs, β; δβ::Float64, χ::Int, n::Int, save_tag = "", hx = -3.04438) + dat = Dict("L"=>n, "δβ"=>δβ, "β"=>β, "χ"=>χ, "ρ"=>ρ, "errs"=>errs, "hx"=>hx) + save(DATA_DIR * "$(save_tag)L$(n)_χ$(χ)_step$(round(δβ,digits=3))_$(round(β,digits=3)).jld2", dat) +end + +function expect_bmps(n::Int, nsteps::Int; hx=-3.04438, δβ = 0.01, χ::Int=4, MPS_message_rank::Int = 10, save_tag = "", load_tag = "") +end + +function evolve_bmps(n::Int, nsteps::Int; hx=-3.04438, δβ = 0.01, use_gpu::Bool = true, χ::Int=4, MPS_message_rank::Int = 10, save_tag = "") g = named_grid((n,n)) - s = siteinds(g, "Pauli") + s = siteinds("Pauli", g) ρ = identitytensornetworkstate(ComplexF64, g, s) + evolve_bmps(ρ, n, nsteps; β=0, hx=hx, δβ=δβ, use_gpu = use_gpu, χ=χ, MPS_message_rank = MPS_message_rank, save_tag = save_tag) + +end + +function evolve_bmps(ρ::TensorNetworkState, n::Int, nsteps::Int; β = 0, hx=-3.04438, δβ = 0.01, use_gpu::Bool=true, χ::Int=4, MPS_message_rank::Int = 10, save_tag = "") + g = ρ.tensornetwork.data_graph.underlying_graph + s = siteinds(ρ) ITensors.disable_warn_order() - δβ = 0.01 J = -1 - # #Do a custom 4-way edge coloring then Trotterise the Hamiltonian into commuting groups - ec1 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 1:2:(n-1)] for i in 1:n]) - ec2 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 2:2:(n-1)] for i in 1:n]) - ec3 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 1:2:(n-1)] for i in 1:n]) - ec4 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 2:2:(n-1)] for i in 1:n]) - ec = [ec1, ec2, ec3, ec4] + ec = prep_edges(n, g) + apply_kwargs = (; maxdim = χ, cutoff = 1e-12) - @assert length(reduce(vcat, ec)) == length(edges(g)) - apply_kwargs = (; maxdim = 4, cutoff = 1e-12) + two_qubit_gates = Dict(pair=>adapt(datatype(ρ), TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s)) for pair=vcat(ec...)) + + single_qubit_gates = Dict(v=>adapt(datatype(ρ), TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s)) for v=vertices(g)) - β = 0 for i in 1:nsteps #Apply the singsite rotations half way - for v in vertices(g) - gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) - gate = adapt(datatype(ρ), gate) - setindex_preserve_graph!(ρ, normalize(apply(gate, ρ[v])), v) - end + apply_single_qubit_layer!(ρ, single_qubit_gates) - #Apply the two site rotations, use a boundary MPS cache to apply them (need to run column or row wise depending on the gates) - for (k, colored_edges) in enumerate(ec) - - #Only if you want to use GPU to do boundary MPS - if use_gpu - ρ_gpu =CUDA.cu(ρ) - ρρ = TN.BoundaryMPSCache(ρ_gpu, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) - else - ρρ = TN.BoundaryMPSCache(ρ, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) - end - ρρ = TN.update(ρρ) - TN.update_partitions!(ρρ, collect(TN.partitionvertices(TN.supergraph(ρρ)))) - - for pair in colored_edges - gate = TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s) - gate = adapt(datatype(ρ), gate) - envs = TN.incoming_messages(ρρ, [src(pair), dst(pair)]) - envs = adapt(datatype(ρ)).(envs) - ρv1, ρv2 = ITensorNetworks.full_update_bp(gate, TN.tensornetwork(ρ), [src(pair), dst(pair)]; envs, apply_kwargs...) - ρ[src(pair)], ρ[dst(pair)] = normalize(ρv1), normalize(ρv2) - end - end + @time apply_two_qubit_layer!(ρ, ec, two_qubit_gates; MPS_message_rank = MPS_message_rank, use_gpu = use_gpu, apply_kwargs...) - - for v in vertices(g) - gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) - gate = adapt(datatype(ρ), gate) - setindex_preserve_graph!(ρ, normalize(apply(gate, ρ[v])), v) - end + apply_single_qubit_layer!(ρ, single_qubit_gates) β += δβ - if use_gpu - sz_doubled = TN.expect(ρ, ("X", [(2,2)]); alg = "boundarymps", mps_bond_dimension = MPS_message_rank) - else - sz_doubled = TN.expect(CUDA.cu(ρ), ("X", [(2,2)]); alg = "boundarymps", mps_bond_dimension = MPS_message_rank) - end - - println("Inverse Temperature is $β") - println("Bond dimension of PEPO $(ITensorNetworks.maxlinkdim(ρ))") + println("Inverse Temperature is $(2*β)"); flush(stdout) - println("Expectation value at beta = $(2*β) is $(sz_doubled)") + intermediate_save(ρ,β; χ=χ,n=n,MPS_message_rank = MPS_message_rank, δβ=δβ, hx=hx, save_tag = save_tag) end +end - +function evolve_bp(n::Int, nsteps::Int; hx=-3.04438, δβ = 0.01, use_gpu::Bool=true, χ::Int=4, save_tag = "") + g = named_grid((n,n)) + s = siteinds("Pauli", g) + ρ = identitytensornetworkstate(ComplexF64, g, s) + evolve_bp(ρ, n, nsteps; β=0, hx=hx, δβ=δβ, use_gpu = use_gpu, χ=χ, save_tag = save_tag) end -main() \ No newline at end of file +function evolve_bp(ρ::TensorNetworkState, n::Int, nsteps::Int; β = 0, hx=-3.04438, δβ = 0.01, use_gpu::Bool=true, χ::Int=4, save_tag = "") + g = ρ.tensornetwork.data_graph.underlying_graph + s = siteinds(ρ) + ITensors.disable_warn_order() + + J = -1 + + ec = prep_edges(n, g) + apply_kwargs = (; maxdim = χ, cutoff = 1e-12) + + ρρ = BeliefPropagationCache(ρ) + + two_qubit_gates = [adapt(datatype(ρ), TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s)) for pair=vcat(ec...)] + + single_qubit_gates = [adapt(datatype(ρ), TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s)) for v=vertices(g)] + + layer = vcat(single_qubit_gates, two_qubit_gates, single_qubit_gates) + for i in 1:nsteps + + @time ρρ, errs = apply_gates(layer, ρρ; apply_kwargs, verbose = false) + + β += δβ + + println("Inverse Temperature is $(2*β)"); flush(stdout) + println("X_center: $(expect(ρρ, ("X", [(n÷2,n÷2)])))") + intermediate_save_bp(ρρ,errs,β; χ=χ,n=n, δβ=δβ, hx=hx, save_tag = save_tag) + end +end diff --git a/PEPO_2d_ising_cpu.jl b/PEPO_2d_ising_cpu.jl index 1e9f288..5ea6e6a 100644 --- a/PEPO_2d_ising_cpu.jl +++ b/PEPO_2d_ising_cpu.jl @@ -46,12 +46,10 @@ function main() n = 6 g = named_grid((n,n)) - s = siteinds(g, "Pauli") + s = siteinds("Pauli", g) ρ = identitytensornetworkstate(ComplexF64, g, s) ITensors.disable_warn_order() - use_gpu = true - δβ = 0.01 hx = -3.1 J = -1 diff --git a/Project.toml b/Project.toml index aad56cd..1584c8f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,14 +1,16 @@ name = "TensorNetworkQuantumSimulator" uuid = "4de3b72a-362e-43dd-83ff-3f381eda9f9c" -authors = ["JoeyT1994 ", "MSRudolph ","and contributors"] -description = "A Julia package for simulating quantum circuits and dynamics with tensor networks of near-arbritrary topology." license = "MIT" -version = "0.1.02" +authors = ["JoeyT1994 ", "MSRudolph ", "and contributors"] +description = "A Julia package for simulating quantum circuits and dynamics with tensor networks of near-arbritrary topology." +version = "0.1.2" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a" +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" +DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a" +Debugger = "31a5f54b-26ea-5ae9-a837-f05ce5417438" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" EinExprs = "b1794770-133b-4de1-afb4-526377e9f4c5" GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73" @@ -16,7 +18,10 @@ Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensorNetworks = "2919e153-833c-4bdc-8836-1ea460a35fc7" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" +NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605" NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19" PauliPropagation = "293282d5-3c99-4fb6-92d0-fd3280a19750" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" @@ -30,8 +35,10 @@ TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" [compat] Adapt = "4.3.0" -DataGraphs = "0.2.9" +CUDA = "5.9.2" Combinatorics = "1.0.3" +DataGraphs = "0.2.9" +Debugger = "0.7.15" Dictionaries = "0.4" EinExprs = "0.6.4" GraphRecipes = "0.5.13" @@ -39,7 +46,10 @@ Graphs = "1.8.0" ITensorMPS = "0.3.17" ITensorNetworks = "0.14.2" ITensors = "0.9" +JLD2 = "0.6.2" LinearAlgebra = "1.11.0" +MKL = "0.9.0" +NPZ = "0.4.3" NamedGraphs = "0.7" PauliPropagation = "0.3.0" Revise = "3.8.0" @@ -52,10 +62,9 @@ TensorOperations = "5.2" TypeParameterAccessors = "0.3.10" julia = "1.10" - [extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "Random"] \ No newline at end of file +test = ["Test", "Random"] From 09db0c6faee96db78828cb215bcd3f37222dc9c4 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Fri, 31 Oct 2025 20:59:01 -0400 Subject: [PATCH 04/36] expect_bmps --- PEPO_2d_ising.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index c9749c5..953b140 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -103,7 +103,14 @@ function intermediate_save_bp(ρ, errs, β; δβ::Float64, χ::Int, n::Int, save save(DATA_DIR * "$(save_tag)L$(n)_χ$(χ)_step$(round(δβ,digits=3))_$(round(β,digits=3)).jld2", dat) end -function expect_bmps(n::Int, nsteps::Int; hx=-3.04438, δβ = 0.01, χ::Int=4, MPS_message_rank::Int = 10, save_tag = "", load_tag = "") +function expect_bmps(dat::Dict; obs = "X", MPS_message_rank::Int = 10, save_tag = "") + all_verts = collect(vertices(dat["sqrtρ"][1].tensornetwork.data_graph.underlying_graph)) + expect_vals = zeros(length(all_verts), length(dat["sqrtρ"])) + for i=1:length(dat["sqrtρ"]) + @time expect_vals[:,i] = real.(TN.expect(dat["sqrtρ"][i], [(obs, [v]) for v=all_verts]; alg = "boundarymps", mps_bond_dimension = MPS_message_rank)) + save(DATA_DIR * "$(save_tag)L$(n)_χ$(dat["χ"])_D$(MPS_message_rank)_step$(dat["δβ"])_$(dat["β"][i]).jld2", Dict(obs=>expect_vals[:,i], "verts"=>all_verts, "hx"=>dat["hx"], "β"=>dat["β"][i], χ=>[dat["χ"],maxlinkdim(dat["sqrtρ"][i])], mps_rank=>MPS_message_rank, "δβ"=>dat["δβ"], "L"=>dat["L"])) + end + all_verts, expect_vals end function evolve_bmps(n::Int, nsteps::Int; hx=-3.04438, δβ = 0.01, use_gpu::Bool = true, χ::Int=4, MPS_message_rank::Int = 10, save_tag = "") From 56db07312e5aa239fffc5dab1872430fc9269d28 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Fri, 31 Oct 2025 21:42:45 -0400 Subject: [PATCH 05/36] fixed typo expect_bmps --- PEPO_2d_ising.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index 953b140..13c4012 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -108,7 +108,7 @@ function expect_bmps(dat::Dict; obs = "X", MPS_message_rank::Int = 10, save_tag expect_vals = zeros(length(all_verts), length(dat["sqrtρ"])) for i=1:length(dat["sqrtρ"]) @time expect_vals[:,i] = real.(TN.expect(dat["sqrtρ"][i], [(obs, [v]) for v=all_verts]; alg = "boundarymps", mps_bond_dimension = MPS_message_rank)) - save(DATA_DIR * "$(save_tag)L$(n)_χ$(dat["χ"])_D$(MPS_message_rank)_step$(dat["δβ"])_$(dat["β"][i]).jld2", Dict(obs=>expect_vals[:,i], "verts"=>all_verts, "hx"=>dat["hx"], "β"=>dat["β"][i], χ=>[dat["χ"],maxlinkdim(dat["sqrtρ"][i])], mps_rank=>MPS_message_rank, "δβ"=>dat["δβ"], "L"=>dat["L"])) + save(DATA_DIR * "$(save_tag)L$(dat["L"])_χ$(dat["χ"])_D$(MPS_message_rank)_step$(dat["δβ"])_$(dat["β"][i]).jld2", Dict(obs=>expect_vals[:,i], "verts"=>all_verts, "hx"=>dat["hx"], "β"=>dat["β"][i], χ=>[dat["χ"],maxlinkdim(dat["sqrtρ"][i])], mps_rank=>MPS_message_rank, "δβ"=>dat["δβ"], "L"=>dat["L"])) end all_verts, expect_vals end From 14bb36ea8194037a179d4e4c0a8ece706931dc2c Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Fri, 31 Oct 2025 21:47:19 -0400 Subject: [PATCH 06/36] adding print statements --- PEPO_2d_ising.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index 13c4012..c804871 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -71,6 +71,7 @@ function apply_two_qubit_layer!(ρ::TensorNetworkState, ec::Array, gates::Dict; for (k, colored_edges) in enumerate(ec) #Only if you want to use GPU to do boundary MPS + println("Starting boundary MPS cache") if use_gpu ρ_gpu =CUDA.cu(ρ) ρρ = TN.BoundaryMPSCache(ρ_gpu, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) @@ -80,8 +81,9 @@ function apply_two_qubit_layer!(ρ::TensorNetworkState, ec::Array, gates::Dict; @time ρρ = TN.update(ρρ) TN.update_partitions!(ρρ, collect(TN.partitionvertices(TN.supergraph(ρρ)))) + println("Starting two-qubit gates") for pair in colored_edges - apply_two_qubit_gate!(ρ,ρρ, gates[pair], pair; apply_kwargs...) + @time apply_two_qubit_gate!(ρ,ρρ, gates[pair], pair; apply_kwargs...) end end end From 8ae63054f055156805c67a5dc896e2b9b2027af5 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sun, 2 Nov 2025 09:25:53 -0500 Subject: [PATCH 07/36] added time checks, use_gpu for expect_bmps --- PEPO_2d_ising.jl | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index c804871..ef3832e 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -73,18 +73,20 @@ function apply_two_qubit_layer!(ρ::TensorNetworkState, ec::Array, gates::Dict; #Only if you want to use GPU to do boundary MPS println("Starting boundary MPS cache") if use_gpu - ρ_gpu =CUDA.cu(ρ) - ρρ = TN.BoundaryMPSCache(ρ_gpu, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) + @time ρ_gpu =CUDA.cu(ρ) + @time ρρ = TN.BoundaryMPSCache(ρ_gpu, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) else - ρρ = TN.BoundaryMPSCache(ρ, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) + @time ρρ = TN.BoundaryMPSCache(ρ, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) end @time ρρ = TN.update(ρρ) - TN.update_partitions!(ρρ, collect(TN.partitionvertices(TN.supergraph(ρρ)))) + @time TN.update_partitions!(ρρ, collect(TN.partitionvertices(TN.supergraph(ρρ)))) println("Starting two-qubit gates") - for pair in colored_edges - @time apply_two_qubit_gate!(ρ,ρρ, gates[pair], pair; apply_kwargs...) - end + @time begin + for pair in colored_edges + apply_two_qubit_gate!(ρ,ρρ, gates[pair], pair; apply_kwargs...) + end + end end end @@ -105,11 +107,16 @@ function intermediate_save_bp(ρ, errs, β; δβ::Float64, χ::Int, n::Int, save save(DATA_DIR * "$(save_tag)L$(n)_χ$(χ)_step$(round(δβ,digits=3))_$(round(β,digits=3)).jld2", dat) end -function expect_bmps(dat::Dict; obs = "X", MPS_message_rank::Int = 10, save_tag = "") +function expect_bmps(dat::Dict; obs = "X", MPS_message_rank::Int = 10, save_tag = "", use_gpu::Bool = true) all_verts = collect(vertices(dat["sqrtρ"][1].tensornetwork.data_graph.underlying_graph)) expect_vals = zeros(length(all_verts), length(dat["sqrtρ"])) for i=1:length(dat["sqrtρ"]) - @time expect_vals[:,i] = real.(TN.expect(dat["sqrtρ"][i], [(obs, [v]) for v=all_verts]; alg = "boundarymps", mps_bond_dimension = MPS_message_rank)) + if use_gpu + sqrtρ = CUDA.cu(dat["sqrtρ"][i]) + else + sqrtρ = dat["sqrtρ"][i] + end + @time expect_vals[:,i] = real.(TN.expect(dat["sqrtρ"][i], [(obs, [v]) for v=all_verts]; alg = "boundarymps", mps_bond_dimension = MPS_message_rank)) save(DATA_DIR * "$(save_tag)L$(dat["L"])_χ$(dat["χ"])_D$(MPS_message_rank)_step$(dat["δβ"])_$(dat["β"][i]).jld2", Dict(obs=>expect_vals[:,i], "verts"=>all_verts, "hx"=>dat["hx"], "β"=>dat["β"][i], χ=>[dat["χ"],maxlinkdim(dat["sqrtρ"][i])], mps_rank=>MPS_message_rank, "δβ"=>dat["δβ"], "L"=>dat["L"])) end all_verts, expect_vals From 1f44ba2ef3c3c24380ad4a630a14f99f4b72113b Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Mon, 3 Nov 2025 21:49:51 -0500 Subject: [PATCH 08/36] fixed small bugs --- PEPO_2d_ising.jl | 32 ++++++++++++++++---------------- Project.toml | 4 ++++ 2 files changed, 20 insertions(+), 16 deletions(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index ef3832e..1f90976 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -104,20 +104,22 @@ end function intermediate_save_bp(ρ, errs, β; δβ::Float64, χ::Int, n::Int, save_tag = "", hx = -3.04438) dat = Dict("L"=>n, "δβ"=>δβ, "β"=>β, "χ"=>χ, "ρ"=>ρ, "errs"=>errs, "hx"=>hx) + dat["X"] = expect(ρ, [("X", [v]) for v=vertices(network(ρ).tensornetwork.data_graph.underlying_graph)]) save(DATA_DIR * "$(save_tag)L$(n)_χ$(χ)_step$(round(δβ,digits=3))_$(round(β,digits=3)).jld2", dat) end -function expect_bmps(dat::Dict; obs = "X", MPS_message_rank::Int = 10, save_tag = "", use_gpu::Bool = true) +function expect_bmps(dat::Dict; obs = "X", MPS_message_rank::Int = 10, save_tag = "", use_gpu::Bool = true, start_i::Int = 1) all_verts = collect(vertices(dat["sqrtρ"][1].tensornetwork.data_graph.underlying_graph)) - expect_vals = zeros(length(all_verts), length(dat["sqrtρ"])) - for i=1:length(dat["sqrtρ"]) + expect_vals = zeros(length(all_verts), length(dat["sqrtρ"])-start_i+1) + for i=start_i:length(dat["sqrtρ"]) if use_gpu sqrtρ = CUDA.cu(dat["sqrtρ"][i]) else sqrtρ = dat["sqrtρ"][i] end - @time expect_vals[:,i] = real.(TN.expect(dat["sqrtρ"][i], [(obs, [v]) for v=all_verts]; alg = "boundarymps", mps_bond_dimension = MPS_message_rank)) - save(DATA_DIR * "$(save_tag)L$(dat["L"])_χ$(dat["χ"])_D$(MPS_message_rank)_step$(dat["δβ"])_$(dat["β"][i]).jld2", Dict(obs=>expect_vals[:,i], "verts"=>all_verts, "hx"=>dat["hx"], "β"=>dat["β"][i], χ=>[dat["χ"],maxlinkdim(dat["sqrtρ"][i])], mps_rank=>MPS_message_rank, "δβ"=>dat["δβ"], "L"=>dat["L"])) + @time expect_vals[:,i-start_i+1] = real.(TN.expect(sqrtρ, [(obs, [v]) for v=all_verts]; alg = "boundarymps", mps_bond_dimension = MPS_message_rank)) + save(DATA_DIR * "$(save_tag)L$(dat["L"])_χ$(dat["χ"])_D$(MPS_message_rank)_step$(dat["δβ"])_$(dat["β"][i]).jld2", Dict(obs=>expect_vals[:,i-start_i+1], "verts"=>all_verts, "hx"=>dat["hx"], "β"=>dat["β"][i], "χ"=>[dat["χ"],maxlinkdim(dat["sqrtρ"][i])], "mps_rank"=>MPS_message_rank, "δβ"=>dat["δβ"], "L"=>dat["L"])) + flush(stdout) end all_verts, expect_vals end @@ -155,7 +157,7 @@ function evolve_bmps(ρ::TensorNetworkState, n::Int, nsteps::Int; β = 0, hx=-3. β += δβ - println("Inverse Temperature is $(2*β)"); flush(stdout) + println("Inverse Temperature is $(β)"); flush(stdout) intermediate_save(ρ,β; χ=χ,n=n,MPS_message_rank = MPS_message_rank, δβ=δβ, hx=hx, save_tag = save_tag) end @@ -165,24 +167,23 @@ function evolve_bp(n::Int, nsteps::Int; hx=-3.04438, δβ = 0.01, use_gpu::Bool= g = named_grid((n,n)) s = siteinds("Pauli", g) ρ = identitytensornetworkstate(ComplexF64, g, s) - evolve_bp(ρ, n, nsteps; β=0, hx=hx, δβ=δβ, use_gpu = use_gpu, χ=χ, save_tag = save_tag) + ρρ = BeliefPropagationCache(ρ) + evolve_bp(ρρ, n, nsteps; β=0, hx=hx, δβ=δβ, use_gpu = use_gpu, χ=χ, save_tag = save_tag) end -function evolve_bp(ρ::TensorNetworkState, n::Int, nsteps::Int; β = 0, hx=-3.04438, δβ = 0.01, use_gpu::Bool=true, χ::Int=4, save_tag = "") - g = ρ.tensornetwork.data_graph.underlying_graph - s = siteinds(ρ) +function evolve_bp(ρρ::BeliefPropagationCache, n::Int, nsteps::Int; β = 0, hx=-3.04438, δβ = 0.01, use_gpu::Bool=true, χ::Int=4, save_tag = "") + g = network(ρρ).tensornetwork.data_graph.underlying_graph + s = siteinds(network(ρρ)) ITensors.disable_warn_order() J = -1 ec = prep_edges(n, g) apply_kwargs = (; maxdim = χ, cutoff = 1e-12) - - ρρ = BeliefPropagationCache(ρ) - two_qubit_gates = [adapt(datatype(ρ), TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s)) for pair=vcat(ec...)] + two_qubit_gates = [adapt(datatype(network(ρρ)), TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s)) for pair=vcat(ec...)] - single_qubit_gates = [adapt(datatype(ρ), TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s)) for v=vertices(g)] + single_qubit_gates = [adapt(datatype(network(ρρ)), TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s)) for v=vertices(g)] layer = vcat(single_qubit_gates, two_qubit_gates, single_qubit_gates) for i in 1:nsteps @@ -191,8 +192,7 @@ function evolve_bp(ρ::TensorNetworkState, n::Int, nsteps::Int; β = 0, hx=-3.04 β += δβ - println("Inverse Temperature is $(2*β)"); flush(stdout) - println("X_center: $(expect(ρρ, ("X", [(n÷2,n÷2)])))") + println("Inverse Temperature is $(β)"); flush(stdout) intermediate_save_bp(ρρ,errs,β; χ=χ,n=n, δβ=δβ, hx=hx, save_tag = save_tag) end end diff --git a/Project.toml b/Project.toml index 1584c8f..936ca01 100644 --- a/Project.toml +++ b/Project.toml @@ -24,6 +24,7 @@ MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605" NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19" PauliPropagation = "293282d5-3c99-4fb6-92d0-fd3280a19750" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SimpleGraphAlgorithms = "41400c72-0c58-5c16-8579-4ecbce768449" SimpleGraphConverter = "205b04f2-f585-4877-a239-566270b3f673" @@ -32,6 +33,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" +UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] Adapt = "4.3.0" @@ -52,6 +54,7 @@ MKL = "0.9.0" NPZ = "0.4.3" NamedGraphs = "0.7" PauliPropagation = "0.3.0" +Plots = "1.41.1" Revise = "3.8.0" SimpleGraphAlgorithms = "0.6.0" SimpleGraphConverter = "0.1.0" @@ -60,6 +63,7 @@ Statistics = "1.11.1" StatsBase = "0.34.4" TensorOperations = "5.2" TypeParameterAccessors = "0.3.10" +UnicodePlots = "3.8.1" julia = "1.10" [extras] From 8576d8351bee3275c9f314010ccd8befacaf0acd Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Thu, 6 Nov 2025 15:49:18 -0500 Subject: [PATCH 09/36] Project.toml --- .gitignore | 4 +++- Project.toml | 16 +++++++++++++++- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index 3a2a9d2..47ba7fc 100644 --- a/.gitignore +++ b/.gitignore @@ -6,4 +6,6 @@ Manifest.toml .vscode/ -.DS_Store \ No newline at end of file +.DS_Store +.* +*.ipynb \ No newline at end of file diff --git a/Project.toml b/Project.toml index 14f9dcd..800537c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,23 +1,30 @@ name = "TensorNetworkQuantumSimulator" uuid = "4de3b72a-362e-43dd-83ff-3f381eda9f9c" license = "MIT" +version = "0.2.2" authors = ["JoeyT1994 ", "MSRudolph ", "and contributors"] description = "A Julia package for quantum simulation with tensor networks of near-arbritrary topology." -version = "0.2.2" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" EinExprs = "b1794770-133b-4de1-afb4-526377e9f4c5" +Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8" GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19" PauliPropagation = "293282d5-3c99-4fb6-92d0-fd3280a19750" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SimpleGraphAlgorithms = "41400c72-0c58-5c16-8579-4ecbce768449" SimpleGraphConverter = "205b04f2-f585-4877-a239-566270b3f673" @@ -32,14 +39,21 @@ Adapt = "4.3.0" Combinatorics = "1.0.3" Dictionaries = "0.4" EinExprs = "0.6.4" +Format = "1.3.7" GraphRecipes = "0.5.13" Graphs = "1.8.0" +IJulia = "1.32.1" ITensorMPS = "0.3.17" ITensors = "0.9" +JLD2 = "0.6.2" KrylovKit = "0.10.2" +LaTeXStrings = "1.4.0" LinearAlgebra = "1.11.0" +LsqFit = "0.15.1" NamedGraphs = "0.7" PauliPropagation = "0.3.0" +Plots = "1.41.1" +PythonPlot = "1.0.6" Revise = "3.8.0" SimpleGraphAlgorithms = "0.6.0" SimpleGraphConverter = "0.1.0" From 29634a87b607ee8012a312fc98a64ad6a3b7c6c9 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Thu, 6 Nov 2025 16:00:11 -0500 Subject: [PATCH 10/36] PEPO_2d_ising no more ITensorNetworks --- PEPO_2d_ising.jl | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index 1f90976..bfad484 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -1,14 +1,9 @@ using TensorNetworkQuantumSimulator const TN = TensorNetworkQuantumSimulator -using ITensorNetworks -const ITN = ITensorNetworks -using ITensors -using ITensors: @OpName_str, @SiteType_str, Algorithm, datatype - -using ITensorNetworks: AbstractBeliefPropagationCache, IndsNetwork, setindex_preserve_graph! -using NamedGraphs -using NamedGraphs: edges +using ITensors: @OpName_str, @SiteType_str, Algorithm, datatype, ITensors + +using NamedGraphs: NamedGraphs, edges, NamedEdge using Graphs const NG = NamedGraphs const G = Graphs @@ -62,7 +57,7 @@ end # apply layer of single qubit gates function apply_single_qubit_layer!(ρ::TensorNetworkState, gates::Dict) for v=keys(gates) - setindex_preserve_graph!(ρ, normalize(apply(gates[v], ρ[v])), v) + setindex_preserve_graph!(ρ, normalize(ITensors.apply(gates[v], ρ[v])), v) end end @@ -93,8 +88,9 @@ end function apply_two_qubit_gate!(ρ::TensorNetworkState,ρρ::TN.BoundaryMPSCache, gate::ITensor, pair::NamedEdge; apply_kwargs...) envs = TN.incoming_messages(ρρ, [src(pair), dst(pair)]) envs = adapt(datatype(ρ)).(envs) - ρv1, ρv2 = ITensorNetworks.full_update_bp(gate, TN.tensornetwork(ρ), [src(pair), dst(pair)]; envs, apply_kwargs...) - ρ[src(pair)], ρ[dst(pair)] = normalize(ρv1), normalize(ρv2) + ρv1, ρv2 = TN.full_update(gate, ρ, [src(pair), dst(pair)]; envs, print_fidelity_loss = true, apply_kwargs...) + TN.setindex_preserve!(ρ, normalize(ρv1), src(pair)) + TN.setindex_preserve!(ρ, normalize(ρv2), dst(pair)) end function intermediate_save(sqrtρ, β; δβ::Float64, χ::Int, n::Int, MPS_message_rank, save_tag = "", hx = -3.04438) From 1608c31acc28c5fd511a5d01ad4709c2faf07213 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Thu, 6 Nov 2025 20:16:34 -0500 Subject: [PATCH 11/36] fixed dependencies in PEPO_2d_ising.jl --- PEPO_2d_ising.jl | 2 +- Project.toml | 8 +++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index bfad484..c4b80c2 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -42,7 +42,7 @@ function ITensors.op( return mat end -function prep_edges(n::Int, g::AbstractNamedGraph) +function prep_edges(n::Int, g::NamedGraphs.AbstractNamedGraph) # #Do a custom 4-way edge coloring then Trotterize the Hamiltonian into commuting groups ec1 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 1:2:(n-1)] for i in 1:n]) ec2 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 2:2:(n-1)] for i in 1:n]) diff --git a/Project.toml b/Project.toml index 25159ca..d06b647 100644 --- a/Project.toml +++ b/Project.toml @@ -1,12 +1,13 @@ name = "TensorNetworkQuantumSimulator" uuid = "4de3b72a-362e-43dd-83ff-3f381eda9f9c" license = "MIT" -version = "0.2.2" authors = ["JoeyT1994 ", "MSRudolph ", "and contributors"] description = "A Julia package for quantum simulation with tensor networks of near-arbritrary topology." +version = "0.2.2" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a" Debugger = "31a5f54b-26ea-5ae9-a837-f05ce5417438" @@ -23,6 +24,8 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" +MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" +NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605" NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19" PauliPropagation = "293282d5-3c99-4fb6-92d0-fd3280a19750" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" @@ -39,6 +42,7 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] Adapt = "4.3.0" +CUDA = "5.9.2" Combinatorics = "1.0.3" DataGraphs = "0.2.9" Debugger = "0.7.15" @@ -55,6 +59,8 @@ KrylovKit = "0.10.2" LaTeXStrings = "1.4.0" LinearAlgebra = "1.11.0" LsqFit = "0.15.1" +MKL = "0.9.0" +NPZ = "0.4.3" NamedGraphs = "0.7" PauliPropagation = "0.3.0" Plots = "1.41.1" From 5a43d66bad091e3759ce9ccaeb9f9944bc01c878 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Thu, 6 Nov 2025 20:27:40 -0500 Subject: [PATCH 12/36] fixed another bug --- PEPO_2d_ising.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index c4b80c2..a0885ef 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -85,7 +85,7 @@ function apply_two_qubit_layer!(ρ::TensorNetworkState, ec::Array, gates::Dict; end end -function apply_two_qubit_gate!(ρ::TensorNetworkState,ρρ::TN.BoundaryMPSCache, gate::ITensor, pair::NamedEdge; apply_kwargs...) +function apply_two_qubit_gate!(ρ::TensorNetworkState,ρρ::TN.BoundaryMPSCache, gate::ITensors.ITensor, pair::NamedEdge; apply_kwargs...) envs = TN.incoming_messages(ρρ, [src(pair), dst(pair)]) envs = adapt(datatype(ρ)).(envs) ρv1, ρv2 = TN.full_update(gate, ρ, [src(pair), dst(pair)]; envs, print_fidelity_loss = true, apply_kwargs...) From 7dfad0065383c65110e86efc791dba15675b045a Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Thu, 6 Nov 2025 20:37:09 -0500 Subject: [PATCH 13/36] cluster folder --- cluster/expect-corrected.jl | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 cluster/expect-corrected.jl diff --git a/cluster/expect-corrected.jl b/cluster/expect-corrected.jl new file mode 100644 index 0000000..938358e --- /dev/null +++ b/cluster/expect-corrected.jl @@ -0,0 +1,2 @@ +function test() +end \ No newline at end of file From bcc34b73b1cff1dc897232d0ef0c8a696b7dde36 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Thu, 6 Nov 2025 20:57:55 -0500 Subject: [PATCH 14/36] started implementing expect-corrected --- Project.toml | 2 + cluster/expect-corrected.jl | 266 +++++++++++++++++++++++++++++++++++- 2 files changed, 266 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 800537c..101cdfd 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ EinExprs = "b1794770-133b-4de1-afb4-526377e9f4c5" Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8" GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97" IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" @@ -42,6 +43,7 @@ EinExprs = "0.6.4" Format = "1.3.7" GraphRecipes = "0.5.13" Graphs = "1.8.0" +HyperDualNumbers = "4.0.10" IJulia = "1.32.1" ITensorMPS = "0.3.17" ITensors = "0.9" diff --git a/cluster/expect-corrected.jl b/cluster/expect-corrected.jl index 938358e..137eb90 100644 --- a/cluster/expect-corrected.jl +++ b/cluster/expect-corrected.jl @@ -1,2 +1,264 @@ -function test() -end \ No newline at end of file +using TensorNetworkQuantumSimulator +using NamedGraphs +using NamedGraphs: AbstractGraph, NamedGraph, AbstractNamedGraph +using NamedGraphs.GraphsExtensions: add_edge +using Graphs + +const G = Graphs +const NG = NamedGraphs +const TN = TensorNetworkQuantumSimulator +using TensorNetworkQuantumSimulator: BoundaryMPSCache +using HyperDualNumbers +using Adapt: adapt +using Dictionaries +using ITensors: inds, onehot, dag, commonind, op + +# Cluster/loop/CC weights +function cluster_weights(bpc::BeliefPropagationCache, clusters::Vector, egs::Vector{<:AbstractNamedGraph}, interaction_graph; rescale::Bool = true) + logZbp = logscalar(bpc) + if rescale + bpc = TN.rescale(bpc) + end + + isempty(egs) && return [0], [[logZbp]], [[1]] + + circuit_lengths = sort(unique([c.weight for c=clusters])) + + # calculate weight of each generalized loop first + wts = TN.weights(bpc, egs) + + logZs = Array{Array}(undef, length(circuit_lengths) + 1) + logZs[1] = [logZbp] + + coeffs = Array{Array}(undef, length(circuit_lengths) + 1) + coeffs[1] = [1] + + # now calculate contribution to logZ from each cluster + for (cl_i, cl)=enumerate(circuit_lengths) + clusters_cl = filter(c->c.weight==cl, clusters) + logZs[cl_i + 1] = [prod([prod(fill(wts[l],c.multiplicities[l])) for l=c.loop_ids]) for c=clusters_cl] + # weird bug in HyperDualNumbers where this doesn't work... + # sum_ws = sum([TN.ursell_function(c, interaction_graph) * prod([wts[l]^c.multiplicities[l] for l=c.loop_ids]) for c=clusters_cl]) + coeffs[cl_i + 1] = [TN.ursell_function(c, interaction_graph) for c=clusters_cl] + end + + return vcat([0],circuit_lengths), logZs, coeffs +end + +function cc_weights(bpc::BeliefPropagationCache, regions::Vector, counting_nums::Dict; rescale::Bool=true) + logZbp = logscalar(bpc) + if rescale + bpc = TN.rescale(bpc) + end + + use_g = findall(gg->counting_nums[gg] != 0, regions) + egs = [induced_subgraph(bpc.partitioned_tensornetwork.partitioned_graph, gg)[1] for gg=regions[use_g]] + + isempty(egs) && return [0], [[logZbp]], [[1]] + + # calculate weight of each generalized loop first + wts = TN.full_weights(bpc, egs) + + return logZbp, log.(wts), [counting_nums[gg] for gg=regions[use_g]] +end + + +# Loop series expansion + +function terms_to_scalar(numerator_numerator_terms, numerator_denominator_terms, denominator_numerator_terms, denominator_denominator_terms) + return exp(sum(log.(numerator_numerator_terms)) - sum(log.(numerator_denominator_terms)) - sum(log.(denominator_numerator_terms)) + sum(log.(denominator_denominator_terms))) +end + +#Project spins on sites v1 and v2 to v1_val (1 = up, 2 = down) and v2_val +function project!(ψIψ::BeliefPropagationCache, v1, v2, v1_val::Int64 = 1, v2_val::Int64=1) + s1 = only(inds(only(ITN.factors(ψIψ, [(v1, "operator")])); plev = 0)) + s2 = only(inds(only(ITN.factors(ψIψ, [(v2, "operator")])); plev = 0)) + ITensorNetworks.@preserve_graph ψIψ[(v1, "operator")] = onehot(s1 => v1_val) * dag(onehot(s1' => v1_val)) + ITensorNetworks.@preserve_graph ψIψ[(v2, "operator")] = onehot(s2 => v2_val) * dag(onehot(s2' => v2_val)) +end + +#Log scalar contraction of bpc +function logscalar(bpc::BeliefPropagationCache) + nums, denoms = TN.scalar_factors_quotient(bpc) + return sum(log.(nums)) - sum(log.(denoms)) +end + +function cumulative_weights(bpc::BeliefPropagationCache, egs::Vector{<:AbstractNamedGraph}) + isempty(egs) && [1] + circuit_lengths = sort(unique(length.(edges.(egs)))) + outs = [] + for cl in circuit_lengths + egs_cl = filter(eg -> length(edges(eg)) == cl, egs) + sum_ws = sum(TN.weights(bpc, egs_cl)) + push!(outs, sum_ws) + end + + # first one is the BP contribution + outs = vcat([1.0], outs) + return cumsum(outs) +end + +function compute_ps!(ψIψ::BeliefPropagationCache, v1, v2, v1_val, v2_val, egs::Vector{<:AbstractNamedGraph}; track::Bool = true, cache_update_kwargs...) + project!(ψIψ, v1, v2, v1_val, v2_val) + + if track + ψIψ, bp_diffs = updatecache(ψIψ; track = true, cache_update_kwargs...) + else + ψIψ = updatecache(ψIψ; track = false, cache_update_kwargs...) + bp_diffs = [] + end + + p_bp = exp(logscalar(ψIψ)) + + ψIψ = ITensorNetworks.rescale(ψIψ) + cfes = cumulative_weights(ψIψ, egs) + + return [p_bp*cfe for cfe in cfes], bp_diffs +end + +function scalar_cumul_loop(bp_cache::BeliefPropagationCache, egs) + zbp = exp(logscalar(bp_cache)) + bp_cache = ITN.rescale(bp_cache) + + cfes = cumulative_weights(bp_cache, egs) + + return [zbp * cfe for cfe=cfes] +end + +function zz_bp_loopcorrect_connected(ψIψ::BeliefPropagationCache, verts, egs::Vector{<:AbstractNamedGraph}; cache_update_kwargs...) + z_expects = [z_bp_loopcorrect(ψIψ, v, egs; cache_update_kwargs...)[1] for v=verts] + zz_expect = zz_bp_loopcorrect(ψIψ, verts...,egs; cache_update_kwargs...) + + zz = (zz_expect[1,1] .+ zz_expect[2,2] .- (zz_expect[1,2] .+ zz_expect[2,1])) ./sum(zz_expect) + vcat([zz], [(z[1] .- z[2]) ./ (z[1] .+ z[2]) for z=z_expects]) +end + +# compute z expectation value, but return cumulative weights rather than just total loop-corrected value +function z_bp_loopcorrect(ψIψ::BeliefPropagationCache, v, egs::Vector{<:AbstractNamedGraph}; cache_update_kwargs...) + ψUpψ = TN.insert_observable(ψIψ, ("Proj0", [v])) + ψUpψ, diffs_up = updatecache(ψUpψ; track = true, cache_update_kwargs...) + p_up = scalar_cumul_loop(ψUpψ, egs) + ψDnψ = TN.insert_observable(ψIψ, ("Proj1", [v])) + ψDnψ, diffs_dn = updatecache(ψDnψ; track = true, cache_update_kwargs...) + p_dn = scalar_cumul_loop(ψDnψ, egs) + + denominator = scalar_cumul_loop(ψIψ, egs) + return [p_up, p_dn, denominator], [diffs_up, diffs_dn] +end + +#Compute zz with loop correction and all bells and whistles +function zz_bp_loopcorrect(ψIψ::BeliefPropagationCache, v1, v2, egs::Vector{<:AbstractNamedGraph}; cache_update_kwargs...) + probs = [compute_ps!(copy(ψIψ), v1, v2, i, j, egs; track = false, cache_update_kwargs...)[1] for i=1:2, j=1:2] +end + +function cluster_twopoint_correlator(ψIψ::BeliefPropagationCache, obs, max_weight::Int) + ng = ITN.partitioned_graph(ψIψ) + op_strings, verts, _ = TN.collectobservable(obs) + @assert length(verts)<=2 + clusters, egs, interaction_graph = enumerate_connected_clusters_twopoint(ng, verts..., max_weight) + + cluster_twopoint_correlator(ψIψ, obs, clusters, egs, interaction_graph) +end + +function cluster_twopoint_correlator(ψIψ::BeliefPropagationCache, obs, clusters, egs, interaction_graph) + # rescale BEFORE inserting operator + ψIψ = TN.rescale(ψIψ) + op_strings, verts, _ = TN.collectobservable(obs) + @assert length(verts)<=2 + + ψIψ_vs = [ψIψ[(v, "operator")] for v in verts] + sinds = + [commonind(ψIψ[(v, "ket")], ψIψ_vs[i]) for (i, v) in enumerate(verts)] + + + coeffs = [Hyper(0,1,0,0),Hyper(0,0,1,0)] + operators = [ψIψ[(v, "operator")].tensor[1,1] * adapt(typeof(ψIψ[(v, "operator")]))(Hyper(1,0,0,0) * op("I", sinds[i]) + coeffs[i] * op(op_strings[i], sinds[i])) for (i, v) in enumerate(verts)] + ψOψ = ITN.update_factors(ψIψ, Dictionary([(v, "operator") for v in verts], operators)) + + cluster_weights(ψOψ, clusters, egs, interaction_graph; rescale = false) +end + +function cluster_onepoint(ψIψ::BeliefPropagationCache, obs, clusters, egs, interaction_graph) + ψOψ, logZbp = prep_insertion(ψIψ, obs) + lens, logZs, coeffs = cluster_weights(ψOψ, clusters, egs, interaction_graph; rescale = false) + logZs[1] .+= logZbp + lens, logZs, coeffs +end + +# copied over from ITensorNetworks, + +# Also return the amount rescaled by +function rescale_partitions_norms( + bpc::ITN.AbstractBeliefPropagationCache, + partitions::Vector; + verts::Vector=vertices(bpc, partitions), +) + bpc = copy(bpc) + tn = ITN.tensornetwork(bpc) + + # not sure why this is done in two steps + norms = map(v -> inv(ITN.norm(tn[v])), verts) + ITN.scale!(bpc, Dictionary(verts, norms)) + + vertices_weights = Dictionary() + for pv in partitions + pv_vs = filter(v -> v ∈ verts, vertices(bpc, pv)) + isempty(pv_vs) && continue + + vn = ITN.region_scalar(bpc, pv) + s = one(vn)#sign(vn) #isreal(vn) ? sign(vn) : one(vn) + vn = s * vn^(-inv(oftype(vn, length(pv_vs)))) + set!(vertices_weights, first(pv_vs), s*vn) + for v in pv_vs[2:length(pv_vs)] + set!(vertices_weights, v, vn) + end + end + + ITN.scale!(bpc, vertices_weights) + + return bpc, Dictionary(verts,norms), vertices_weights +end + +# cluster cumulant expansion + +function cc_twopoint_correlator(ψIψ::BeliefPropagationCache, obs, regions::Vector, counting_nums::Dict) + # rescale BEFORE inserting operator + ψIψ = TN.rescale(ψIψ) + op_strings, verts, _ = TN.collectobservable(obs) + @assert length(verts)==2 + + ψIψ_vs = [ψIψ[(v, "operator")] for v in verts] + sinds = + [commonind(ψIψ[(v, "ket")], ψIψ_vs[i]) for (i, v) in enumerate(verts)] + + coeffs = [Hyper(0,1,0,0),Hyper(0,0,1,0)] + operators = [ψIψ[(v, "operator")].tensor[1,1] * adapt(typeof(ψIψ[(v, "operator")]))(Hyper(1,0,0,0) * op("I", sinds[i]) + coeffs[i] * op(op_strings[i], sinds[i])) for (i, v) in enumerate(verts)] + ψOψ = ITN.update_factors(ψIψ, Dictionary([(v, "operator") for v in verts], operators)) + + cc_weights(ψOψ, regions, counting_nums; rescale=false) +end + +function cc_onepoint(ψIψ::BeliefPropagationCache, obs, regions::Vector, counting_nums::Dict) + ψOψ, logZbp = prep_insertion(ψIψ, obs) + lz, logZs, coeffs = cc_weights(ψOψ, regions, counting_nums; rescale=false) + logZbp + lz, logZs, coeffs +end + +function prep_insertion(ψIψ::BeliefPropagationCache, obs) + # rescale BEFORE inserting operator + ψIψ = ITN.rescale_messages(ψIψ) + op_strings, verts, _ = TN.collectobservable(obs) + @assert length(verts)<=2 + + ψIψ_vs = [ψIψ[(v, "operator")] for v in verts] + sinds = + [commonind(ψIψ[(v, "ket")], ψIψ_vs[i]) for (i, v) in enumerate(verts)] + + coeffs = [Hyper(0,1,0,0),Hyper(0,0,1,0)] + operators = [adapt(typeof(ψIψ[(v, "operator")]))(Hyper(1,0,0,0) * op("I", sinds[i]) + coeffs[i] * op(op_strings[i], sinds[i])) for (i, v) in enumerate(verts)] + ψOψ = ITN.update_factors(ψIψ, Dictionary([(v, "operator") for v in verts], operators)) + + ψOψ, norms, vertices_weights = rescale_partitions_norms(ψOψ, collect(ITN.partitions(ψIψ))) + return ψOψ, sum([log(1/(norms[(v,tag)] * vertices_weights[(v,tag)])) for v=verts, tag = ["bra", "ket", "operator"]]) +end + From bf30deeed4f8d26fe110cc8d5f3d4a29db0434fe Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Fri, 7 Nov 2025 13:43:30 -0500 Subject: [PATCH 15/36] fixed typos in compatibility --- PEPO_2d_ising.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index a0885ef..e1c0abd 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -57,7 +57,7 @@ end # apply layer of single qubit gates function apply_single_qubit_layer!(ρ::TensorNetworkState, gates::Dict) for v=keys(gates) - setindex_preserve_graph!(ρ, normalize(ITensors.apply(gates[v], ρ[v])), v) + TN.setindex_preserve!(ρ, normalize(ITensors.apply(gates[v], ρ[v])), v) end end @@ -100,12 +100,12 @@ end function intermediate_save_bp(ρ, errs, β; δβ::Float64, χ::Int, n::Int, save_tag = "", hx = -3.04438) dat = Dict("L"=>n, "δβ"=>δβ, "β"=>β, "χ"=>χ, "ρ"=>ρ, "errs"=>errs, "hx"=>hx) - dat["X"] = expect(ρ, [("X", [v]) for v=vertices(network(ρ).tensornetwork.data_graph.underlying_graph)]) + dat["X"] = expect(ρ, [("X", [v]) for v=vertices(network(ρ).tensornetwork.graph)]) save(DATA_DIR * "$(save_tag)L$(n)_χ$(χ)_step$(round(δβ,digits=3))_$(round(β,digits=3)).jld2", dat) end function expect_bmps(dat::Dict; obs = "X", MPS_message_rank::Int = 10, save_tag = "", use_gpu::Bool = true, start_i::Int = 1) - all_verts = collect(vertices(dat["sqrtρ"][1].tensornetwork.data_graph.underlying_graph)) + all_verts = collect(vertices(dat["sqrtρ"][1].tensornetwork.graph)) expect_vals = zeros(length(all_verts), length(dat["sqrtρ"])-start_i+1) for i=start_i:length(dat["sqrtρ"]) if use_gpu @@ -130,7 +130,7 @@ function evolve_bmps(n::Int, nsteps::Int; hx=-3.04438, δβ = 0.01, use_gpu::Boo end function evolve_bmps(ρ::TensorNetworkState, n::Int, nsteps::Int; β = 0, hx=-3.04438, δβ = 0.01, use_gpu::Bool=true, χ::Int=4, MPS_message_rank::Int = 10, save_tag = "") - g = ρ.tensornetwork.data_graph.underlying_graph + g = ρ.tensornetwork.graph s = siteinds(ρ) ITensors.disable_warn_order() @@ -168,7 +168,7 @@ function evolve_bp(n::Int, nsteps::Int; hx=-3.04438, δβ = 0.01, use_gpu::Bool= end function evolve_bp(ρρ::BeliefPropagationCache, n::Int, nsteps::Int; β = 0, hx=-3.04438, δβ = 0.01, use_gpu::Bool=true, χ::Int=4, save_tag = "") - g = network(ρρ).tensornetwork.data_graph.underlying_graph + g = network(ρρ).tensornetwork.graph s = siteinds(network(ρρ)) ITensors.disable_warn_order() From bec5ebf89460083fa56d33757a5c278d911b50bc Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Fri, 7 Nov 2025 13:44:36 -0500 Subject: [PATCH 16/36] joey's version --- Finitetemp_2d_ising_fullupdate.jl | 127 ++++++++++++++++++++++++++++++ 1 file changed, 127 insertions(+) create mode 100644 Finitetemp_2d_ising_fullupdate.jl diff --git a/Finitetemp_2d_ising_fullupdate.jl b/Finitetemp_2d_ising_fullupdate.jl new file mode 100644 index 0000000..3d93d74 --- /dev/null +++ b/Finitetemp_2d_ising_fullupdate.jl @@ -0,0 +1,127 @@ +using TensorNetworkQuantumSimulator +const TN = TensorNetworkQuantumSimulator + +using ITensors: @OpName_str, @SiteType_str, Algorithm, datatype, ITensors + +using NamedGraphs: NamedGraphs, edges, NamedEdge +using Graphs +const NG = NamedGraphs +const G = Graphs +using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph +using NamedGraphs.GraphsExtensions: add_edges, add_vertices + +using Random +using TOML + +using Base.Threads +using MKL +using LinearAlgebra +using NPZ + +using CUDA + +using Adapt +using Dictionaries + +BLAS.set_num_threads(min(6, Sys.CPU_THREADS)) +println("Julia is using "*string(nthreads())) +println("BLAS is using "*string(BLAS.get_num_threads())) + +#Gate : rho -> rho .X. With this defined, expect(sqrt_rho, X; alg) = Tr(sqrt_rho . X sqrt_rho) / Tr(sqrt_rho sqrt_rho) = Tr(rho . X) / Tr(rho) +function ITensors.op( + ::OpName"X", ::SiteType"Pauli" + ) + mat = zeros(ComplexF64, 4,4) + mat[1,2] = 1 + mat[2,1] = 1 + mat[3,4] = im + mat[4,3] = -im + return mat +end + +function main() + + n = 6 + g = named_grid((n,n)) + #Pauli inds run over identity, X, Y, Z + s = siteinds("Pauli", g) + ρ = identitytensornetworkstate(ComplexF64, g, s) + ITensors.disable_warn_order() + + use_gpu = true + + δβ = 0.01 + hx = -3.1 + J = -1 + + # #Do a custom 4-way edge coloring then Trotterise the Hamiltonian into commuting groups + ec1 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 1:2:(n-1)] for i in 1:n]) + ec2 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 2:2:(n-1)] for i in 1:n]) + ec3 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 1:2:(n-1)] for i in 1:n]) + ec4 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 2:2:(n-1)] for i in 1:n]) + ec = [ec1, ec2, ec3, ec4] + + @assert length(reduce(vcat, ec)) == length(edges(g)) + nsteps = 50 + apply_kwargs = (; maxdim = 4, cutoff = 1e-12) + + MPS_message_rank = 10 + + β = 0 + for i in 1:nsteps + #Apply the singsite rotations half way + for v in vertices(g) + gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) + gate = adapt(datatype(ρ), gate) + TN.setindex_preserve!(ρ, normalize(ITensors.apply(gate, ρ[v])), v) + end + + #Apply the two site rotations, use a boundary MPS cache to apply them (need to run column or row wise depending on the gates) + for (k, colored_edges) in enumerate(ec) + + #Only if you want to use GPU to do boundary MPS + if use_gpu + ρ_gpu =CUDA.cu(ρ) + ρρ = TN.BoundaryMPSCache(ρ_gpu, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) + else + ρρ = TN.BoundaryMPSCache(ρ, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) + end + ρρ = TN.update(ρρ) + TN.update_partitions!(ρρ, collect(TN.partitionvertices(TN.supergraph(ρρ)))) + + for pair in colored_edges + gate = TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s) + gate = adapt(datatype(ρ), gate) + envs = TN.incoming_messages(ρρ, [src(pair), dst(pair)]) + envs = adapt(datatype(ρ)).(envs) + ρv1, ρv2 = TN.full_update(gate, ρ, [src(pair), dst(pair)]; envs, print_fidelity_loss = true, apply_kwargs...) + TN.setindex_preserve!(ρ, normalize(ρv1), src(pair)) + TN.setindex_preserve!(ρ, normalize(ρv2), dst(pair)) + end + end + + + for v in vertices(g) + gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) + gate = adapt(datatype(ρ), gate) + TN.setindex_preserve!(ρ, normalize(ITensors.apply(gate, ρ[v])), v) + end + + β += δβ + + if use_gpu + sz_doubled = TN.expect(ρ, ("X", [(2,2)]); alg = "boundarymps", mps_bond_dimension = MPS_message_rank) + else + sz_doubled = TN.expect(CUDA.cu(ρ), ("X", [(2,2)]); alg = "boundarymps", mps_bond_dimension = MPS_message_rank) + end + + println("Inverse Temperature is $β") + println("Bond dimension of PEPO $(TN.maxvirtualdim(ρ))") + + println("Expectation value at beta = $(2*β) is $(sz_doubled)") + end + + +end + +main() \ No newline at end of file From aa7fdcb3ebbaccaba3db2896cab83e27aca2098f Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sat, 8 Nov 2025 10:05:59 -0500 Subject: [PATCH 17/36] gpu for bp --- Finitetemp_2d_ising_fullupdate.jl | 2 +- PEPO_2d_ising.jl | 5 ++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/Finitetemp_2d_ising_fullupdate.jl b/Finitetemp_2d_ising_fullupdate.jl index 3d93d74..7577c5e 100644 --- a/Finitetemp_2d_ising_fullupdate.jl +++ b/Finitetemp_2d_ising_fullupdate.jl @@ -48,7 +48,7 @@ function main() ρ = identitytensornetworkstate(ComplexF64, g, s) ITensors.disable_warn_order() - use_gpu = true + use_gpu = false δβ = 0.01 hx = -3.1 diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index e1c0abd..412e82f 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -176,7 +176,10 @@ function evolve_bp(ρρ::BeliefPropagationCache, n::Int, nsteps::Int; β = 0, hx ec = prep_edges(n, g) apply_kwargs = (; maxdim = χ, cutoff = 1e-12) - + + if use_gpu + ρρ = CUDA.cu(ρρ) + end two_qubit_gates = [adapt(datatype(network(ρρ)), TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s)) for pair=vcat(ec...)] single_qubit_gates = [adapt(datatype(network(ρρ)), TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s)) for v=vertices(g)] From 5a30bb9a578532293f8d37af44bda469963606bc Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sat, 8 Nov 2025 12:12:47 -0500 Subject: [PATCH 18/36] time evolution kicked Ising --- .gitignore | 1 + cluster/time_evolution_tfim.jl | 77 ++++++++++++++++++++++++++++++++++ 2 files changed, 78 insertions(+) create mode 100644 cluster/time_evolution_tfim.jl diff --git a/.gitignore b/.gitignore index fec2ffa..6de8f23 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ *.jl.*.cov *.jl.cov *.jl.mem +*.~ /docs/Manifest*.toml /docs/build/ diff --git a/cluster/time_evolution_tfim.jl b/cluster/time_evolution_tfim.jl new file mode 100644 index 0000000..725bda0 --- /dev/null +++ b/cluster/time_evolution_tfim.jl @@ -0,0 +1,77 @@ +using TensorNetworkQuantumSimulator +const TN = TensorNetworkQuantumSimulator + +using NamedGraphs.NamedGraphGenerators: named_grid +using Statistics +using Base.Threads +using ITensors + +function get_columns(L) + pairs = [] + for i=1:L, j=1:L + push!(pairs, [[(i,j),(i,k)] for k=j+1:L]...) + end + pairs +end + +function main(L, χ, bmps_ranks; nl::Int=20,θh = 0) + ITensors.disable_warn_order() + g = named_grid((L,L); periodic =false) + nq = length(vertices(g)) + + #Define the gate parameters + J = pi / 4 + + layer = [] + + Rx_layer = [("Rx", [v], θh) for v in vertices(g)] + ec = edge_color(g, 4) + + Rzz_layer = [] + for edge_group in ec + append!(Rzz_layer, ("Rzz", pair, -2*J) for pair in edge_group) + end + + layer = vcat(Rx_layer, Rzz_layer) + + pairs = get_columns(L) + verts = reshape([(i,j) for i=1:L,j=1:L],L^2) + + # the initial state (all up, use Float 32 precision) + ψ0 = tensornetworkstate(ComplexF32, v -> "↑", g, "S=1/2") + + # max bond dimension for the TN + apply_kwargs = (maxdim = χ, cutoff = 1.0e-12, normalize_tensors = false) + + # create the BP cache representing the square of the tensor network + ψ_bpc = BeliefPropagationCache(ψ0) + + # an array to keep track of expectations taken via two methods + bpc_states = [] + bmps_expects_z = [zeros(ComplexF64,L,L,nl) for r=bmps_ranks] + bmps_expects_zz = [zeros(ComplexF64,length(pairs),nl) for r=bmps_ranks] + + # evolve! (First step takes long due to compilation) + for l in 1:nl + println("Layer $l") + + t1 = @timed ψ_bpc, errors = + apply_gates(layer, ψ_bpc; apply_kwargs, verbose = false) + + push!(bpc_states, copy(ψ_bpc)) + + #Boundary MPS expectation + ψ = network(ψ_bpc) + + @threads for r_i=1:length(bmps_ranks) + bmps_expects_zz[r_i][:,l] = expect(ψ, [(["Z","Z"], pair) for pair=pairs]; alg="boundarymps", mps_bond_dimension=bmps_ranks[r_i]) + z = expect(ψ, [("Z", [v]) for v=verts]; alg="boundarymps", mps_bond_dimension=bmps_ranks[r_i]) + bmps_expects_z[r_i][:,:,l] = reshape(z, (L,L)) + end + + println(" Took time: $(t1.time) [s]. Max bond dimension: $(maxvirtualdim(ψ_bpc))") + println(" Maximum Gate error for layer was $(maximum(errors))"); flush(stdout) + + end + return bpc_states, bmps_expects_z, bmps_expects_zz +end From 2c99aa09eab8b87944ebb7dba05ef83aa2afe56e Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sat, 8 Nov 2025 15:17:56 -0500 Subject: [PATCH 19/36] linkdim is now virtualdim --- PEPO_2d_ising.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index 412e82f..0168b7c 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -114,7 +114,7 @@ function expect_bmps(dat::Dict; obs = "X", MPS_message_rank::Int = 10, save_tag sqrtρ = dat["sqrtρ"][i] end @time expect_vals[:,i-start_i+1] = real.(TN.expect(sqrtρ, [(obs, [v]) for v=all_verts]; alg = "boundarymps", mps_bond_dimension = MPS_message_rank)) - save(DATA_DIR * "$(save_tag)L$(dat["L"])_χ$(dat["χ"])_D$(MPS_message_rank)_step$(dat["δβ"])_$(dat["β"][i]).jld2", Dict(obs=>expect_vals[:,i-start_i+1], "verts"=>all_verts, "hx"=>dat["hx"], "β"=>dat["β"][i], "χ"=>[dat["χ"],maxlinkdim(dat["sqrtρ"][i])], "mps_rank"=>MPS_message_rank, "δβ"=>dat["δβ"], "L"=>dat["L"])) + save(DATA_DIR * "$(save_tag)L$(dat["L"])_χ$(dat["χ"])_D$(MPS_message_rank)_step$(dat["δβ"])_$(dat["β"][i]).jld2", Dict(obs=>expect_vals[:,i-start_i+1], "verts"=>all_verts, "hx"=>dat["hx"], "β"=>dat["β"][i], "χ"=>[dat["χ"],maxvirtualdim(dat["sqrtρ"][i])], "mps_rank"=>MPS_message_rank, "δβ"=>dat["δβ"], "L"=>dat["L"])) flush(stdout) end all_verts, expect_vals From dd68ee8ced3e06862e5e29cde7bf2164f39aa699 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sat, 8 Nov 2025 15:59:36 -0500 Subject: [PATCH 20/36] update Project, add error tracking --- Project.toml | 14 +++++++++++++- cluster/time_evolution_tfim.jl | 9 +++++---- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 800537c..d06b647 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,16 @@ name = "TensorNetworkQuantumSimulator" uuid = "4de3b72a-362e-43dd-83ff-3f381eda9f9c" license = "MIT" -version = "0.2.2" authors = ["JoeyT1994 ", "MSRudolph ", "and contributors"] description = "A Julia package for quantum simulation with tensor networks of near-arbritrary topology." +version = "0.2.2" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" +DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a" +Debugger = "31a5f54b-26ea-5ae9-a837-f05ce5417438" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" EinExprs = "b1794770-133b-4de1-afb4-526377e9f4c5" Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8" @@ -21,6 +24,8 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" +MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" +NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605" NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19" PauliPropagation = "293282d5-3c99-4fb6-92d0-fd3280a19750" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" @@ -33,10 +38,14 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" +UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] Adapt = "4.3.0" +CUDA = "5.9.2" Combinatorics = "1.0.3" +DataGraphs = "0.2.9" +Debugger = "0.7.15" Dictionaries = "0.4" EinExprs = "0.6.4" Format = "1.3.7" @@ -50,6 +59,8 @@ KrylovKit = "0.10.2" LaTeXStrings = "1.4.0" LinearAlgebra = "1.11.0" LsqFit = "0.15.1" +MKL = "0.9.0" +NPZ = "0.4.3" NamedGraphs = "0.7" PauliPropagation = "0.3.0" Plots = "1.41.1" @@ -62,6 +73,7 @@ Statistics = "1.11.1" StatsBase = "0.34.4" TensorOperations = "5.2" TypeParameterAccessors = "0.3.10" +UnicodePlots = "3.8.1" julia = "1.10" [extras] diff --git a/cluster/time_evolution_tfim.jl b/cluster/time_evolution_tfim.jl index 725bda0..2fbebfb 100644 --- a/cluster/time_evolution_tfim.jl +++ b/cluster/time_evolution_tfim.jl @@ -47,7 +47,8 @@ function main(L, χ, bmps_ranks; nl::Int=20,θh = 0) ψ_bpc = BeliefPropagationCache(ψ0) # an array to keep track of expectations taken via two methods - bpc_states = [] + bpc_states = Array{BeliefPropagationCache}(undef, nl) + errs = zeros(nl) bmps_expects_z = [zeros(ComplexF64,L,L,nl) for r=bmps_ranks] bmps_expects_zz = [zeros(ComplexF64,length(pairs),nl) for r=bmps_ranks] @@ -58,7 +59,7 @@ function main(L, χ, bmps_ranks; nl::Int=20,θh = 0) t1 = @timed ψ_bpc, errors = apply_gates(layer, ψ_bpc; apply_kwargs, verbose = false) - push!(bpc_states, copy(ψ_bpc)) + bpc_states[l] = copy(ψ_bpc) #Boundary MPS expectation ψ = network(ψ_bpc) @@ -71,7 +72,7 @@ function main(L, χ, bmps_ranks; nl::Int=20,θh = 0) println(" Took time: $(t1.time) [s]. Max bond dimension: $(maxvirtualdim(ψ_bpc))") println(" Maximum Gate error for layer was $(maximum(errors))"); flush(stdout) - + errs[l] = maximum(errors) end - return bpc_states, bmps_expects_z, bmps_expects_zz + return bpc_states, bmps_expects_z, bmps_expects_zz, errs end From 402a5b829c20269dfc84c8958b9ddd9dd7614259 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sat, 8 Nov 2025 20:35:06 -0500 Subject: [PATCH 21/36] use_gpu option for boundary mps --- cluster/time_evolution_tfim.jl | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/cluster/time_evolution_tfim.jl b/cluster/time_evolution_tfim.jl index 2fbebfb..4175d0c 100644 --- a/cluster/time_evolution_tfim.jl +++ b/cluster/time_evolution_tfim.jl @@ -5,6 +5,7 @@ using NamedGraphs.NamedGraphGenerators: named_grid using Statistics using Base.Threads using ITensors +using CUDA function get_columns(L) pairs = [] @@ -14,7 +15,7 @@ function get_columns(L) pairs end -function main(L, χ, bmps_ranks; nl::Int=20,θh = 0) +function main(L, χ, bmps_ranks; nl::Int=20,θh = 0, use_gpu = true) ITensors.disable_warn_order() g = named_grid((L,L); periodic =false) nq = length(vertices(g)) @@ -62,11 +63,15 @@ function main(L, χ, bmps_ranks; nl::Int=20,θh = 0) bpc_states[l] = copy(ψ_bpc) #Boundary MPS expectation - ψ = network(ψ_bpc) - - @threads for r_i=1:length(bmps_ranks) - bmps_expects_zz[r_i][:,l] = expect(ψ, [(["Z","Z"], pair) for pair=pairs]; alg="boundarymps", mps_bond_dimension=bmps_ranks[r_i]) - z = expect(ψ, [("Z", [v]) for v=verts]; alg="boundarymps", mps_bond_dimension=bmps_ranks[r_i]) + if use_gpu + ψ = CUDA.cu(network(ψ_bpc)) + else + ψ = network(ψ_bpc) + end + + for r_i=1:length(bmps_ranks) + @time bmps_expects_zz[r_i][:,l] = expect(ψ, [(["Z","Z"], pair) for pair=pairs]; alg="boundarymps", mps_bond_dimension=bmps_ranks[r_i]) + @time z = expect(ψ, [("Z", [v]) for v=verts]; alg="boundarymps", mps_bond_dimension=bmps_ranks[r_i]) bmps_expects_z[r_i][:,:,l] = reshape(z, (L,L)) end From 2bd6ccbf25cb5d48d436a3d81c4385c7759ecaa6 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Mon, 10 Nov 2025 20:38:40 -0500 Subject: [PATCH 22/36] tracking BP message diffs --- PEPO_2d_ising.jl | 4 ++-- src/Apply/full_update.jl | 1 + .../abstractbeliefpropagationcache.jl | 22 ++++++++++++++----- 3 files changed, 20 insertions(+), 7 deletions(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index 0168b7c..4fd8ae3 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -176,7 +176,7 @@ function evolve_bp(ρρ::BeliefPropagationCache, n::Int, nsteps::Int; β = 0, hx ec = prep_edges(n, g) apply_kwargs = (; maxdim = χ, cutoff = 1e-12) - + bp_update_kwargs = (; maxiter = 50, tolerance = 1e-8, verbose = true) if use_gpu ρρ = CUDA.cu(ρρ) end @@ -187,7 +187,7 @@ function evolve_bp(ρρ::BeliefPropagationCache, n::Int, nsteps::Int; β = 0, hx layer = vcat(single_qubit_gates, two_qubit_gates, single_qubit_gates) for i in 1:nsteps - @time ρρ, errs = apply_gates(layer, ρρ; apply_kwargs, verbose = false) + @time ρρ, errs = apply_gates(layer, ρρ; bp_update_kwargs = bp_update_kwargs, apply_kwargs) β += δβ diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl index 5ec9a67..e5f331d 100644 --- a/src/Apply/full_update.jl +++ b/src/Apply/full_update.jl @@ -148,6 +148,7 @@ function optimise_p_q( q_cur, info = linsolve( M_p_tilde_partial, b_tilde_vec, q_cur; isposdef = envisposdef, ishermitian = false ) + println(info) end fend = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0 diff --git a/src/MessagePassing/abstractbeliefpropagationcache.jl b/src/MessagePassing/abstractbeliefpropagationcache.jl index 9780c94..46a593e 100644 --- a/src/MessagePassing/abstractbeliefpropagationcache.jl +++ b/src/MessagePassing/abstractbeliefpropagationcache.jl @@ -190,20 +190,32 @@ More generic interface for update, with default params """ function update(alg::Algorithm"bp", bpc::AbstractBeliefPropagationCache) compute_error = !isnothing(alg.kwargs.tolerance) + if !compute_error + println("ALERT, not computing error") + end if isnothing(alg.kwargs.maxiter) error("You need to specify a number of iterations for BP!") end bpc = copy(bpc) + diffs = zeros(alg.kwargs.maxiter) + tot_iter = alg.kwargs.maxiter for i in 1:alg.kwargs.maxiter diff = compute_error ? Ref(0.0) : nothing update_iteration!(alg, bpc, alg.kwargs.edge_sequence; (update_diff!) = diff) - if compute_error && (diff.x / length(alg.kwargs.edge_sequence)) <= alg.kwargs.tolerance - if alg.kwargs.verbose - println("BP converged to desired precision after $i iterations.") - end - break + if compute_error + diffs[i] = diff.x + if (diff.x / length(alg.kwargs.edge_sequence)) <= alg.kwargs.tolerance + if alg.kwargs.verbose + println("BP converged to desired precision after $i iterations.") + end + tot_iter = i + break + end end end + if compute_error && alg.kwargs.verbose + println("Diffs during message passing: $(diffs[1:tot_iter])") + end return bpc end From 48ba53994d01cf99e020aab66d6d0b531592dc85 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sat, 15 Nov 2025 23:25:29 -0500 Subject: [PATCH 23/36] coeffs and op_strings in free_energy --- .../abstractbeliefpropagationcache.jl | 44 +++++++++++++------ src/MessagePassing/beliefpropagationcache.jl | 5 +-- src/MessagePassing/loopcorrection.jl | 40 ++++++++++++----- 3 files changed, 61 insertions(+), 28 deletions(-) diff --git a/src/MessagePassing/abstractbeliefpropagationcache.jl b/src/MessagePassing/abstractbeliefpropagationcache.jl index 9780c94..ec5d70e 100644 --- a/src/MessagePassing/abstractbeliefpropagationcache.jl +++ b/src/MessagePassing/abstractbeliefpropagationcache.jl @@ -18,9 +18,9 @@ function rescale_vertices!( return not_implemented() end -function vertex_scalar(bp_cache::AbstractBeliefPropagationCache, vertex) +function vertex_scalar(bp_cache::AbstractBeliefPropagationCache, vertex; op_strings::Function = v->"I", coeffs::Function = v->1, use_epsilon::Bool = false) incoming_ms = incoming_messages(bp_cache, vertex) - state = bp_factors(bp_cache, vertex) + state = bp_factors(bp_cache, vertex; op_strings = op_strings, coeffs = coeffs, use_epsilon = use_epsilon) contract_list = [state; incoming_ms] sequence = contraction_sequence(contract_list; alg = "optimal") return contract(contract_list; sequence)[] @@ -115,8 +115,8 @@ function edge_scalars( return map(e -> edge_scalar(bp_cache, e; kwargs...), edges) end -function scalar_factors_quotient(bp_cache::AbstractBeliefPropagationCache) - return vertex_scalars(bp_cache), edge_scalars(bp_cache) +function scalar_factors_quotient(bp_cache::AbstractBeliefPropagationCache; kwargs...) + return vertex_scalars(bp_cache; kwargs...), edge_scalars(bp_cache) end function incoming_messages( @@ -239,21 +239,39 @@ function Adapt.adapt_structure(to, bpc::AbstractBeliefPropagationCache) return bpc end -function freenergy(bp_cache::AbstractBeliefPropagationCache) - numerator_terms, denominator_terms = scalar_factors_quotient(bp_cache) - if any(t -> real(t) < 0, numerator_terms) - numerator_terms = complex.(numerator_terms) +function adapt_complex(t) + if typeof(t)<:Hyper + return Hyper{ComplexF64}(t.value,t.epsilon1,t.epsilon2,t.epsilon12) + else + return complex(t) end - if any(t -> real(t) < 0, denominator_terms) - denominator_terms = complex.(denominator_terms) +end + +function get_real_part(t) + if typeof(t)<:Hyper + return real(t.value) + else + return real(t) + end +end + +function free_energy(bp_cache::AbstractBeliefPropagationCache; op_strings::Function = v->"I", coeffs::Function = v->1) + numerator_terms, denominator_terms = scalar_factors_quotient(bp_cache; op_strings = op_strings, coeffs = coeffs, use_epsilon = true) + + # Skip this piece for now + if any(t -> get_real_part(t) < 0, numerator_terms) + numerator_terms = adapt_complex.(numerator_terms) + end + if any(t -> get_real_part(t) < 0, denominator_terms) + denominator_terms = adapt_complex.(denominator_terms) end any(iszero, denominator_terms) && return -Inf return sum(log.(numerator_terms)) - sum(log.((denominator_terms))) end -function partitionfunction(bp_cache::AbstractBeliefPropagationCache) - return exp(freenergy(bp_cache)) +function partitionfunction(bp_cache::AbstractBeliefPropagationCache; op_strings::Function = v->"I", coeffs::Function = v->1) + return exp(free_energy(bp_cache; op_strings = op_strings, coeffs = coeffs)) end function rescale_messages!(bp_cache::AbstractBeliefPropagationCache, edge::AbstractEdge) @@ -278,4 +296,4 @@ function rescale(bpc::AbstractBeliefPropagationCache, args...; kwargs...) bpc = copy(bpc) rescale!(bpc, args...; kwargs...) return bpc -end +end \ No newline at end of file diff --git a/src/MessagePassing/beliefpropagationcache.jl b/src/MessagePassing/beliefpropagationcache.jl index 7da4c33..b42bc85 100644 --- a/src/MessagePassing/beliefpropagationcache.jl +++ b/src/MessagePassing/beliefpropagationcache.jl @@ -86,13 +86,12 @@ end function rescale_vertices!( bpc::BeliefPropagationCache, - vertices::Vector - ) + vertices::Vector; kwargs...) tn = network(bpc) for v in vertices vn = vertex_scalar(bpc, v) - s = isreal(vn) ? sign(vn) : one(vn) + s = isreal(vn) ? sign(vn) : one(vn) if tn isa TensorNetworkState setindex_preserve!(tn, tn[v] * s * inv(sqrt(vn)), v) elseif tn isa TensorNetwork diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl index bbbe7a4..266b03c 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/loopcorrection.jl @@ -1,4 +1,5 @@ using NamedGraphs.GraphsExtensions: boundary_edges +using HyperDualNumbers function loopcorrected_partitionfunction( bp_cache::BeliefPropagationCache, @@ -63,34 +64,49 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) end #Get the all edges incident to the region specified by the vector of edges passed +# Optionally pass in a vector of vertices in the region, to handle the corner case where the region is just one vertex. function NamedGraphs.GraphsExtensions.boundary_edges( bpc::BeliefPropagationCache, - es::Vector{<:NamedEdge}, + es::Vector{<:NamedEdge}; + vs::Vector = [] ) - vs = unique(vcat(src.(es), dst.(es))) + + if isempty(vs) + vs = unique(vcat(src.(es), dst.(es))) + end bpes = NamedEdge[] for v in vs - incoming_es = NamedGraphs.GraphsExtensions.boundary_edges(bpc, [v]; dir = :in) + incoming_es = 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) +#Compute the contraction of the bp configuration specified by the edge induced subgraph eg. Insert I + epsilon O on up to two sites. +function weight(bpc::BeliefPropagationCache, eg; project_out::Bool = true, op_strings::Function = v->"I", coeffs::Function = v->1, rescales = Dictionary(1 for v=vertices(eg))) vs = collect(vertices(eg)) es = collect(edges(eg)) - bpc, antiprojectors = sim_edgeinduced_subgraph(bpc, eg) + + if project_out + bpc, antiprojectors = sim_edgeinduced_subgraph(bpc, eg) + end incoming_ms = - ITensor[message(bpc, e) for e in boundary_edges(bpc, es)] - local_tensors = reduce(vcat, [bp_factors(bpc, v) for v in vs]) - ts = [incoming_ms; local_tensors; antiprojectors] + ITensor[message(bpc, e) for e in boundary_edges(bpc, es; vs = vs)] + local_tensors = reduce(vcat, bp_factors(bpc, vs; op_strings = op_strings, coeffs = coeffs, use_epsilon = true)) + + if project_out + ts = [incoming_ms; local_tensors; antiprojectors] + else + ts = [incoming_ms; local_tensors] + end + seq = any(hasqns.(ts)) ? contraction_sequence(ts; alg = "optimal") : contraction_sequence(ts; alg = "einexpr", optimizer = Greedy()) - return contract(ts; sequence = seq)[] + output = contract(ts; sequence = seq)[] + return output / prod([rescales[v] for v=vs]) end #Vectorized version of weight -function weights(bpc::BeliefPropagationCache, egs) - return [weight(bpc, eg) for eg in egs] +function weights(bpc::BeliefPropagationCache, egs; kwargs...) + return [weight(bpc, eg; kwargs...) for eg in egs] end From ef1d6231285878e85454c0b00a4edcc25e17368c Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sat, 15 Nov 2025 23:26:44 -0500 Subject: [PATCH 24/36] norm_factors with kwargs --- src/TensorNetworks/tensornetworkstate.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl index 410278d..0483faa 100644 --- a/src/TensorNetworks/tensornetworkstate.jl +++ b/src/TensorNetworks/tensornetworkstate.jl @@ -1,4 +1,5 @@ using ITensors: random_itensor +using HyperDualNumbers #TODO: Make this show() nicely. struct TensorNetworkState{V} <: AbstractTensorNetwork{V} @@ -39,7 +40,7 @@ function Base.setindex!(tns::TensorNetworkState, value::ITensor, v) return tns end -function norm_factors(tns::TensorNetworkState, verts::Vector; op_strings::Function = v -> "I") +function norm_factors(tns::TensorNetworkState, verts::Vector; op_strings::Function = v -> "I", coeffs::Function = v->1, use_epsilon::Bool = false) factors = ITensor[] for v in verts sinds = siteinds(tns, v) @@ -47,17 +48,18 @@ function norm_factors(tns::TensorNetworkState, verts::Vector; op_strings::Functi tnv_dag = dag(prime(tnv)) if op_strings(v) == "I" tnv_dag = replaceinds(tnv_dag, prime.(sinds), sinds) - append!(factors, ITensor[tnv, tnv_dag]) + append!(factors, ITensor[coeffs(v) * tnv, tnv_dag]) else - op = adapt(datatype(tnv))(ITensors.op(op_strings(v), only(sinds))) + op = use_epsilon ? Hyper(1,0,0,0) * ITensors.op("I", only(sinds)) + coeffs(v) * ITensors.op(op_strings(v), only(sinds)) : coeffs(v) * ITensors.op(op_strings(v), only(sinds)) append!(factors, ITensor[tnv, tnv_dag, op]) end end return factors end -norm_factors(tns::TensorNetworkState, v) = norm_factors(tns, [v]) -bp_factors(tns::TensorNetworkState, v) = norm_factors(tns, v) +norm_factors(tns::TensorNetworkState, v; kwargs...) = norm_factors(tns, [v]; kwargs...) +bp_factors(tns::TensorNetworkState, v; kwargs...) = norm_factors(tns, v; kwargs...) +bp_factors(tns::TensorNetworkState, verts::Vector; kwargs...) = norm_factors(tns, verts; kwargs...) function default_message(tns::TensorNetworkState, edge::AbstractEdge) linds = virtualinds(tns, edge) From 7f7346391f9e6dd992810a6de0eb251d04b19bc3 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sat, 15 Nov 2025 23:29:58 -0500 Subject: [PATCH 25/36] cluster expansions --- Project.toml | 6 + cluster/expect-corrected.jl | 285 ++++--------- src/MessagePassing/clustercorrections.jl | 235 +++++++++++ .../cumulant-clustercorrections.jl | 394 ++++++++++++++++++ src/TensorNetworkQuantumSimulator.jl | 3 + src/graph_enumeration.jl | 108 +++++ 6 files changed, 818 insertions(+), 213 deletions(-) create mode 100644 src/MessagePassing/clustercorrections.jl create mode 100644 src/MessagePassing/cumulant-clustercorrections.jl create mode 100644 src/graph_enumeration.jl diff --git a/Project.toml b/Project.toml index 2ed5a3b..25a06b7 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,7 @@ Debugger = "31a5f54b-26ea-5ae9-a837-f05ce5417438" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" EinExprs = "b1794770-133b-4de1-afb4-526377e9f4c5" Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8" +GraphIO = "aa1b3936-2fda-51b9-ab35-c553d3a640a2" GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97" @@ -21,6 +22,7 @@ IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +JuliaSyntaxHighlighting = "ac6e5ff7-fb65-4e79-a425-ec3bc9c03011" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -30,6 +32,7 @@ NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605" NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19" PauliPropagation = "293282d5-3c99-4fb6-92d0-fd3280a19750" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SimpleGraphAlgorithms = "41400c72-0c58-5c16-8579-4ecbce768449" @@ -50,6 +53,7 @@ Debugger = "0.7.15" Dictionaries = "0.4" EinExprs = "0.6.4" Format = "1.3.7" +GraphIO = "0.7.1" GraphRecipes = "0.5.13" Graphs = "1.8.0" HyperDualNumbers = "4.0.10" @@ -57,6 +61,7 @@ IJulia = "1.32.1" ITensorMPS = "0.3.17" ITensors = "0.9" JLD2 = "0.6.2" +JuliaSyntaxHighlighting = "0.1.0" KrylovKit = "0.10.2" LaTeXStrings = "1.4.0" LinearAlgebra = "1.11.0" @@ -66,6 +71,7 @@ NPZ = "0.4.3" NamedGraphs = "0.7" PauliPropagation = "0.3.0" Plots = "1.41.1" +ProgressMeter = "1.11.0" PythonPlot = "1.0.6" Revise = "3.8.0" SimpleGraphAlgorithms = "0.6.0" diff --git a/cluster/expect-corrected.jl b/cluster/expect-corrected.jl index 137eb90..beb28fb 100644 --- a/cluster/expect-corrected.jl +++ b/cluster/expect-corrected.jl @@ -13,252 +13,111 @@ using Adapt: adapt using Dictionaries using ITensors: inds, onehot, dag, commonind, op -# Cluster/loop/CC weights -function cluster_weights(bpc::BeliefPropagationCache, clusters::Vector, egs::Vector{<:AbstractNamedGraph}, interaction_graph; rescale::Bool = true) - logZbp = logscalar(bpc) - if rescale - bpc = TN.rescale(bpc) +function prep_insertions(obs) + if isnothing(obs) + return (coeffs = identity, op_strings = v->"I") + end + op_strings, verts, _ = TN.collectobservable(obs) + @assert length(verts) <= 2 + + function hyper_coeff(v) + if v==verts[1] + return Hyper(0,1,0,0) + elseif length(verts)==2 && v==verts[2] + return Hyper(0,0,1,0) + else + return 1 + end + end + + function insertion_operator(v) + if v==verts[1] + return op_strings[1] + elseif length(verts)==2 && v==verts[2] + return op_strings[2] + else + return "I" + end end + return (coeffs = hyper_coeff, op_strings = insertion_operator) +end + +# Cluster/cluster cumulant weights +function cluster_weights(bpc::BeliefPropagationCache, clusters::Vector, egs::Vector{<:AbstractNamedGraph}, interaction_graph; obs = nothing) + kwargs = prep_insertions(obs) + + logZbp = TN.free_energy(bpc; kwargs...) + println("MADE IT HERE"); flush(stdout) isempty(egs) && return [0], [[logZbp]], [[1]] circuit_lengths = sort(unique([c.weight for c=clusters])) + # Rescale the messages, but deal with the vertices separately + TN.rescale_messages!(bpc) + vns = Dictionary(TN.vertex_scalar(bpc, v; use_epsilon = true, kwargs...) for v=vertices(network(bpc).tensornetwork.graph)) + # calculate weight of each generalized loop first - wts = TN.weights(bpc, egs) - + wts = TN.weights(bpc, egs; rescales = vns, kwargs...) + logZs = Array{Array}(undef, length(circuit_lengths) + 1) logZs[1] = [logZbp] coeffs = Array{Array}(undef, length(circuit_lengths) + 1) coeffs[1] = [1] - + println("MADE IT HERE TOO"); flush(stdout) # now calculate contribution to logZ from each cluster for (cl_i, cl)=enumerate(circuit_lengths) clusters_cl = filter(c->c.weight==cl, clusters) logZs[cl_i + 1] = [prod([prod(fill(wts[l],c.multiplicities[l])) for l=c.loop_ids]) for c=clusters_cl] - # weird bug in HyperDualNumbers where this doesn't work... - # sum_ws = sum([TN.ursell_function(c, interaction_graph) * prod([wts[l]^c.multiplicities[l] for l=c.loop_ids]) for c=clusters_cl]) coeffs[cl_i + 1] = [TN.ursell_function(c, interaction_graph) for c=clusters_cl] end return vcat([0],circuit_lengths), logZs, coeffs end -function cc_weights(bpc::BeliefPropagationCache, regions::Vector, counting_nums::Dict; rescale::Bool=true) - logZbp = logscalar(bpc) - if rescale - bpc = TN.rescale(bpc) - end +function cc_weights(bpc::BeliefPropagationCache, regions::Vector, counting_nums::Dict; obs = nothing, rescale::Bool = false) + kwargs = prep_insertions(obs) + use_g = findall(gg->counting_nums[gg] != 0, regions) - egs = [induced_subgraph(bpc.partitioned_tensornetwork.partitioned_graph, gg)[1] for gg=regions[use_g]] + egs = [induced_subgraph(network(bpc).tensornetwork.graph, gg)[1] for gg=regions[use_g]] - isempty(egs) && return [0], [[logZbp]], [[1]] - - # calculate weight of each generalized loop first - wts = TN.full_weights(bpc, egs) - - return logZbp, log.(wts), [counting_nums[gg] for gg=regions[use_g]] -end - - -# Loop series expansion + isempty(egs) && return logZbp, [], [] -function terms_to_scalar(numerator_numerator_terms, numerator_denominator_terms, denominator_numerator_terms, denominator_denominator_terms) - return exp(sum(log.(numerator_numerator_terms)) - sum(log.(numerator_denominator_terms)) - sum(log.(denominator_numerator_terms)) + sum(log.(denominator_denominator_terms))) -end - -#Project spins on sites v1 and v2 to v1_val (1 = up, 2 = down) and v2_val -function project!(ψIψ::BeliefPropagationCache, v1, v2, v1_val::Int64 = 1, v2_val::Int64=1) - s1 = only(inds(only(ITN.factors(ψIψ, [(v1, "operator")])); plev = 0)) - s2 = only(inds(only(ITN.factors(ψIψ, [(v2, "operator")])); plev = 0)) - ITensorNetworks.@preserve_graph ψIψ[(v1, "operator")] = onehot(s1 => v1_val) * dag(onehot(s1' => v1_val)) - ITensorNetworks.@preserve_graph ψIψ[(v2, "operator")] = onehot(s2 => v2_val) * dag(onehot(s2' => v2_val)) -end - -#Log scalar contraction of bpc -function logscalar(bpc::BeliefPropagationCache) - nums, denoms = TN.scalar_factors_quotient(bpc) - return sum(log.(nums)) - sum(log.(denoms)) -end - -function cumulative_weights(bpc::BeliefPropagationCache, egs::Vector{<:AbstractNamedGraph}) - isempty(egs) && [1] - circuit_lengths = sort(unique(length.(edges.(egs)))) - outs = [] - for cl in circuit_lengths - egs_cl = filter(eg -> length(edges(eg)) == cl, egs) - sum_ws = sum(TN.weights(bpc, egs_cl)) - push!(outs, sum_ws) - end - - # first one is the BP contribution - outs = vcat([1.0], outs) - return cumsum(outs) -end - -function compute_ps!(ψIψ::BeliefPropagationCache, v1, v2, v1_val, v2_val, egs::Vector{<:AbstractNamedGraph}; track::Bool = true, cache_update_kwargs...) - project!(ψIψ, v1, v2, v1_val, v2_val) - - if track - ψIψ, bp_diffs = updatecache(ψIψ; track = true, cache_update_kwargs...) + # Rescale the messages, but deal with the vertices separately + if rescale + TN.rescale_messages!(bpc) + vns = Dictionary(TN.vertex_scalar(bpc, v; use_epsilon = true, kwargs...) for v=vertices(network(bpc).tensornetwork.graph)) else - ψIψ = updatecache(ψIψ; track = false, cache_update_kwargs...) - bp_diffs = [] + vns = Dictionary(1 for v=vertices(network(bpc).tensornetwork.graph)) end + + # calculate weight of each cluster first + wts = TN.weights(bpc, egs; rescales = vns, project_out = false, kwargs...) - p_bp = exp(logscalar(ψIψ)) - - ψIψ = ITensorNetworks.rescale(ψIψ) - cfes = cumulative_weights(ψIψ, egs) - - return [p_bp*cfe for cfe in cfes], bp_diffs -end - -function scalar_cumul_loop(bp_cache::BeliefPropagationCache, egs) - zbp = exp(logscalar(bp_cache)) - bp_cache = ITN.rescale(bp_cache) - - cfes = cumulative_weights(bp_cache, egs) - - return [zbp * cfe for cfe=cfes] -end - -function zz_bp_loopcorrect_connected(ψIψ::BeliefPropagationCache, verts, egs::Vector{<:AbstractNamedGraph}; cache_update_kwargs...) - z_expects = [z_bp_loopcorrect(ψIψ, v, egs; cache_update_kwargs...)[1] for v=verts] - zz_expect = zz_bp_loopcorrect(ψIψ, verts...,egs; cache_update_kwargs...) - - zz = (zz_expect[1,1] .+ zz_expect[2,2] .- (zz_expect[1,2] .+ zz_expect[2,1])) ./sum(zz_expect) - vcat([zz], [(z[1] .- z[2]) ./ (z[1] .+ z[2]) for z=z_expects]) + return log.(wts), [counting_nums[gg] for gg=regions[use_g]] end -# compute z expectation value, but return cumulative weights rather than just total loop-corrected value -function z_bp_loopcorrect(ψIψ::BeliefPropagationCache, v, egs::Vector{<:AbstractNamedGraph}; cache_update_kwargs...) - ψUpψ = TN.insert_observable(ψIψ, ("Proj0", [v])) - ψUpψ, diffs_up = updatecache(ψUpψ; track = true, cache_update_kwargs...) - p_up = scalar_cumul_loop(ψUpψ, egs) - ψDnψ = TN.insert_observable(ψIψ, ("Proj1", [v])) - ψDnψ, diffs_dn = updatecache(ψDnψ; track = true, cache_update_kwargs...) - p_dn = scalar_cumul_loop(ψDnψ, egs) - - denominator = scalar_cumul_loop(ψIψ, egs) - return [p_up, p_dn, denominator], [diffs_up, diffs_dn] -end - -#Compute zz with loop correction and all bells and whistles -function zz_bp_loopcorrect(ψIψ::BeliefPropagationCache, v1, v2, egs::Vector{<:AbstractNamedGraph}; cache_update_kwargs...) - probs = [compute_ps!(copy(ψIψ), v1, v2, i, j, egs; track = false, cache_update_kwargs...)[1] for i=1:2, j=1:2] -end - -function cluster_twopoint_correlator(ψIψ::BeliefPropagationCache, obs, max_weight::Int) - ng = ITN.partitioned_graph(ψIψ) +# onepoint or twopoint connected correlation function, using cluster cumulant expansion +function cc_correlation(bpc::BeliefPropagationCache, regions::Vector, counting_nums::Dict, obs) + logZs, cnums = cc_weights(bpc, regions, counting_nums; obs = obs) op_strings, verts, _ = TN.collectobservable(obs) - @assert length(verts)<=2 - clusters, egs, interaction_graph = enumerate_connected_clusters_twopoint(ng, verts..., max_weight) - - cluster_twopoint_correlator(ψIψ, obs, clusters, egs, interaction_graph) -end - -function cluster_twopoint_correlator(ψIψ::BeliefPropagationCache, obs, clusters, egs, interaction_graph) - # rescale BEFORE inserting operator - ψIψ = TN.rescale(ψIψ) - op_strings, verts, _ = TN.collectobservable(obs) - @assert length(verts)<=2 - - ψIψ_vs = [ψIψ[(v, "operator")] for v in verts] - sinds = - [commonind(ψIψ[(v, "ket")], ψIψ_vs[i]) for (i, v) in enumerate(verts)] - - - coeffs = [Hyper(0,1,0,0),Hyper(0,0,1,0)] - operators = [ψIψ[(v, "operator")].tensor[1,1] * adapt(typeof(ψIψ[(v, "operator")]))(Hyper(1,0,0,0) * op("I", sinds[i]) + coeffs[i] * op(op_strings[i], sinds[i])) for (i, v) in enumerate(verts)] - ψOψ = ITN.update_factors(ψIψ, Dictionary([(v, "operator") for v in verts], operators)) - - cluster_weights(ψOψ, clusters, egs, interaction_graph; rescale = false) -end - -function cluster_onepoint(ψIψ::BeliefPropagationCache, obs, clusters, egs, interaction_graph) - ψOψ, logZbp = prep_insertion(ψIψ, obs) - lens, logZs, coeffs = cluster_weights(ψOψ, clusters, egs, interaction_graph; rescale = false) - logZs[1] .+= logZbp - lens, logZs, coeffs -end - -# copied over from ITensorNetworks, - -# Also return the amount rescaled by -function rescale_partitions_norms( - bpc::ITN.AbstractBeliefPropagationCache, - partitions::Vector; - verts::Vector=vertices(bpc, partitions), -) - bpc = copy(bpc) - tn = ITN.tensornetwork(bpc) - - # not sure why this is done in two steps - norms = map(v -> inv(ITN.norm(tn[v])), verts) - ITN.scale!(bpc, Dictionary(verts, norms)) - - vertices_weights = Dictionary() - for pv in partitions - pv_vs = filter(v -> v ∈ verts, vertices(bpc, pv)) - isempty(pv_vs) && continue - - vn = ITN.region_scalar(bpc, pv) - s = one(vn)#sign(vn) #isreal(vn) ? sign(vn) : one(vn) - vn = s * vn^(-inv(oftype(vn, length(pv_vs)))) - set!(vertices_weights, first(pv_vs), s*vn) - for v in pv_vs[2:length(pv_vs)] - set!(vertices_weights, v, vn) + if length(verts)==1 + return sum(logZs .* cnums).epsilon1 + else + return sum(logZs .* cnums).epsilon12 end - end - - ITN.scale!(bpc, vertices_weights) - - return bpc, Dictionary(verts,norms), vertices_weights -end - -# cluster cumulant expansion - -function cc_twopoint_correlator(ψIψ::BeliefPropagationCache, obs, regions::Vector, counting_nums::Dict) - # rescale BEFORE inserting operator - ψIψ = TN.rescale(ψIψ) - op_strings, verts, _ = TN.collectobservable(obs) - @assert length(verts)==2 - - ψIψ_vs = [ψIψ[(v, "operator")] for v in verts] - sinds = - [commonind(ψIψ[(v, "ket")], ψIψ_vs[i]) for (i, v) in enumerate(verts)] - - coeffs = [Hyper(0,1,0,0),Hyper(0,0,1,0)] - operators = [ψIψ[(v, "operator")].tensor[1,1] * adapt(typeof(ψIψ[(v, "operator")]))(Hyper(1,0,0,0) * op("I", sinds[i]) + coeffs[i] * op(op_strings[i], sinds[i])) for (i, v) in enumerate(verts)] - ψOψ = ITN.update_factors(ψIψ, Dictionary([(v, "operator") for v in verts], operators)) - - cc_weights(ψOψ, regions, counting_nums; rescale=false) -end - -function cc_onepoint(ψIψ::BeliefPropagationCache, obs, regions::Vector, counting_nums::Dict) - ψOψ, logZbp = prep_insertion(ψIψ, obs) - lz, logZs, coeffs = cc_weights(ψOψ, regions, counting_nums; rescale=false) - logZbp + lz, logZs, coeffs end -function prep_insertion(ψIψ::BeliefPropagationCache, obs) - # rescale BEFORE inserting operator - ψIψ = ITN.rescale_messages(ψIψ) +# onepoint or twopoint connected correlation function, using cluster expansion +function cluster_correlation(bpc::BeliefPropagationCache, clusters::Vector, egs::Vector, interaction_graph, obs) + cluster_wts, logZs, ursells = cluster_weights(bpc, clusters, egs, interaction_graph; obs = obs) op_strings, verts, _ = TN.collectobservable(obs) - @assert length(verts)<=2 - - ψIψ_vs = [ψIψ[(v, "operator")] for v in verts] - sinds = - [commonind(ψIψ[(v, "ket")], ψIψ_vs[i]) for (i, v) in enumerate(verts)] - - coeffs = [Hyper(0,1,0,0),Hyper(0,0,1,0)] - operators = [adapt(typeof(ψIψ[(v, "operator")]))(Hyper(1,0,0,0) * op("I", sinds[i]) + coeffs[i] * op(op_strings[i], sinds[i])) for (i, v) in enumerate(verts)] - ψOψ = ITN.update_factors(ψIψ, Dictionary([(v, "operator") for v in verts], operators)) - - ψOψ, norms, vertices_weights = rescale_partitions_norms(ψOψ, collect(ITN.partitions(ψIψ))) - return ψOψ, sum([log(1/(norms[(v,tag)] * vertices_weights[(v,tag)])) for v=verts, tag = ["bra", "ket", "operator"]]) -end - + cumul_dat = cumsum([sum([logZs[i][j] * ursells[i][j] for j=1:length(logZs[i])]) for i=1:length(logZs)]) + if length(verts)==1 + return cluster_wts, [d.epsilon1 for d=cumul_dat] + else + return cluster_wts, [d.epsilon12 for d=cumul_dat] + end +end \ No newline at end of file diff --git a/src/MessagePassing/clustercorrections.jl b/src/MessagePassing/clustercorrections.jl new file mode 100644 index 0000000..e4ff304 --- /dev/null +++ b/src/MessagePassing/clustercorrections.jl @@ -0,0 +1,235 @@ + +using NamedGraphs +using NamedGraphs: AbstractNamedGraph +using ProgressMeter + +include("../graph_enumeration.jl") + +struct Cluster + loop_ids::Vector{Int} + multiplicities::Dict{Int, Int} + weight::Int + total_loops::Int +end + +struct Loop + vertices::Vector{Int} + edges::Vector{Tuple{Int,Int}} + weight::Int +end + + +function canonical_cluster_signature(cluster::Cluster) + items = [(loop_id, cluster.multiplicities[loop_id]) for loop_id in sort(cluster.loop_ids)] + return (tuple(items...), cluster.weight) +end + +function build_interaction_graph(loops::Vector{Loop}) + """Build interaction graph with optimizations for speed.""" + interaction_graph = Dict{Int, Vector{Int}}() + n_loops = length(loops) + + println(" Building optimized interaction graph for $n_loops loops...") + flush(stdout) + progress = Progress(n_loops, dt=0.1, desc="Building graph: ", color=:green, barlen=50) + + # Optimization 1: Pre-compute vertex sets once + vertex_sets = [Set(loop.vertices) for loop in loops] + + # Optimization 2: Build vertex-to-loops mapping for faster lookup + vertex_to_loops = Dict{Int, Vector{Int}}() + for (i, loop) in enumerate(loops) + for vertex in loop.vertices + if !haskey(vertex_to_loops, vertex) + vertex_to_loops[vertex] = Int[] + end + push!(vertex_to_loops[vertex], i) + end + end + + for i in 1:n_loops + interaction_graph[i] = unique(vcat([vertex_to_loops[v] for v=loops[i].vertices]...)) + + next!(progress) + end + + return interaction_graph +end + +""" +Enumerate connected clusters using DFS starting from loops supported on target site. +Connectivity is guaranteed by growing through the interaction graph. +Courtesy of Frank Zhang and Siddhant Midha +""" +function dfs_enumerate_clusters_from_supported(all_loops::Vector{Loop}, supported_loop_ids::Vector{Int}, max_weight::Int, interaction_graph::Dict{Int, Vector{Int}}; verbose::Bool = false) + clusters = Cluster[] + seen_clusters = Set{Tuple}() + cluster_count = 0 + + verbose && println(" Starting DFS cluster enumeration...") + verbose && println(" Supported loops: $(length(supported_loop_ids)), Max weight: $max_weight") + + # Progress tracking + last_report_time = time() + last_cluster_count = 0 + + # DFS to grow clusters starting from each supported loop + function dfs_grow_cluster(current_cluster::Vector{Int}, current_weight::Int, + has_supported::Bool) + + # If we've found a valid cluster (has supported loop), record it + if has_supported && current_weight >= 1 + # Create cluster with multiplicities + multiplicities = Dict{Int, Int}() + for loop_id in current_cluster + multiplicities[loop_id] = get(multiplicities, loop_id, 0) + 1 + end + + cluster = Cluster( + collect(keys(multiplicities)), + multiplicities, + current_weight, + length(current_cluster) + ) + + # Avoid duplicates using canonical signature + signature = canonical_cluster_signature(cluster) + if !(signature in seen_clusters) + push!(seen_clusters, signature) + push!(clusters, cluster) + cluster_count += 1 + + # Progress reporting every 2 seconds + current_time = time() + if current_time - last_report_time >= 2.0 + new_clusters = cluster_count - last_cluster_count + verbose && println(" Found $cluster_count clusters (+$new_clusters in last 2s)") + last_report_time = current_time + last_cluster_count = cluster_count + end + end + end + + # Stop if we've reached max weight + if current_weight >= max_weight + return + end + + # Find candidate loops to add (adjacent loops or multiplicities) + candidate_loops = Set{Int}() + + if isempty(current_cluster) + # Start with supported loops only + for loop_id in supported_loop_ids + if all_loops[loop_id].weight <= max_weight - current_weight + push!(candidate_loops, loop_id) + end + end + else + # Add loops connected to current cluster via interaction graph + for loop_id in current_cluster + # Add connected loops (touching loops) + for neighbor_id in get(interaction_graph, loop_id, Int[]) + if all_loops[neighbor_id].weight <= max_weight - current_weight + push!(candidate_loops, neighbor_id) + end + end + # Allow multiplicity increases (same loop added again) + if all_loops[loop_id].weight <= max_weight - current_weight + push!(candidate_loops, loop_id) + end + end + end + + # Try each candidate loop + for loop_id in candidate_loops + loop_weight = all_loops[loop_id].weight + new_weight = current_weight + loop_weight + + if new_weight <= max_weight + new_cluster = copy(current_cluster) + push!(new_cluster, loop_id) + new_has_supported = has_supported || (loop_id in supported_loop_ids) + + # Continue DFS (connectivity guaranteed by interaction graph) + dfs_grow_cluster(new_cluster, new_weight, new_has_supported) + end + end + end + + # Start DFS with empty cluster + dfs_grow_cluster(Int[], 0, false) + + verbose && println(" DFS enumeration completed: $cluster_count total clusters found") + return clusters +end + +""" +Build all clusters on named graph ng, up to a given weight. Optionally, must be supported on the vertices must_contain, in which case those vertices can be leaves + +This is overkill as it finds ALL subgraphs first, but my other implementation had bugs +""" +function enumerate_clusters(ng::NamedGraph, max_weight::Int; min_v::Int = 4, triangle_free::Bool = true, must_contain = [], min_deg::Int = 2, verbose::Bool = false) + g = ng.position_graph + ordered_indices = ng.vertices.ordered_indices + + verbose && println("Step 1: find embedded generalized loops") + subgraphs = generate_embedded_graphs(g, max_weight; min_v = min_v, triangle_free = triangle_free, min_deg = min_deg, leaf_vertices = [ng.vertices.index_positions[v] for v=must_contain]) + + # convert into form of LoopEnumeration.jl + loops = [Loop(sort(unique(vcat([[e[1],e[2]] for e=subg]...))), subg, length(subg)) for subg=subgraphs] + + verbose && println("Found $(length(loops)) loops") + + verbose && println("Step 2: Building interaction graph...") + interaction_graph = build_interaction_graph(loops) + + # DFS cluster enumeration + verbose && println("Step 3: DFS cluster enumeration...") + if isempty(must_contain) + supported_loops = [1:length(loops);] + else + supported_loops = findall(el->all(l->ng.vertices.index_positions[l] in el.vertices, must_contain), loops) + verbose && println("$(length(supported_loops)) supported...") + end + + all_clusters = dfs_enumerate_clusters_from_supported(loops, supported_loops, max_weight, interaction_graph, verbose = verbose) + verbose && println("Found $(length(all_clusters)) connected clusters") + + # converting loops into NamedGraphs, for use in tensor_weights + return all_clusters, [generalized_loop_named(l, ordered_indices) for l=loops], interaction_graph + +end + +""" +Convert from Loop into NamedGraph +""" +function generalized_loop_named(loop::Loop, ordered_indices) + g = NamedGraph(ordered_indices[loop.vertices]) + for e=loop.edges + add_edge!(g, ordered_indices[e[1]], ordered_indices[e[2]]) + end + g +end + +function ursell_function(cluster::Cluster, adj::Dict) + """ + Compute the Ursell function φ(W) for a connected cluster W. + """ + total_loops = cluster.total_loops + + if length(cluster.loop_ids) > 2 + for i=1:length(cluster.loop_ids) + for j=1:i-1 + if !(cluster.loop_ids[i] in adj[cluster.loop_ids[j]]) + error("Only implemented clusters corresponding to complete graphs for now, but got $(cluster)") + end + end + end + end + + no_vertices = sum(values(cluster.multiplicities)) + denominator = prod(factorial.(values(cluster.multiplicities))) + numerator = (-1)^(no_vertices - 1)* factorial(no_vertices - 1) + return numerator / denominator +end \ No newline at end of file diff --git a/src/MessagePassing/cumulant-clustercorrections.jl b/src/MessagePassing/cumulant-clustercorrections.jl new file mode 100644 index 0000000..7c1b3bd --- /dev/null +++ b/src/MessagePassing/cumulant-clustercorrections.jl @@ -0,0 +1,394 @@ +using Graphs: nv, induced_subgraph, is_connected, connected_components +using NamedGraphs +using NamedGraphs: AbstractNamedGraph, position_graph +using ProgressMeter +using Dictionaries + +include("../graph_enumeration.jl") + +const RegionKey = NTuple{N,Int} where {N} + +""" + to_key(vs::AbstractVector{Int})::RegionKey +Convert a collection of vertex IDs to a canonical, sorted tuple key. +""" +function to_key(vs::AbstractVector{Int})::RegionKey + return Tuple(sort!(collect(vs))) +end + +""" + key_intersection(a::RegionKey, b::RegionKey)::Vector{Int} +Return the vertex list of the intersection of two region keys. +""" +function key_intersection(a::RegionKey, b::RegionKey)::Vector{Int} + sa = Set(a); sb = Set(b) + return sort!(collect(intersect(sa, sb))) +end + +""" + is_loopful(g::SimpleGraph, key::RegionKey)::Bool +Check if the induced subgraph on `key` is connected and has at least one cycle. +For connected component `H`, loopfulness is `ne(H) - nv(H) + 1 > 0`. +""" +function is_loopful(g::SimpleGraph, key::RegionKey)::Bool + if length(key) < 3 + return false + end + vs = collect(key) + h, _ = induced_subgraph(g, vs) + # ensure connected (defensive; we try to keep regions connected elsewhere) + if nv(h) == 0 + return false + end + if !is_connected(h) + return false + end + return ne(h) - nv(h) + 1 > 0 +end + +""" + is_proper_subset(a::RegionKey, b::RegionKey)::Bool +Return true if region `a` is a proper subset of `b`. +""" +function is_proper_subset(a::RegionKey, b::RegionKey)::Bool + if length(a) >= length(b) + return false + end + sb = Set(b) + @inbounds for x in a + if x ∉ sb + return false + end + end + return true +end + +""" + induced_components(g::SimpleGraph, vs::AbstractVector{Int})::Vector{RegionKey} +Return connected components (as RegionKey) of the induced subgraph on `vs`. +""" +function induced_components(g::SimpleGraph, vs::AbstractVector{Int})::Vector{RegionKey} + if isempty(vs) + return RegionKey[] + end + h, vmap = induced_subgraph(g, vs) # h has vertices 1..nv(h), vmap maps h->g + comps = connected_components(h) # comps as vectors of 1..nv(h) + # map back to original vertex IDs via `vs` + return [to_key(vs[c]) for c in comps] +end + +# --- Maximal regions under inclusion ------------------------------------------ + +""" + maximal_regions(regions::Set{RegionKey})::Set{RegionKey} +Select the inclusion-maximal regions from a set. +""" +function maximal_regions(regions::Set{RegionKey})::Set{RegionKey} + keys = collect(regions) + sort!(keys; by=length) # small to large + maximal = Set{RegionKey}() + for i in eachindex(keys) + a = keys[i] + is_sub = false + for j in (i+1):length(keys) + b = keys[j] + if is_proper_subset(a, b) + is_sub = true + break + end + end + if !is_sub + push!(maximal, a) + end + end + return maximal +end + +# --- Close under intersections (with connected components) --------------------- +# Note that keeping the connected components of graphs with >1 component is unnecessary and is included here as legacy. +# use close_under_intersections_connected instead +""" + close_under_intersections(g::SimpleGraph, seed::Set{RegionKey}; loop_only::Bool=true)::Set{RegionKey} +Given `seed` regions, iteratively add intersections (split into connected components), +optionally keeping only loopful components; stop when no new regions appear. Optionally, only keep if component contains `must_contain` vertices. +""" +function close_under_intersections(g::SimpleGraph, seed::Set{RegionKey}; loop_only::Bool=true, must_contain = []) + R = Set(seed) + changed = true + while changed + changed = false + keys = collect(R) + for i in 1:length(keys)-1 + a = keys[i] + for j in (i+1):length(keys) + b = keys[j] + X = key_intersection(a, b) + comps = induced_components(g, X) + + for comp in comps + if loop_only && !is_loopful(g, comp) + continue + end + if !isempty(must_contain) && intersect(must_contain, comp) != must_contain + continue + end + if comp ∉ R + push!(R, comp) + changed = true + end + end + end + end + end + return R +end + +# --- Close under intersections (with connected components) --------------------- + +""" + close_under_intersections(g::SimpleGraph, seed::Set{RegionKey}; loop_only::Bool=true)::Set{RegionKey} +Given `seed` regions, iteratively add intersections. Only keep connected components. +optionally keeping only loopful components; stop when no new regions appear. Optionally, only keep if component contains `must_contain` vertices. +""" +function close_under_intersections_connected(g::SimpleGraph, seed::Set{RegionKey}; loop_only::Bool=true, must_contain = []) + R = Set(seed) + changed = true + while changed + changed = false + keys = collect(R) + for i in 1:length(keys)-1 + a = keys[i] + for j in (i+1):length(keys) + b = keys[j] + X = key_intersection(a, b) + + # must be connected and nonempty + if isempty(X) || !is_connected(induced_subgraph(g, X)[1]) + continue + end + + # must contain the required vertices + if !isempty(must_contain) && intersect(must_contain, X) != must_contain + continue + end + + comp = to_key(X) + # must be loopy, if loop_only + if loop_only && !is_loopful(g, comp) + continue + end + + if comp ∉ R + push!(R, comp) + changed = true + end + end + end + end + return R +end + +# --- Counting numbers (top-down Möbius) --------------------------------------- + +""" + counting_numbers(regions::Set{RegionKey})::Dict{RegionKey,Int} +Compute inclusion–exclusion counting numbers c(r): set c=1 for maximals, +then for other regions c(r) = 1 - sum_{a ⊃ r} c(a). +""" +function counting_numbers(regions::Set{RegionKey},maximals::Set{RegionKey})::Dict{RegionKey,Int} + R = collect(regions) + # sort supersets first (decreasing size) + sort!(R; by=r -> (-length(r), r)) + c = Dict{RegionKey,Int}() + + for r=maximals + c[r] = 1 + end + + # fill others (now supersets are guaranteed to have c set already) + for r in R + if haskey(c, r) + continue + end + s = 0 + for a in R + if is_proper_subset(r, a) + s += get(c, a, 0) + end + end + c[r] = 1 - s + end + return c +end + + +""" +Find all subgraphs of g that contain both u and v, up to C vertices, +and have no other leaves +""" +function vertex_walks_up_to_C_regions(g::AbstractGraph, u::Integer, v::Integer, C::Int; buffer::Int = C) + walks = Set{RegionKey}() + + stack = [(u, Set(), -1,false,0)] # (current vertex, path, previous vertex, has_both,num_steps) + + while !isempty(stack) + node, path, prev_node, has_both,num_steps = pop!(stack) + if node==v + has_both = true + + end + + # contains both u and v, and completed a cycle, or a path from u to v + if has_both && (node==v || node in path) + push!(walks, to_key([vv for vv=union(Set([node]), path)])) + end + + push!(path, node) + @assert length(path) <= C + if num_steps==C + buffer + continue + end + for w in neighbors(g, node) + if w==prev_node + # don't backtrack + continue + end + + if length(path) < C || (w in path) + push!(stack, (w, copy(path), node, has_both, num_steps + 1)) + end + end + + end + return walks +end + +""" +Build clusters on graph g out of the regions regs +""" +function build_clusters(g::SimpleGraph, regs::Set; loop_only::Bool=true, must_contain = [], smart::Bool=true, verbose::Bool=false) + verbose && println("Finding maximal"); flush(stdout) + @time Rmax = maximal_regions(regs) + verbose && println("Finding intersections of $(length(Rmax)) regions"); flush(stdout) + if smart + R = close_under_intersections_connected(g, Rmax;loop_only=loop_only, must_contain = unique(must_contain)) + else + R = close_under_intersections(g, Rmax;loop_only=loop_only, must_contain = unique(must_contain)) + end + verbose && println("Finding counting numbers"); flush(stdout) + @time c = counting_numbers(R, Rmax) + return R, Rmax, c +end + +""" +Maps regions to NamedGraphs +""" +function map_regions_named(R::Set, Rmax::Set, c::Dict, vs_dict::Dictionary) + # map back to names + R = [map(v -> vs_dict[v], set) for set in R] + Rmax = [map(v -> vs_dict[v], set) for set in Rmax] + c = Dict(map(v -> vs_dict[v], key) => val for (key, val) in c) + return R, Rmax, c +end + +# prune branches except those ending at keep_vertices +function prune_cc(g::AbstractGraph, regions::Vector, counting_nums::Dict; keep_vertices = []) + counting_graphs = Dict() + for r=regions + if counting_nums[r] != 0 + eg = induced_subgraph(g, r)[1] + pb = Tuple(sort(prune_branches(eg, keep_vertices))) + # pg = induced_subgraph(eg, prune_branches(eg, keep_vertices))[1] + if haskey(counting_graphs, pb) + counting_graphs[pb] += counting_nums[r] + else + counting_graphs[pb] = counting_nums[r] + end + end + end + counting_graphs +end + +""" + build_region_family_correlation(g::SimpleGraph, u::Int, v::Int, C::Int) +Return (R, Rmax, c) where: + R :: Set{RegionKey} — full region family closed under intersections + Rmax :: Set{RegionKey} — maximal regions used as seeds + c :: Dict{RegionKey,Int} — counting numbers for all r ∈ R +""" +function build_region_family_correlation(g::SimpleGraph, u::Int, v::Int, C::Int; buffer::Int=C, smart::Bool=true, verbose::Bool=false) + verbose && println("Finding graphs"); flush(stdout) + @time regs = vertex_walks_up_to_C_regions(g,u,v,C; buffer = buffer) + R,Rmax,c = build_clusters(g, regs; loop_only = false, must_contain = [u,v], smart=smart, verbose = verbose) +end + +""" + Function to enumerate all generalized regions on a NamedGraph up to size Cluster_size, containing u and v. All other vertices have degree geq 2. + Counting numbers are found via top-down Möbius inversion. + For one-point function, just set u=v. + Returns (R, Rmax, c) where: + R :: Vector{Vector{T}} — full region family closed under intersections + Rmax :: Vector{Vector{T}} — maximal regions (largest generalizd loops) used as seeds and not subsets of any other regions + c :: Dict{Vector{T},Int} — counting numbers for all r ∈ R +""" +function build_region_family_correlation(ng::NamedGraph, u, v, Cluster_size::Int; buffer::Int = Cluster_size, smart::Bool=true, prune::Bool=true,verbose::Bool=false) + g, vs_dict = position_graph(ng), Dictionary([i for i in 1:nv(ng)], collect(vertices(ng))) + mapped_u, mapped_v = ng.vertices.index_positions[u], ng.vertices.index_positions[v] + R, Rmax, c = build_region_family_correlation(g, mapped_u, mapped_v, Cluster_size; buffer = buffer, smart=smart,verbose=verbose) + R, Rmax, c = map_regions_named(R, Rmax, c, vs_dict) + if prune + c = prune_cc(ng,R,c; keep_vertices=[u,v]) + return collect(keys(c)), Rmax, c + else + return R, Rmax, c + end +end + +""" +Generate graphs up to isomorphism and then embed in the larger graph g. +max_v is max number of vertices +min_v is min number of vertices (e.g. 4 on square lattice) +""" +function generate_embedded_leafless_graphs(g::AbstractGraph, max_v::Int; min_v::Int=4, triangle_free::Bool = true, min_deg::Int = 2) + @assert min_deg >= 2 + k = maximum([degree(g,v) for v=vertices(g)]) + mygraphs = [generate_graphs(no_vertices, k*no_vertices; triangle_free = triangle_free, min_deg = min_deg, max_deg = k, connected = true) + for no_vertices=min_v:max_v] + + # now embed each one + embeddings = [[embed_graphs(g, subg) for subg = mygs] for mygs=mygraphs] + subgraphs = vcat(vcat(embeddings...)...) + + # make into loopful induced subgraphs only + # Since they're connected and min_deg >= 2, will by definition be loopful + + regions = unique([to_key(unique(vcat([[e[1], e[2]] for e=subg]...))) for subg = subgraphs]) + return Set{RegionKey}(regions) +end + +""" + build_region_family(g::SimpleGraph, C::Int) +Return (R, Rmax, c) where: + R :: Set{RegionKey} — full region family closed under intersections + Rmax :: Set{RegionKey} — maximal regions used as seeds + c :: Dict{RegionKey,Int} — counting numbers for all r ∈ R +""" +function build_region_family(g::SimpleGraph, C::Int; min_deg::Int=2, min_v::Int=4,triangle_free::Bool=true, smart::Bool=true, verbose::Bool=false) + verbose && println("Finding graphs") + @time regs = generate_embedded_leafless_graphs(g, C; min_deg = min_deg, min_v=min_v,triangle_free=triangle_free) + build_clusters(g, regs; loop_only = true,smart=smart, verbose=verbose) +end + +""" + Function to enumerate all generalized regions on a NamedGraph up to size Cluster_size. + Counting numbers are found via top-down Möbius inversion. + Returns (R, Rmax, c) where: + R :: Vector{Vector{T}} — full region family closed under intersections + Rmax :: Vector{Vector{T}} — maximal regions (largest generalizd loops) used as seeds and not subsets of any other regions + c :: Dict{Vector{T},Int} — counting numbers for all r ∈ R +""" +function build_region_family(ng::NamedGraph, Cluster_size::Int; min_deg::Int=2, min_v::Int=4,triangle_free::Bool=true, smart::Bool=true, verbose::Bool=false) + g, vs_dict = position_graph(ng), Dictionary([i for i in 1:nv(ng)], collect(vertices(ng))) + R, Rmax, c = build_region_family(g, Cluster_size; min_deg = min_deg, min_v=min_v,triangle_free=triangle_free, smart=smart, verbose=verbose) + map_regions_named(R, Rmax, c, vs_dict) +end \ No newline at end of file diff --git a/src/TensorNetworkQuantumSimulator.jl b/src/TensorNetworkQuantumSimulator.jl index eb59923..9a29fb1 100644 --- a/src/TensorNetworkQuantumSimulator.jl +++ b/src/TensorNetworkQuantumSimulator.jl @@ -15,7 +15,10 @@ include("MessagePassing/abstractbeliefpropagationcache.jl") include("MessagePassing/beliefpropagationcache.jl") include("MessagePassing/boundarympscache.jl") include("MessagePassing/loopcorrection.jl") +include("MessagePassing/clustercorrections.jl") +include("MessagePassing/cumulant-clustercorrections.jl") include("graph_ops.jl") +include("graph_enumeration.jl") include("utils.jl") include("Apply/apply_gates.jl") diff --git a/src/graph_enumeration.jl b/src/graph_enumeration.jl new file mode 100644 index 0000000..a8cc524 --- /dev/null +++ b/src/graph_enumeration.jl @@ -0,0 +1,108 @@ +using Graphs: loadgraphs +using GraphIO +using Graphs.Experimental +using Graphs.SimpleGraphs +using NamedGraphs # from TensorNetworkQuantumSimulator +using StatsBase +using ProgressMeter + +# Generate non-isomorphic graphs with geng +# on square lattice, choose triangle_free = true to speed things up +function generate_graphs(n::Int, max_e::Int; min_e::Int=0, triangle_free::Bool = true, min_deg::Int=1, max_deg::Int=4, connected::Bool = true) + # Build geng command + if connected + if triangle_free + cmd = `geng -c -t -d$(min_deg) -D$(max_deg) $n $(min_e):$(max_e)` + else + cmd = `geng -c -d$(min_deg) -D$(max_deg) $n $(min_e):$(max_e)` + end + elseif triangle_free + cmd = `geng -t -d$(min_deg) -D$(max_deg) $n $(min_e):$(max_e)` + else + cmd = `geng -d$(min_deg) -D$(max_deg) $n $(min_e):$(max_e)` + end + + graphs = SimpleGraph[] + try open(cmd, "r") do io + for g in loadgraphs(io, GraphIO.Graph6.Graph6Format()) + push!(graphs, g[2]) + end + end + catch e + println("Couldn't find graphs with these parameters: $(cmd)") + end + + # also + return graphs +end + +function map_edges(subg::SimpleGraph, iso_map::Vector) + iso_map_dict = Dict(v[2]=>v[1] for v=iso_map) + [(min(iso_map_dict[src(e)],iso_map_dict[dst(e)]),max(iso_map_dict[src(e)],iso_map_dict[dst(e)])) for e=edges(subg)] +end + +function embed_graphs(g::AbstractGraph, subg::SimpleGraph) + all_embeddings = unique(sort.([map_edges(subg, iso_map) for iso_map=all_subgraphisomorph(g, subg)])) +end + +function is_valid_graph(g; leaf_vertices = []) + vertex_counts = countmap(vcat([[e[1],e[2]] for e=g]...)) + for (v,k)=vertex_counts + if !(v in leaf_vertices) && k < 2 + return false + end + end + return true +end + +""" +Generate graphs up to isomorphism and then embed in the larger graph g. +max_weight is max number of edges +min_v is min number of vertices (e.g. 4 on square lattice) +""" +function generate_embedded_graphs(g::AbstractGraph, max_weight::Int; min_v::Int=4, triangle_free::Bool = true, min_deg::Int = 2, leaf_vertices = []) + + k = maximum([degree(g,v) for v=vertices(g)]) + # max edge weight in max_weight, so the max number of vertices is max_weight-1 + mygraphs = [generate_graphs(no_vertices, max_weight; min_e=no_vertices-1, triangle_free = triangle_free, min_deg = min_deg, max_deg = k, connected = true) for no_vertices=min_v:max_weight+1] + + # only try to embed graphs with at most length(leaf_vertices) leaves + mygraphs = [filter(subg->count(isone, [degree(subg,v) for v=vertices(subg)])<=length(leaf_vertices), mygs) for mygs=mygraphs] + + # now embed each one + embeddings = [[embed_graphs(g, subg) for subg = mygs] for mygs=mygraphs] + subgraphs = filter(g->is_valid_graph(g; leaf_vertices = leaf_vertices), vcat(vcat(embeddings...)...)) +end + +""" + prune_branches(g::AbstractGraph, keep_vertices) + +Return a new graph obtained from `g` by pruning away leaf branches +(iteratively) except those that terminate at any vertex in `keep_vertices`. + +Returns the vertices in the pruned graph +""" +function prune_branches(g::AbstractGraph, keep_vertices) + keep_set = Set(keep_vertices) + alive = Dict(v=>true for v=vertices(g)) + + changed = true + while changed + changed = false + # compute degree within the induced alive-subgraph (count alive neighbors) + for v=vertices(g) + if alive[v] + cnt = sum([alive[u] for u=neighbors(g,v)]) + + # remove if it's a leaf (degree 1) or isolated (degree 0), + # and not in keep_set. + if cnt <= 1 && !(v in keep_set) + alive[v] = false + changed = true + end + end + end + end + + return [v for v=vertices(g) if alive[v]] +end From 63549486c7aed8fa97516f2613e353199e341c07 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sat, 15 Nov 2025 23:50:07 -0500 Subject: [PATCH 26/36] bp message diffs --- PEPO_2d_ising.jl | 6 +++--- src/MessagePassing/abstractbeliefpropagationcache.jl | 5 +++++ src/MessagePassing/beliefpropagationcache.jl | 5 ++++- 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl index 4fd8ae3..3a61f65 100644 --- a/PEPO_2d_ising.jl +++ b/PEPO_2d_ising.jl @@ -109,7 +109,7 @@ function expect_bmps(dat::Dict; obs = "X", MPS_message_rank::Int = 10, save_tag expect_vals = zeros(length(all_verts), length(dat["sqrtρ"])-start_i+1) for i=start_i:length(dat["sqrtρ"]) if use_gpu - sqrtρ = CUDA.cu(dat["sqrtρ"][i]) + sqrtρ = adapt(CUDA.CuArray{ComplexF64})(dat["sqrtρ"][i]) else sqrtρ = dat["sqrtρ"][i] end @@ -178,10 +178,10 @@ function evolve_bp(ρρ::BeliefPropagationCache, n::Int, nsteps::Int; β = 0, hx apply_kwargs = (; maxdim = χ, cutoff = 1e-12) bp_update_kwargs = (; maxiter = 50, tolerance = 1e-8, verbose = true) if use_gpu - ρρ = CUDA.cu(ρρ) + ρρ = adapt(CUDA.CuArray{ComplexF64})(ρρ) end two_qubit_gates = [adapt(datatype(network(ρρ)), TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s)) for pair=vcat(ec...)] - + single_qubit_gates = [adapt(datatype(network(ρρ)), TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s)) for v=vertices(g)] layer = vcat(single_qubit_gates, two_qubit_gates, single_qubit_gates) diff --git a/src/MessagePassing/abstractbeliefpropagationcache.jl b/src/MessagePassing/abstractbeliefpropagationcache.jl index 46a593e..b9ff541 100644 --- a/src/MessagePassing/abstractbeliefpropagationcache.jl +++ b/src/MessagePassing/abstractbeliefpropagationcache.jl @@ -199,6 +199,7 @@ function update(alg::Algorithm"bp", bpc::AbstractBeliefPropagationCache) bpc = copy(bpc) diffs = zeros(alg.kwargs.maxiter) tot_iter = alg.kwargs.maxiter + success = false for i in 1:alg.kwargs.maxiter diff = compute_error ? Ref(0.0) : nothing update_iteration!(alg, bpc, alg.kwargs.edge_sequence; (update_diff!) = diff) @@ -208,12 +209,16 @@ function update(alg::Algorithm"bp", bpc::AbstractBeliefPropagationCache) if alg.kwargs.verbose println("BP converged to desired precision after $i iterations.") end + success = true tot_iter = i break end end end if compute_error && alg.kwargs.verbose + if !success + println("Did not converge.") + end println("Diffs during message passing: $(diffs[1:tot_iter])") end return bpc diff --git a/src/MessagePassing/beliefpropagationcache.jl b/src/MessagePassing/beliefpropagationcache.jl index 7da4c33..ba81ec4 100644 --- a/src/MessagePassing/beliefpropagationcache.jl +++ b/src/MessagePassing/beliefpropagationcache.jl @@ -13,10 +13,13 @@ struct BeliefPropagationCache{V, N <: AbstractTensorNetwork{V}, M <: Union{ITens end #TODO: Take `dot` without precontracting the messages to allow scaling to more complex messages +# we shouldn't let this be negative, it doesn't make sense function message_diff(message_a::ITensor, message_b::ITensor) n_a, n_b = norm(message_a), norm(message_b) f = abs2(dot(message_a, message_b) / (n_a * n_b)) - return 1 - f + + # or do abs(1-f)? + return max(0,1 - f) end messages(bp_cache::BeliefPropagationCache) = bp_cache.messages From ad05ebe4a55058fe7e8e2fdf3ba85cff5da6f784 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sun, 16 Nov 2025 10:59:08 -0500 Subject: [PATCH 27/36] added comments --- cluster/expect-corrected.jl | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/cluster/expect-corrected.jl b/cluster/expect-corrected.jl index beb28fb..c301bc5 100644 --- a/cluster/expect-corrected.jl +++ b/cluster/expect-corrected.jl @@ -42,13 +42,14 @@ function prep_insertions(obs) return (coeffs = hyper_coeff, op_strings = insertion_operator) end -# Cluster/cluster cumulant weights +""" +Cluster expansion. See clustercorrections.jl +""" function cluster_weights(bpc::BeliefPropagationCache, clusters::Vector, egs::Vector{<:AbstractNamedGraph}, interaction_graph; obs = nothing) kwargs = prep_insertions(obs) logZbp = TN.free_energy(bpc; kwargs...) - println("MADE IT HERE"); flush(stdout) isempty(egs) && return [0], [[logZbp]], [[1]] circuit_lengths = sort(unique([c.weight for c=clusters])) @@ -65,7 +66,7 @@ function cluster_weights(bpc::BeliefPropagationCache, clusters::Vector, egs::Vec coeffs = Array{Array}(undef, length(circuit_lengths) + 1) coeffs[1] = [1] - println("MADE IT HERE TOO"); flush(stdout) + # now calculate contribution to logZ from each cluster for (cl_i, cl)=enumerate(circuit_lengths) clusters_cl = filter(c->c.weight==cl, clusters) @@ -75,7 +76,10 @@ function cluster_weights(bpc::BeliefPropagationCache, clusters::Vector, egs::Vec return vcat([0],circuit_lengths), logZs, coeffs end - + +""" +Cluster cumulant expansion. See cumulant-clustercorrections.jl +""" function cc_weights(bpc::BeliefPropagationCache, regions::Vector, counting_nums::Dict; obs = nothing, rescale::Bool = false) kwargs = prep_insertions(obs) @@ -99,7 +103,9 @@ function cc_weights(bpc::BeliefPropagationCache, regions::Vector, counting_nums: return log.(wts), [counting_nums[gg] for gg=regions[use_g]] end -# onepoint or twopoint connected correlation function, using cluster cumulant expansion +""" +onepoint or twopoint connected correlation function, using cluster cumulant expansion +""" function cc_correlation(bpc::BeliefPropagationCache, regions::Vector, counting_nums::Dict, obs) logZs, cnums = cc_weights(bpc, regions, counting_nums; obs = obs) op_strings, verts, _ = TN.collectobservable(obs) @@ -110,7 +116,9 @@ function cc_correlation(bpc::BeliefPropagationCache, regions::Vector, counting_n end end -# onepoint or twopoint connected correlation function, using cluster expansion +""" +onepoint or twopoint connected correlation function, using cluster expansion +""" function cluster_correlation(bpc::BeliefPropagationCache, clusters::Vector, egs::Vector, interaction_graph, obs) cluster_wts, logZs, ursells = cluster_weights(bpc, clusters, egs, interaction_graph; obs = obs) op_strings, verts, _ = TN.collectobservable(obs) From f6c46548dc0828c428bf69d93249833fd2424f45 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sun, 16 Nov 2025 11:02:10 -0500 Subject: [PATCH 28/36] removed finite temp evolution --- Finitetemp_2d_ising_fullupdate.jl | 127 ------------------- PEPO_2d_ising.jl | 197 ------------------------------ PEPO_2d_ising_cpu.jl | 117 ------------------ 3 files changed, 441 deletions(-) delete mode 100644 Finitetemp_2d_ising_fullupdate.jl delete mode 100644 PEPO_2d_ising.jl delete mode 100644 PEPO_2d_ising_cpu.jl diff --git a/Finitetemp_2d_ising_fullupdate.jl b/Finitetemp_2d_ising_fullupdate.jl deleted file mode 100644 index 7577c5e..0000000 --- a/Finitetemp_2d_ising_fullupdate.jl +++ /dev/null @@ -1,127 +0,0 @@ -using TensorNetworkQuantumSimulator -const TN = TensorNetworkQuantumSimulator - -using ITensors: @OpName_str, @SiteType_str, Algorithm, datatype, ITensors - -using NamedGraphs: NamedGraphs, edges, NamedEdge -using Graphs -const NG = NamedGraphs -const G = Graphs -using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph -using NamedGraphs.GraphsExtensions: add_edges, add_vertices - -using Random -using TOML - -using Base.Threads -using MKL -using LinearAlgebra -using NPZ - -using CUDA - -using Adapt -using Dictionaries - -BLAS.set_num_threads(min(6, Sys.CPU_THREADS)) -println("Julia is using "*string(nthreads())) -println("BLAS is using "*string(BLAS.get_num_threads())) - -#Gate : rho -> rho .X. With this defined, expect(sqrt_rho, X; alg) = Tr(sqrt_rho . X sqrt_rho) / Tr(sqrt_rho sqrt_rho) = Tr(rho . X) / Tr(rho) -function ITensors.op( - ::OpName"X", ::SiteType"Pauli" - ) - mat = zeros(ComplexF64, 4,4) - mat[1,2] = 1 - mat[2,1] = 1 - mat[3,4] = im - mat[4,3] = -im - return mat -end - -function main() - - n = 6 - g = named_grid((n,n)) - #Pauli inds run over identity, X, Y, Z - s = siteinds("Pauli", g) - ρ = identitytensornetworkstate(ComplexF64, g, s) - ITensors.disable_warn_order() - - use_gpu = false - - δβ = 0.01 - hx = -3.1 - J = -1 - - # #Do a custom 4-way edge coloring then Trotterise the Hamiltonian into commuting groups - ec1 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 1:2:(n-1)] for i in 1:n]) - ec2 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 2:2:(n-1)] for i in 1:n]) - ec3 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 1:2:(n-1)] for i in 1:n]) - ec4 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 2:2:(n-1)] for i in 1:n]) - ec = [ec1, ec2, ec3, ec4] - - @assert length(reduce(vcat, ec)) == length(edges(g)) - nsteps = 50 - apply_kwargs = (; maxdim = 4, cutoff = 1e-12) - - MPS_message_rank = 10 - - β = 0 - for i in 1:nsteps - #Apply the singsite rotations half way - for v in vertices(g) - gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) - gate = adapt(datatype(ρ), gate) - TN.setindex_preserve!(ρ, normalize(ITensors.apply(gate, ρ[v])), v) - end - - #Apply the two site rotations, use a boundary MPS cache to apply them (need to run column or row wise depending on the gates) - for (k, colored_edges) in enumerate(ec) - - #Only if you want to use GPU to do boundary MPS - if use_gpu - ρ_gpu =CUDA.cu(ρ) - ρρ = TN.BoundaryMPSCache(ρ_gpu, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) - else - ρρ = TN.BoundaryMPSCache(ρ, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) - end - ρρ = TN.update(ρρ) - TN.update_partitions!(ρρ, collect(TN.partitionvertices(TN.supergraph(ρρ)))) - - for pair in colored_edges - gate = TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s) - gate = adapt(datatype(ρ), gate) - envs = TN.incoming_messages(ρρ, [src(pair), dst(pair)]) - envs = adapt(datatype(ρ)).(envs) - ρv1, ρv2 = TN.full_update(gate, ρ, [src(pair), dst(pair)]; envs, print_fidelity_loss = true, apply_kwargs...) - TN.setindex_preserve!(ρ, normalize(ρv1), src(pair)) - TN.setindex_preserve!(ρ, normalize(ρv2), dst(pair)) - end - end - - - for v in vertices(g) - gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) - gate = adapt(datatype(ρ), gate) - TN.setindex_preserve!(ρ, normalize(ITensors.apply(gate, ρ[v])), v) - end - - β += δβ - - if use_gpu - sz_doubled = TN.expect(ρ, ("X", [(2,2)]); alg = "boundarymps", mps_bond_dimension = MPS_message_rank) - else - sz_doubled = TN.expect(CUDA.cu(ρ), ("X", [(2,2)]); alg = "boundarymps", mps_bond_dimension = MPS_message_rank) - end - - println("Inverse Temperature is $β") - println("Bond dimension of PEPO $(TN.maxvirtualdim(ρ))") - - println("Expectation value at beta = $(2*β) is $(sz_doubled)") - end - - -end - -main() \ No newline at end of file diff --git a/PEPO_2d_ising.jl b/PEPO_2d_ising.jl deleted file mode 100644 index 3a61f65..0000000 --- a/PEPO_2d_ising.jl +++ /dev/null @@ -1,197 +0,0 @@ -using TensorNetworkQuantumSimulator -const TN = TensorNetworkQuantumSimulator - -using ITensors: @OpName_str, @SiteType_str, Algorithm, datatype, ITensors - -using NamedGraphs: NamedGraphs, edges, NamedEdge -using Graphs -const NG = NamedGraphs -const G = Graphs -using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph -using NamedGraphs.GraphsExtensions: add_edges, add_vertices - -using Random -using TOML - -using Base.Threads -using MKL -using LinearAlgebra -using NPZ - -using CUDA - -using Adapt -using Dictionaries -using JLD2 - -BLAS.set_num_threads(min(6, Sys.CPU_THREADS)) -println("Julia is using "*string(nthreads())) -println("BLAS is using "*string(BLAS.get_num_threads())) - -DATA_DIR = "/mnt/ceph/users/gsommers/data/" - -#Gate : rho -> rho .X -function ITensors.op( - ::OpName"X", ::SiteType"Pauli" - ) - mat = zeros(ComplexF64, 4,4) - mat[1,2] = 1 - mat[2,1] = 1 - mat[3,4] = im - mat[4,3] = -im - return mat -end - -function prep_edges(n::Int, g::NamedGraphs.AbstractNamedGraph) - # #Do a custom 4-way edge coloring then Trotterize the Hamiltonian into commuting groups - ec1 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 1:2:(n-1)] for i in 1:n]) - ec2 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 2:2:(n-1)] for i in 1:n]) - ec3 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 1:2:(n-1)] for i in 1:n]) - ec4 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 2:2:(n-1)] for i in 1:n]) - ec = [ec1, ec2, ec3, ec4] - - @assert length(reduce(vcat, ec)) == length(edges(g)) - ec -end - -# apply layer of single qubit gates -function apply_single_qubit_layer!(ρ::TensorNetworkState, gates::Dict) - for v=keys(gates) - TN.setindex_preserve!(ρ, normalize(ITensors.apply(gates[v], ρ[v])), v) - end -end - -#Apply the two site rotations, use a boundary MPS cache to apply them (need to run column or row wise depending on the gates) -function apply_two_qubit_layer!(ρ::TensorNetworkState, ec::Array, gates::Dict; MPS_message_rank::Int, use_gpu::Bool=true, apply_kwargs...) - for (k, colored_edges) in enumerate(ec) - - #Only if you want to use GPU to do boundary MPS - println("Starting boundary MPS cache") - if use_gpu - @time ρ_gpu =CUDA.cu(ρ) - @time ρρ = TN.BoundaryMPSCache(ρ_gpu, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) - else - @time ρρ = TN.BoundaryMPSCache(ρ, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) - end - @time ρρ = TN.update(ρρ) - @time TN.update_partitions!(ρρ, collect(TN.partitionvertices(TN.supergraph(ρρ)))) - - println("Starting two-qubit gates") - @time begin - for pair in colored_edges - apply_two_qubit_gate!(ρ,ρρ, gates[pair], pair; apply_kwargs...) - end - end - end -end - -function apply_two_qubit_gate!(ρ::TensorNetworkState,ρρ::TN.BoundaryMPSCache, gate::ITensors.ITensor, pair::NamedEdge; apply_kwargs...) - envs = TN.incoming_messages(ρρ, [src(pair), dst(pair)]) - envs = adapt(datatype(ρ)).(envs) - ρv1, ρv2 = TN.full_update(gate, ρ, [src(pair), dst(pair)]; envs, print_fidelity_loss = true, apply_kwargs...) - TN.setindex_preserve!(ρ, normalize(ρv1), src(pair)) - TN.setindex_preserve!(ρ, normalize(ρv2), dst(pair)) -end - -function intermediate_save(sqrtρ, β; δβ::Float64, χ::Int, n::Int, MPS_message_rank, save_tag = "", hx = -3.04438) - dat = Dict("L"=>n, "δβ"=>δβ, "β"=>β, "χ"=>χ, "sqrtρ"=>sqrtρ, "mps_rank"=>MPS_message_rank, "hx"=>hx) - save(DATA_DIR * "$(save_tag)L$(n)_χ$(χ)_D$(MPS_message_rank)_step$(round(δβ,digits=3))_$(round(β,digits=3)).jld2", dat) -end - -function intermediate_save_bp(ρ, errs, β; δβ::Float64, χ::Int, n::Int, save_tag = "", hx = -3.04438) - dat = Dict("L"=>n, "δβ"=>δβ, "β"=>β, "χ"=>χ, "ρ"=>ρ, "errs"=>errs, "hx"=>hx) - dat["X"] = expect(ρ, [("X", [v]) for v=vertices(network(ρ).tensornetwork.graph)]) - save(DATA_DIR * "$(save_tag)L$(n)_χ$(χ)_step$(round(δβ,digits=3))_$(round(β,digits=3)).jld2", dat) -end - -function expect_bmps(dat::Dict; obs = "X", MPS_message_rank::Int = 10, save_tag = "", use_gpu::Bool = true, start_i::Int = 1) - all_verts = collect(vertices(dat["sqrtρ"][1].tensornetwork.graph)) - expect_vals = zeros(length(all_verts), length(dat["sqrtρ"])-start_i+1) - for i=start_i:length(dat["sqrtρ"]) - if use_gpu - sqrtρ = adapt(CUDA.CuArray{ComplexF64})(dat["sqrtρ"][i]) - else - sqrtρ = dat["sqrtρ"][i] - end - @time expect_vals[:,i-start_i+1] = real.(TN.expect(sqrtρ, [(obs, [v]) for v=all_verts]; alg = "boundarymps", mps_bond_dimension = MPS_message_rank)) - save(DATA_DIR * "$(save_tag)L$(dat["L"])_χ$(dat["χ"])_D$(MPS_message_rank)_step$(dat["δβ"])_$(dat["β"][i]).jld2", Dict(obs=>expect_vals[:,i-start_i+1], "verts"=>all_verts, "hx"=>dat["hx"], "β"=>dat["β"][i], "χ"=>[dat["χ"],maxvirtualdim(dat["sqrtρ"][i])], "mps_rank"=>MPS_message_rank, "δβ"=>dat["δβ"], "L"=>dat["L"])) - flush(stdout) - end - all_verts, expect_vals -end - -function evolve_bmps(n::Int, nsteps::Int; hx=-3.04438, δβ = 0.01, use_gpu::Bool = true, χ::Int=4, MPS_message_rank::Int = 10, save_tag = "") - - g = named_grid((n,n)) - s = siteinds("Pauli", g) - ρ = identitytensornetworkstate(ComplexF64, g, s) - evolve_bmps(ρ, n, nsteps; β=0, hx=hx, δβ=δβ, use_gpu = use_gpu, χ=χ, MPS_message_rank = MPS_message_rank, save_tag = save_tag) - -end - -function evolve_bmps(ρ::TensorNetworkState, n::Int, nsteps::Int; β = 0, hx=-3.04438, δβ = 0.01, use_gpu::Bool=true, χ::Int=4, MPS_message_rank::Int = 10, save_tag = "") - g = ρ.tensornetwork.graph - s = siteinds(ρ) - ITensors.disable_warn_order() - - J = -1 - - ec = prep_edges(n, g) - apply_kwargs = (; maxdim = χ, cutoff = 1e-12) - - two_qubit_gates = Dict(pair=>adapt(datatype(ρ), TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s)) for pair=vcat(ec...)) - - single_qubit_gates = Dict(v=>adapt(datatype(ρ), TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s)) for v=vertices(g)) - - for i in 1:nsteps - #Apply the singsite rotations half way - apply_single_qubit_layer!(ρ, single_qubit_gates) - - @time apply_two_qubit_layer!(ρ, ec, two_qubit_gates; MPS_message_rank = MPS_message_rank, use_gpu = use_gpu, apply_kwargs...) - - apply_single_qubit_layer!(ρ, single_qubit_gates) - - β += δβ - - println("Inverse Temperature is $(β)"); flush(stdout) - - intermediate_save(ρ,β; χ=χ,n=n,MPS_message_rank = MPS_message_rank, δβ=δβ, hx=hx, save_tag = save_tag) - end -end - -function evolve_bp(n::Int, nsteps::Int; hx=-3.04438, δβ = 0.01, use_gpu::Bool=true, χ::Int=4, save_tag = "") - g = named_grid((n,n)) - s = siteinds("Pauli", g) - ρ = identitytensornetworkstate(ComplexF64, g, s) - ρρ = BeliefPropagationCache(ρ) - evolve_bp(ρρ, n, nsteps; β=0, hx=hx, δβ=δβ, use_gpu = use_gpu, χ=χ, save_tag = save_tag) -end - -function evolve_bp(ρρ::BeliefPropagationCache, n::Int, nsteps::Int; β = 0, hx=-3.04438, δβ = 0.01, use_gpu::Bool=true, χ::Int=4, save_tag = "") - g = network(ρρ).tensornetwork.graph - s = siteinds(network(ρρ)) - ITensors.disable_warn_order() - - J = -1 - - ec = prep_edges(n, g) - apply_kwargs = (; maxdim = χ, cutoff = 1e-12) - bp_update_kwargs = (; maxiter = 50, tolerance = 1e-8, verbose = true) - if use_gpu - ρρ = adapt(CUDA.CuArray{ComplexF64})(ρρ) - end - two_qubit_gates = [adapt(datatype(network(ρρ)), TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s)) for pair=vcat(ec...)] - - single_qubit_gates = [adapt(datatype(network(ρρ)), TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s)) for v=vertices(g)] - - layer = vcat(single_qubit_gates, two_qubit_gates, single_qubit_gates) - for i in 1:nsteps - - @time ρρ, errs = apply_gates(layer, ρρ; bp_update_kwargs = bp_update_kwargs, apply_kwargs) - - β += δβ - - println("Inverse Temperature is $(β)"); flush(stdout) - intermediate_save_bp(ρρ,errs,β; χ=χ,n=n, δβ=δβ, hx=hx, save_tag = save_tag) - end -end diff --git a/PEPO_2d_ising_cpu.jl b/PEPO_2d_ising_cpu.jl deleted file mode 100644 index 5ea6e6a..0000000 --- a/PEPO_2d_ising_cpu.jl +++ /dev/null @@ -1,117 +0,0 @@ -using TensorNetworkQuantumSimulator -const TN = TensorNetworkQuantumSimulator - -using ITensorNetworks -const ITN = ITensorNetworks -using ITensors -using ITensors: @OpName_str, @SiteType_str, Algorithm, datatype - -using ITensorNetworks: AbstractBeliefPropagationCache, IndsNetwork, setindex_preserve_graph! -using NamedGraphs -using NamedGraphs: edges -using Graphs -const NG = NamedGraphs -const G = Graphs -using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph -using NamedGraphs.GraphsExtensions: add_edges, add_vertices - -using Random -using TOML - -using Base.Threads -using MKL -using LinearAlgebra -using NPZ - -using Adapt -using Dictionaries - -BLAS.set_num_threads(min(6, Sys.CPU_THREADS)) -println("Julia is using "*string(nthreads())) -println("BLAS is using "*string(BLAS.get_num_threads())) - -#Gate : rho -> rho .X -function ITensors.op( - ::OpName"X", ::SiteType"Pauli" - ) - mat = zeros(ComplexF64, 4,4) - mat[1,2] = 1 - mat[2,1] = 1 - mat[3,4] = im - mat[4,3] = -im - return mat -end - -function main() - - n = 6 - g = named_grid((n,n)) - s = siteinds("Pauli", g) - ρ = identitytensornetworkstate(ComplexF64, g, s) - ITensors.disable_warn_order() - - δβ = 0.01 - hx = -3.1 - J = -1 - - # #Do a custom 4-way edge coloring then Trotterise the Hamiltonian into commuting groups - ec1 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 1:2:(n-1)] for i in 1:n]) - ec2 = reduce(vcat, [[NamedEdge((j, i) => (j+1, i)) for j in 2:2:(n-1)] for i in 1:n]) - ec3 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 1:2:(n-1)] for i in 1:n]) - ec4 = reduce(vcat, [[NamedEdge((i,j) => (i, j+1)) for j in 2:2:(n-1)] for i in 1:n]) - ec = [ec1, ec2, ec3, ec4] - - @assert length(reduce(vcat, ec)) == length(edges(g)) - nsteps = 50 - apply_kwargs = (; maxdim = 4, cutoff = 1e-12) - - MPS_message_rank = 10 - - β = 0 - for i in 1:nsteps - #Apply the singsite rotations half way - for v in vertices(g) - gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) - gate = adapt(datatype(ρ), gate) - setindex_preserve_graph!(ρ, normalize(apply(gate, ρ[v])), v) - end - - #Apply the two site rotations, use a boundary MPS cache to apply them (need to run column or row wise depending on the gates) - for (k, colored_edges) in enumerate(ec) - - ρρ = TN.BoundaryMPSCache(ρ, MPS_message_rank; partition_by = (k== 1 || k == 2) ? "col" : "row", gauge_state = false) - - ρρ = TN.update(ρρ) - TN.update_partitions!(ρρ, collect(TN.partitionvertices(TN.supergraph(ρρ)))) - - for pair in colored_edges - gate = TN.toitensor(("Rzz", [src(pair), dst(pair)], -0.5*im * J * δβ), s) - gate = adapt(datatype(ρ), gate) - envs = TN.incoming_messages(ρρ, [src(pair), dst(pair)]) - envs = adapt(datatype(ρ)).(envs) - ρv1, ρv2 = ITensorNetworks.full_update_bp(gate, TN.tensornetwork(ρ), [src(pair), dst(pair)]; envs, apply_kwargs...) - ρ[src(pair)], ρ[dst(pair)] = normalize(ρv1), normalize(ρv2) - end - end - - - for v in vertices(g) - gate = TN.toitensor(("Rx", [v], -0.25 * im * hx *δβ), s) - gate = adapt(datatype(ρ), gate) - setindex_preserve_graph!(ρ, normalize(apply(gate, ρ[v])), v) - end - - β += δβ - - sz_doubled = TN.expect(ρ, ("X", [(2,2)]); alg = "boundarymps", mps_bond_dimension = MPS_message_rank) - - println("Inverse Temperature is $β") - println("Bond dimension of PEPO $(ITensorNetworks.maxlinkdim(ρ))") - - println("Expectation value at beta = $(2*β) is $(sz_doubled)") - end - - -end - -main() \ No newline at end of file From 36b51b888992d95475b5aa0bfd63074f35b05bbc Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sun, 16 Nov 2025 22:21:12 -0500 Subject: [PATCH 29/36] example with cluster expansions --- examples/clustercorrections.jl | 101 +++++++++++++++++++++++++++ src/TensorNetworkQuantumSimulator.jl | 1 - 2 files changed, 101 insertions(+), 1 deletion(-) create mode 100644 examples/clustercorrections.jl diff --git a/examples/clustercorrections.jl b/examples/clustercorrections.jl new file mode 100644 index 0000000..7c8e14a --- /dev/null +++ b/examples/clustercorrections.jl @@ -0,0 +1,101 @@ +using TensorNetworkQuantumSimulator +const TN = TensorNetworkQuantumSimulator + +using ITensors + +using NamedGraphs +using Graphs +const NG = NamedGraphs +const G = Graphs +using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph +using ProgressMeter + +using LinearAlgebra: norm + +using EinExprs: Greedy + +using Random +Random.seed!(1634) + +include("../cluster/expect-corrected.jl") + +function main(nx,ny) + χ = 3 + ITensors.disable_warn_order() + gs = [ + (named_grid((nx, 1)), "line", 0,-1), + (named_hexagonal_lattice_graph(nx, ny), "hexagonal", 6,11), + (named_grid((nx, ny)), "square", 4,11), + ] + + states = [] + for (g, g_str, smallest_loop_size, wmax) in gs + println("*****************************************") + println("Testing for $g_str lattice with $(NG.nv(g)) vertices") + wmax = min(wmax, NG.nv(g)) + ψ = TN.random_tensornetworkstate(ComplexF32, g, "S=1/2"; bond_dimension = χ) + + ψ = normalize(ψ; alg = "bp") + ψIψ = BeliefPropagationCache(ψ) + ψIψ = update(ψIψ) + + # BP expectation value + v = first(center(g)) + expect_bp = real(expect(ψIψ, ("Z", [v]))) + expect_exact_v = real(expect(ψ, ("Z", [v]); alg = "exact")) + clusters, egs, ig = TN.enumerate_clusters(g, wmax; must_contain=[v], min_deg = 1, min_v = smallest_loop_size) + cluster_wts, expects = cluster_correlation(ψIψ,clusters, egs, ig, ("Z", [v])) + + + regs = Dict() + cnums = Dict() + + cc_wts = [1; smallest_loop_size:wmax;] + for w=cc_wts + regs[w],_,cnums[w]=TN.build_region_family_correlation(g,v,v,w) + end + + expects_cc = Dict() + @showprogress for w=cc_wts + expects_cc[w] = real(cc_correlation(ψIψ,regs[w], cnums[w], ("Z", [v]))) + end + + println("Bp expectation value for Z on site $(v) is $expect_bp") + println("Cluster expansion expectation values: $(cluster_wts), $(real.(expects))") + println("Cluster cumulant expansion: $(cc_wts), $([expects_cc[w] for w=cc_wts])") + println("Exact expectation value is $expect_exact_v") + + println("***********************************") + u = neighbors(g, v)[1] + obs = (["Z","Z"], [u,v]) + expect_exact_u = real(expect(ψ, ("Z", [u]); alg = "exact")) + expect_exact = real(expect(ψ,obs; alg = "exact")) - expect_exact_u * expect_exact_v + println("Calculating connected correlation function between $(v) and $(u)") + + clusters, egs, ig = TN.enumerate_clusters(g, max(1,min(wmax,2*smallest_loop_size)); must_contain=[u,v], min_deg = 1, min_v = 2) + + cluster_wts, expects = cluster_correlation(ψIψ,clusters, egs, ig, obs) + + regs = Dict() + cnums = Dict() + + cc_wts = [2;3:wmax;] + for w=cc_wts + regs[w],_,cnums[w]=TN.build_region_family_correlation(g,u,v,w) + end + + expects_cc = Dict() + @showprogress for w=cc_wts + expects_cc[w] = real(cc_correlation(ψIψ,regs[w], cnums[w], obs)) + end + println("Cluster expansion expectation values: $(cluster_wts), $(real.(expects))") + println("Cluster cumulant expansion: $(cc_wts), $([expects_cc[w] for w=cc_wts])") + + println("Exact expectation value is $expect_exact") + push!(states, ψ) + end + return states +end + +#main(4,4) +#main(3,3) \ No newline at end of file diff --git a/src/TensorNetworkQuantumSimulator.jl b/src/TensorNetworkQuantumSimulator.jl index 9a29fb1..1373e56 100644 --- a/src/TensorNetworkQuantumSimulator.jl +++ b/src/TensorNetworkQuantumSimulator.jl @@ -18,7 +18,6 @@ include("MessagePassing/loopcorrection.jl") include("MessagePassing/clustercorrections.jl") include("MessagePassing/cumulant-clustercorrections.jl") include("graph_ops.jl") -include("graph_enumeration.jl") include("utils.jl") include("Apply/apply_gates.jl") From ce6c49e977d8568c178cce113152edad31221f1e Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Sun, 16 Nov 2025 22:46:38 -0500 Subject: [PATCH 30/36] removed extraneous packages from Project --- Project.toml | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/Project.toml b/Project.toml index 2a734d9..d043c23 100644 --- a/Project.toml +++ b/Project.toml @@ -10,30 +10,21 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a" -Debugger = "31a5f54b-26ea-5ae9-a837-f05ce5417438" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" EinExprs = "b1794770-133b-4de1-afb4-526377e9f4c5" -Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8" GraphIO = "aa1b3936-2fda-51b9-ab35-c553d3a640a2" GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97" -IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" -JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -JuliaSyntaxHighlighting = "ac6e5ff7-fb65-4e79-a425-ec3bc9c03011" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" -LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605" NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19" PauliPropagation = "293282d5-3c99-4fb6-92d0-fd3280a19750" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" -PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SimpleGraphAlgorithms = "41400c72-0c58-5c16-8579-4ecbce768449" SimpleGraphConverter = "205b04f2-f585-4877-a239-566270b3f673" @@ -42,37 +33,27 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" -UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] Adapt = "4.3.0" CUDA = "5.9.2" Combinatorics = "1.0.3" DataGraphs = "0.2.9" -Debugger = "0.7.15" Dictionaries = "0.4" EinExprs = "0.6.4" -Format = "1.3.7" GraphIO = "0.7.1" GraphRecipes = "0.5.13" Graphs = "1.8.0" HyperDualNumbers = "4.0.10" -IJulia = "1.32.1" ITensorMPS = "0.3.17" ITensors = "0.9" -JLD2 = "0.6.2" -JuliaSyntaxHighlighting = "0.1.0" KrylovKit = "0.10.2" -LaTeXStrings = "1.4.0" LinearAlgebra = "1.11.0" -LsqFit = "0.15.1" MKL = "0.9.0" NPZ = "0.4.3" NamedGraphs = "0.7" PauliPropagation = "0.3.0" -Plots = "1.41.1" ProgressMeter = "1.11.0" -PythonPlot = "1.0.6" Revise = "3.8.0" SimpleGraphAlgorithms = "0.6.0" SimpleGraphConverter = "0.1.0" @@ -81,7 +62,6 @@ Statistics = "1.11.1" StatsBase = "0.34.4" TensorOperations = "5.2" TypeParameterAccessors = "0.3.10" -UnicodePlots = "3.8.1" julia = "1.10" [extras] From 00370f5e120b0cad4566fa43fa51b5d64109f979 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Tue, 18 Nov 2025 20:48:46 -0500 Subject: [PATCH 31/36] cleaned up Project.toml --- Project.toml | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/Project.toml b/Project.toml index d043c23..0da6e77 100644 --- a/Project.toml +++ b/Project.toml @@ -7,9 +7,7 @@ version = "0.2.3" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" -DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" EinExprs = "b1794770-133b-4de1-afb4-526377e9f4c5" GraphIO = "aa1b3936-2fda-51b9-ab35-c553d3a640a2" @@ -20,11 +18,8 @@ ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" -NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605" NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19" PauliPropagation = "293282d5-3c99-4fb6-92d0-fd3280a19750" -ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SimpleGraphAlgorithms = "41400c72-0c58-5c16-8579-4ecbce768449" SimpleGraphConverter = "205b04f2-f585-4877-a239-566270b3f673" @@ -36,9 +31,7 @@ TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" [compat] Adapt = "4.3.0" -CUDA = "5.9.2" Combinatorics = "1.0.3" -DataGraphs = "0.2.9" Dictionaries = "0.4" EinExprs = "0.6.4" GraphIO = "0.7.1" @@ -49,11 +42,8 @@ ITensorMPS = "0.3.17" ITensors = "0.9" KrylovKit = "0.10.2" LinearAlgebra = "1.11.0" -MKL = "0.9.0" -NPZ = "0.4.3" NamedGraphs = "0.7" PauliPropagation = "0.3.0" -ProgressMeter = "1.11.0" Revise = "3.8.0" SimpleGraphAlgorithms = "0.6.0" SimpleGraphConverter = "0.1.0" From 1c2f735565618551c11aa4b204af19aaa11009ee Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Tue, 18 Nov 2025 20:54:46 -0500 Subject: [PATCH 32/36] changes from clusters branch --- src/MessagePassing/clustercorrections.jl | 25 +++---------------- .../cumulant-clustercorrections.jl | 24 ++---------------- 2 files changed, 5 insertions(+), 44 deletions(-) diff --git a/src/MessagePassing/clustercorrections.jl b/src/MessagePassing/clustercorrections.jl index e4ff304..63db8cb 100644 --- a/src/MessagePassing/clustercorrections.jl +++ b/src/MessagePassing/clustercorrections.jl @@ -1,9 +1,5 @@ - using NamedGraphs using NamedGraphs: AbstractNamedGraph -using ProgressMeter - -include("../graph_enumeration.jl") struct Cluster loop_ids::Vector{Int} @@ -13,8 +9,8 @@ struct Cluster end struct Loop - vertices::Vector{Int} - edges::Vector{Tuple{Int,Int}} + vertices::Vector + edges::Vector weight::Int end @@ -31,7 +27,6 @@ function build_interaction_graph(loops::Vector{Loop}) println(" Building optimized interaction graph for $n_loops loops...") flush(stdout) - progress = Progress(n_loops, dt=0.1, desc="Building graph: ", color=:green, barlen=50) # Optimization 1: Pre-compute vertex sets once vertex_sets = [Set(loop.vertices) for loop in loops] @@ -50,7 +45,6 @@ function build_interaction_graph(loops::Vector{Loop}) for i in 1:n_loops interaction_graph[i] = unique(vcat([vertex_to_loops[v] for v=loops[i].vertices]...)) - next!(progress) end return interaction_graph @@ -69,10 +63,6 @@ function dfs_enumerate_clusters_from_supported(all_loops::Vector{Loop}, supporte verbose && println(" Starting DFS cluster enumeration...") verbose && println(" Supported loops: $(length(supported_loop_ids)), Max weight: $max_weight") - # Progress tracking - last_report_time = time() - last_cluster_count = 0 - # DFS to grow clusters starting from each supported loop function dfs_grow_cluster(current_cluster::Vector{Int}, current_weight::Int, has_supported::Bool) @@ -97,16 +87,7 @@ function dfs_enumerate_clusters_from_supported(all_loops::Vector{Loop}, supporte if !(signature in seen_clusters) push!(seen_clusters, signature) push!(clusters, cluster) - cluster_count += 1 - - # Progress reporting every 2 seconds - current_time = time() - if current_time - last_report_time >= 2.0 - new_clusters = cluster_count - last_cluster_count - verbose && println(" Found $cluster_count clusters (+$new_clusters in last 2s)") - last_report_time = current_time - last_cluster_count = cluster_count - end + cluster_count += 1 end end diff --git a/src/MessagePassing/cumulant-clustercorrections.jl b/src/MessagePassing/cumulant-clustercorrections.jl index 7c1b3bd..8778d8e 100644 --- a/src/MessagePassing/cumulant-clustercorrections.jl +++ b/src/MessagePassing/cumulant-clustercorrections.jl @@ -1,11 +1,8 @@ using Graphs: nv, induced_subgraph, is_connected, connected_components using NamedGraphs using NamedGraphs: AbstractNamedGraph, position_graph -using ProgressMeter using Dictionaries -include("../graph_enumeration.jl") - const RegionKey = NTuple{N,Int} where {N} """ @@ -46,23 +43,6 @@ function is_loopful(g::SimpleGraph, key::RegionKey)::Bool return ne(h) - nv(h) + 1 > 0 end -""" - is_proper_subset(a::RegionKey, b::RegionKey)::Bool -Return true if region `a` is a proper subset of `b`. -""" -function is_proper_subset(a::RegionKey, b::RegionKey)::Bool - if length(a) >= length(b) - return false - end - sb = Set(b) - @inbounds for x in a - if x ∉ sb - return false - end - end - return true -end - """ induced_components(g::SimpleGraph, vs::AbstractVector{Int})::Vector{RegionKey} Return connected components (as RegionKey) of the induced subgraph on `vs`. @@ -92,7 +72,7 @@ function maximal_regions(regions::Set{RegionKey})::Set{RegionKey} is_sub = false for j in (i+1):length(keys) b = keys[j] - if is_proper_subset(a, b) + if a != b && issubset(a, b) is_sub = true break end @@ -212,7 +192,7 @@ function counting_numbers(regions::Set{RegionKey},maximals::Set{RegionKey})::Dic end s = 0 for a in R - if is_proper_subset(r, a) + if r!=a && issubset(r, a) s += get(c, a, 0) end end From 26a1ec26e3723adbf251cc5f7ca9e3f4e55f6990 Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Tue, 18 Nov 2025 20:58:35 -0500 Subject: [PATCH 33/36] merge from clusters --- src/MessagePassing/loopcorrection.jl | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl index 266b03c..9f36652 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/loopcorrection.jl @@ -64,16 +64,13 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) end #Get the all edges incident to the region specified by the vector of edges passed -# Optionally pass in a vector of vertices in the region, to handle the corner case where the region is just one vertex. function NamedGraphs.GraphsExtensions.boundary_edges( bpc::BeliefPropagationCache, - es::Vector{<:NamedEdge}; - vs::Vector = [] + es::Vector{<:NamedEdge} ) - if isempty(vs) - vs = unique(vcat(src.(es), dst.(es))) - end + vs = unique(vcat(src.(es), dst.(es))) + bpes = NamedEdge[] for v in vs incoming_es = boundary_edges(bpc, [v]; dir = :in) @@ -91,8 +88,10 @@ function weight(bpc::BeliefPropagationCache, eg; project_out::Bool = true, op_st if project_out bpc, antiprojectors = sim_edgeinduced_subgraph(bpc, eg) end - incoming_ms = - ITensor[message(bpc, e) for e in boundary_edges(bpc, es; vs = vs)] + if isempty(es) + incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, vs)] + else + incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, es)] local_tensors = reduce(vcat, bp_factors(bpc, vs; op_strings = op_strings, coeffs = coeffs, use_epsilon = true)) if project_out From 9fb5c49cfdedfbfd476306ab7d860c2fe0d5f45e Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Wed, 19 Nov 2025 10:06:43 -0500 Subject: [PATCH 34/36] cluster cumulant seems buggygit st! --- cluster/expect-corrected.jl | 7 +- cluster/time_evolution_tfim.jl | 83 ------------------- examples/clustercorrections.jl | 10 +-- src/Apply/full_update.jl | 7 +- .../abstractbeliefpropagationcache.jl | 5 +- src/MessagePassing/clustercorrections.jl | 2 +- .../cumulant-clustercorrections.jl | 27 +++++- src/MessagePassing/loopcorrection.jl | 1 + src/TensorNetworkQuantumSimulator.jl | 1 + src/graph_enumeration.jl | 1 - 10 files changed, 36 insertions(+), 108 deletions(-) delete mode 100644 cluster/time_evolution_tfim.jl diff --git a/cluster/expect-corrected.jl b/cluster/expect-corrected.jl index c301bc5..1d9a0ba 100644 --- a/cluster/expect-corrected.jl +++ b/cluster/expect-corrected.jl @@ -1,17 +1,14 @@ using TensorNetworkQuantumSimulator using NamedGraphs -using NamedGraphs: AbstractGraph, NamedGraph, AbstractNamedGraph -using NamedGraphs.GraphsExtensions: add_edge +using NamedGraphs: AbstractNamedGraph using Graphs const G = Graphs const NG = NamedGraphs const TN = TensorNetworkQuantumSimulator -using TensorNetworkQuantumSimulator: BoundaryMPSCache using HyperDualNumbers using Adapt: adapt using Dictionaries -using ITensors: inds, onehot, dag, commonind, op function prep_insertions(obs) if isnothing(obs) @@ -45,7 +42,7 @@ end """ Cluster expansion. See clustercorrections.jl """ -function cluster_weights(bpc::BeliefPropagationCache, clusters::Vector, egs::Vector{<:AbstractNamedGraph}, interaction_graph; obs = nothing) +function cluster_weights(bpc::BeliefPropagationCache, clusters::Vector, egs::Vector, interaction_graph; obs = nothing) kwargs = prep_insertions(obs) diff --git a/cluster/time_evolution_tfim.jl b/cluster/time_evolution_tfim.jl deleted file mode 100644 index 4175d0c..0000000 --- a/cluster/time_evolution_tfim.jl +++ /dev/null @@ -1,83 +0,0 @@ -using TensorNetworkQuantumSimulator -const TN = TensorNetworkQuantumSimulator - -using NamedGraphs.NamedGraphGenerators: named_grid -using Statistics -using Base.Threads -using ITensors -using CUDA - -function get_columns(L) - pairs = [] - for i=1:L, j=1:L - push!(pairs, [[(i,j),(i,k)] for k=j+1:L]...) - end - pairs -end - -function main(L, χ, bmps_ranks; nl::Int=20,θh = 0, use_gpu = true) - ITensors.disable_warn_order() - g = named_grid((L,L); periodic =false) - nq = length(vertices(g)) - - #Define the gate parameters - J = pi / 4 - - layer = [] - - Rx_layer = [("Rx", [v], θh) for v in vertices(g)] - ec = edge_color(g, 4) - - Rzz_layer = [] - for edge_group in ec - append!(Rzz_layer, ("Rzz", pair, -2*J) for pair in edge_group) - end - - layer = vcat(Rx_layer, Rzz_layer) - - pairs = get_columns(L) - verts = reshape([(i,j) for i=1:L,j=1:L],L^2) - - # the initial state (all up, use Float 32 precision) - ψ0 = tensornetworkstate(ComplexF32, v -> "↑", g, "S=1/2") - - # max bond dimension for the TN - apply_kwargs = (maxdim = χ, cutoff = 1.0e-12, normalize_tensors = false) - - # create the BP cache representing the square of the tensor network - ψ_bpc = BeliefPropagationCache(ψ0) - - # an array to keep track of expectations taken via two methods - bpc_states = Array{BeliefPropagationCache}(undef, nl) - errs = zeros(nl) - bmps_expects_z = [zeros(ComplexF64,L,L,nl) for r=bmps_ranks] - bmps_expects_zz = [zeros(ComplexF64,length(pairs),nl) for r=bmps_ranks] - - # evolve! (First step takes long due to compilation) - for l in 1:nl - println("Layer $l") - - t1 = @timed ψ_bpc, errors = - apply_gates(layer, ψ_bpc; apply_kwargs, verbose = false) - - bpc_states[l] = copy(ψ_bpc) - - #Boundary MPS expectation - if use_gpu - ψ = CUDA.cu(network(ψ_bpc)) - else - ψ = network(ψ_bpc) - end - - for r_i=1:length(bmps_ranks) - @time bmps_expects_zz[r_i][:,l] = expect(ψ, [(["Z","Z"], pair) for pair=pairs]; alg="boundarymps", mps_bond_dimension=bmps_ranks[r_i]) - @time z = expect(ψ, [("Z", [v]) for v=verts]; alg="boundarymps", mps_bond_dimension=bmps_ranks[r_i]) - bmps_expects_z[r_i][:,:,l] = reshape(z, (L,L)) - end - - println(" Took time: $(t1.time) [s]. Max bond dimension: $(maxvirtualdim(ψ_bpc))") - println(" Maximum Gate error for layer was $(maximum(errors))"); flush(stdout) - errs[l] = maximum(errors) - end - return bpc_states, bmps_expects_z, bmps_expects_zz, errs -end diff --git a/examples/clustercorrections.jl b/examples/clustercorrections.jl index 7c8e14a..406356b 100644 --- a/examples/clustercorrections.jl +++ b/examples/clustercorrections.jl @@ -8,7 +8,6 @@ using Graphs const NG = NamedGraphs const G = Graphs using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph -using ProgressMeter using LinearAlgebra: norm @@ -56,7 +55,7 @@ function main(nx,ny) end expects_cc = Dict() - @showprogress for w=cc_wts + for w=cc_wts expects_cc[w] = real(cc_correlation(ψIψ,regs[w], cnums[w], ("Z", [v]))) end @@ -85,7 +84,7 @@ function main(nx,ny) end expects_cc = Dict() - @showprogress for w=cc_wts + for w=cc_wts expects_cc[w] = real(cc_correlation(ψIψ,regs[w], cnums[w], obs)) end println("Cluster expansion expectation values: $(cluster_wts), $(real.(expects))") @@ -95,7 +94,4 @@ function main(nx,ny) push!(states, ψ) end return states -end - -#main(4,4) -#main(3,3) \ No newline at end of file +end \ No newline at end of file diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl index e5f331d..174deb0 100644 --- a/src/Apply/full_update.jl +++ b/src/Apply/full_update.jl @@ -107,6 +107,7 @@ function optimise_p_q( nfullupdatesweeps = 10, print_fidelity_loss = false, envisposdef = true, + verbose = false, apply_kwargs..., ) p_cur, q_cur = factorize( @@ -138,17 +139,17 @@ function optimise_p_q( b_vec = b(p, q, o, envs, q_cur) M_p_partial = partial(M_p, envs, q_cur, qs_ind) - p_cur, info = linsolve( + p_cur, info1 = linsolve( M_p_partial, b_vec, p_cur; isposdef = envisposdef, ishermitian = false ) b_tilde_vec = b(p, q, o, envs, p_cur) M_p_tilde_partial = partial(M_p, envs, p_cur, ps_ind) - q_cur, info = linsolve( + q_cur, info2 = linsolve( M_p_tilde_partial, b_tilde_vec, q_cur; isposdef = envisposdef, ishermitian = false ) - println(info) + verbose && println("Linsolve info, iteration $(i): $(info1), $(info2)") end fend = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0 diff --git a/src/MessagePassing/abstractbeliefpropagationcache.jl b/src/MessagePassing/abstractbeliefpropagationcache.jl index ef174af..bd70e05 100644 --- a/src/MessagePassing/abstractbeliefpropagationcache.jl +++ b/src/MessagePassing/abstractbeliefpropagationcache.jl @@ -190,9 +190,6 @@ More generic interface for update, with default params """ function update(alg::Algorithm"bp", bpc::AbstractBeliefPropagationCache) compute_error = !isnothing(alg.kwargs.tolerance) - if !compute_error - println("ALERT, not computing error") - end if isnothing(alg.kwargs.maxiter) error("You need to specify a number of iterations for BP!") end @@ -205,7 +202,7 @@ function update(alg::Algorithm"bp", bpc::AbstractBeliefPropagationCache) update_iteration!(alg, bpc, alg.kwargs.edge_sequence; (update_diff!) = diff) if compute_error diffs[i] = diff.x - if (diff.x / length(alg.kwargs.edge_sequence)) <= alg.kwargs.tolerance + if (diffs[i] / length(alg.kwargs.edge_sequence)) <= alg.kwargs.tolerance if alg.kwargs.verbose println("BP converged to desired precision after $i iterations.") end diff --git a/src/MessagePassing/clustercorrections.jl b/src/MessagePassing/clustercorrections.jl index 63db8cb..20d41ab 100644 --- a/src/MessagePassing/clustercorrections.jl +++ b/src/MessagePassing/clustercorrections.jl @@ -1,5 +1,5 @@ using NamedGraphs -using NamedGraphs: AbstractNamedGraph +using NamedGraphs: AbstractGraph,AbstractNamedGraph struct Cluster loop_ids::Vector{Int} diff --git a/src/MessagePassing/cumulant-clustercorrections.jl b/src/MessagePassing/cumulant-clustercorrections.jl index 8778d8e..822a0f9 100644 --- a/src/MessagePassing/cumulant-clustercorrections.jl +++ b/src/MessagePassing/cumulant-clustercorrections.jl @@ -57,6 +57,18 @@ function induced_components(g::SimpleGraph, vs::AbstractVector{Int})::Vector{Reg return [to_key(vs[c]) for c in comps] end +function is_proper_subset(a::RegionKey, b::RegionKey)::Bool + if length(a) >= length(b) + return false + end + sb = Set(b) + @inbounds for x in a + if x ∉ sb + return false + end + end + return true +end # --- Maximal regions under inclusion ------------------------------------------ """ @@ -72,10 +84,13 @@ function maximal_regions(regions::Set{RegionKey})::Set{RegionKey} is_sub = false for j in (i+1):length(keys) b = keys[j] - if a != b && issubset(a, b) + if length(Set(a)) < length(Set(b)) && issubset(Set(a), Set(b)) + @assert is_proper_subset(a, b) is_sub = true break - end + else + @assert !is_proper_subset(a,b) + end end if !is_sub push!(maximal, a) @@ -192,9 +207,13 @@ function counting_numbers(regions::Set{RegionKey},maximals::Set{RegionKey})::Dic end s = 0 for a in R - if r!=a && issubset(r, a) + # proper subset + if length(Set(r)) < length(Set(a)) && issubset(Set(r), Set(a)) + @assert is_proper_subset(r,a) s += get(c, a, 0) - end + else + @assert !is_proper_subset(r,a) + end end c[r] = 1 - s end diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl index 9f36652..6b26f9c 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/loopcorrection.jl @@ -92,6 +92,7 @@ function weight(bpc::BeliefPropagationCache, eg; project_out::Bool = true, op_st incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, vs)] else incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, es)] + end local_tensors = reduce(vcat, bp_factors(bpc, vs; op_strings = op_strings, coeffs = coeffs, use_epsilon = true)) if project_out diff --git a/src/TensorNetworkQuantumSimulator.jl b/src/TensorNetworkQuantumSimulator.jl index 1373e56..2f69ebb 100644 --- a/src/TensorNetworkQuantumSimulator.jl +++ b/src/TensorNetworkQuantumSimulator.jl @@ -9,6 +9,7 @@ include("TensorNetworks/tensornetwork.jl") include("TensorNetworks/tensornetworkstate.jl") include("TensorNetworks/tensornetworkstate_constructors.jl") include("contraction_sequences.jl") +include("graph_enumeration.jl") include("Forms/bilinearform.jl") include("Forms/quadraticform.jl") include("MessagePassing/abstractbeliefpropagationcache.jl") diff --git a/src/graph_enumeration.jl b/src/graph_enumeration.jl index a8cc524..3c9d915 100644 --- a/src/graph_enumeration.jl +++ b/src/graph_enumeration.jl @@ -4,7 +4,6 @@ using Graphs.Experimental using Graphs.SimpleGraphs using NamedGraphs # from TensorNetworkQuantumSimulator using StatsBase -using ProgressMeter # Generate non-isomorphic graphs with geng # on square lattice, choose triangle_free = true to speed things up From 0356915f710662e36858befed10ad1e9e3d28fba Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Wed, 19 Nov 2025 10:52:21 -0500 Subject: [PATCH 35/36] maybe fixed? --- .../cumulant-clustercorrections.jl | 18 ------------------ src/MessagePassing/loopcorrection.jl | 13 +++++++++---- 2 files changed, 9 insertions(+), 22 deletions(-) diff --git a/src/MessagePassing/cumulant-clustercorrections.jl b/src/MessagePassing/cumulant-clustercorrections.jl index 822a0f9..5a4a55d 100644 --- a/src/MessagePassing/cumulant-clustercorrections.jl +++ b/src/MessagePassing/cumulant-clustercorrections.jl @@ -57,18 +57,6 @@ function induced_components(g::SimpleGraph, vs::AbstractVector{Int})::Vector{Reg return [to_key(vs[c]) for c in comps] end -function is_proper_subset(a::RegionKey, b::RegionKey)::Bool - if length(a) >= length(b) - return false - end - sb = Set(b) - @inbounds for x in a - if x ∉ sb - return false - end - end - return true -end # --- Maximal regions under inclusion ------------------------------------------ """ @@ -85,11 +73,8 @@ function maximal_regions(regions::Set{RegionKey})::Set{RegionKey} for j in (i+1):length(keys) b = keys[j] if length(Set(a)) < length(Set(b)) && issubset(Set(a), Set(b)) - @assert is_proper_subset(a, b) is_sub = true break - else - @assert !is_proper_subset(a,b) end end if !is_sub @@ -209,10 +194,7 @@ function counting_numbers(regions::Set{RegionKey},maximals::Set{RegionKey})::Dic for a in R # proper subset if length(Set(r)) < length(Set(a)) && issubset(Set(r), Set(a)) - @assert is_proper_subset(r,a) s += get(c, a, 0) - else - @assert !is_proper_subset(r,a) end end c[r] = 1 - s diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl index 6b26f9c..82df505 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/loopcorrection.jl @@ -66,11 +66,12 @@ 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} + es::Vector{<:NamedEdge}; vs = [] ) - vs = unique(vcat(src.(es), dst.(es))) - + if isempty(vs) + vs = unique(vcat(src.(es), dst.(es))) + end bpes = NamedEdge[] for v in vs incoming_es = boundary_edges(bpc, [v]; dir = :in) @@ -89,9 +90,13 @@ function weight(bpc::BeliefPropagationCache, eg; project_out::Bool = true, op_st bpc, antiprojectors = sim_edgeinduced_subgraph(bpc, eg) end if isempty(es) - incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, vs)] + + # this is buggy + # incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, vs)] + incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, es; vs=vs)] else incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, es)] + end local_tensors = reduce(vcat, bp_factors(bpc, vs; op_strings = op_strings, coeffs = coeffs, use_epsilon = true)) From d04aefa57eff94652b58d32ede03365e63dae21b Mon Sep 17 00:00:00 2001 From: Grace Sommers Date: Wed, 19 Nov 2025 20:27:18 -0500 Subject: [PATCH 36/36] fixed boundary_edges --- src/MessagePassing/loopcorrection.jl | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl index 82df505..55cc835 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/loopcorrection.jl @@ -66,12 +66,10 @@ 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 = [] - ) + es::Vector{<:NamedEdge}) - if isempty(vs) - vs = unique(vcat(src.(es), dst.(es))) - end + vs = unique(vcat(src.(es), dst.(es))) + bpes = NamedEdge[] for v in vs incoming_es = boundary_edges(bpc, [v]; dir = :in) @@ -90,10 +88,7 @@ function weight(bpc::BeliefPropagationCache, eg; project_out::Bool = true, op_st bpc, antiprojectors = sim_edgeinduced_subgraph(bpc, eg) end if isempty(es) - - # this is buggy - # incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, vs)] - incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, es; vs=vs)] + incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, vs; dir=:in)] else incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, es)]