From 1e754c861f5d6de3496efc7e168f562bfc7563cc Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 27 Oct 2025 11:17:55 -0400 Subject: [PATCH 01/12] QuadarticForm --- src/Backend/boundarympscache.jl | 4 ++-- src/TensorNetworkQuantumSimulator.jl | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Backend/boundarympscache.jl b/src/Backend/boundarympscache.jl index d70c09d..f9cd936 100644 --- a/src/Backend/boundarympscache.jl +++ b/src/Backend/boundarympscache.jl @@ -27,7 +27,7 @@ end default_bp_maxiter(bmps_cache::BoundaryMPSCache) = is_tree(partitions_graph(supergraph(bmps_cache))) ? 1 : 5 function default_message_update_alg(bmps_cache::BoundaryMPSCache) tn = network(bmps_cache) - if tn isa TensorNetworkState || tn isa BilinearForm + if tn isa TensorNetworkState || tn isa BilinearForm || tn isa QuadraticForm return "orthogonal" elseif tn isa ITensorNetwork return "ITensorMPS" @@ -137,7 +137,7 @@ function virtual_index_dimension( end function BoundaryMPSCache( - tn::Union{TensorNetworkState, ITensorNetwork, BilinearForm}, + tn::Union{TensorNetworkState, ITensorNetwork, BilinearForm, QuadraticForm}, mps_bond_dimension::Int; partition_by = "row", gauge_state = true diff --git a/src/TensorNetworkQuantumSimulator.jl b/src/TensorNetworkQuantumSimulator.jl index 32ece60..7b8f9a9 100644 --- a/src/TensorNetworkQuantumSimulator.jl +++ b/src/TensorNetworkQuantumSimulator.jl @@ -6,6 +6,7 @@ include("imports.jl") include("siteinds.jl") include("tensornetworkstate.jl") include("bilinearform.jl") +include("quadraticform.jl") include("Backend/abstractbeliefpropagationcache.jl") include("Backend/beliefpropagationcache.jl") include("Backend/boundarympscache.jl") From cea13460f275924b68b7c39de6aa6bff9a6ceec6 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 27 Oct 2025 11:18:19 -0400 Subject: [PATCH 02/12] QuadarticForm --- src/quadraticform.jl | 57 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 src/quadraticform.jl diff --git a/src/quadraticform.jl b/src/quadraticform.jl new file mode 100644 index 0000000..f351bf3 --- /dev/null +++ b/src/quadraticform.jl @@ -0,0 +1,57 @@ +struct QuadraticForm{V} <: AbstractITensorNetwork{V} + ket::TensorNetworkState{V} + operator::TensorNetworkState{V} +end + +ket(qf::QuadraticForm) = qf.ket +operator(qf::QuadraticForm) = qf.operator +bra(qf::QuadraticForm) = prime(dag(ket(qf))) + +Base.copy(qf::QuadraticForm) = QuadraticForm(copy(qf.ket), copy(qf.operator)) + +#Forward onto the ket +for f in [ + :(ITensorNetworks.underlying_graph), + :(ITensorNetworks.data_graph_type), + :(ITensorNetworks.data_graph), + :(ITensors.datatype), + :(ITensors.NDTensors.scalartype), + :(NamedGraphs.edgeinduced_subgraphs_no_leaves), + ] + @eval begin + function $f(qf::QuadraticForm, args...; kwargs...) + return $f(ket(qf), args...; kwargs...) + end + end +end + +#Constructor, bra is taken to be in the vector space of ket so the dual is taken +function QuadraticForm(ket::TensorNetworkState, f::Function = v -> "I") + sinds = siteinds(ket) + verts = collect(vertices(ket)) + dtype = datatype(ket) + operator_tensors = adapt(dtype).([reduce(prod, ITensor[ITensors.op(f(v), sind) for sind in sinds[v]]) for v in verts]) + operator = TensorNetworkState(verts, operator_tensors) + return QuadraticForm(ket, operator) +end + +function default_message(qf::QuadraticForm, edge::AbstractEdge) + linds = ITensorNetworks.linkinds(qf, edge) + return adapt(datatype(qf))(denseblocks(delta(linds))) +end + +function ITensorNetworks.linkinds(qf::QuadraticForm, edge::NamedEdge) + ket_linds = linkinds(ket(qf), edge) + return Index[ket_linds; linkinds(operator(qf), edge); dag.(prime.(ket_linds))] +end + +function bp_factors(qf::QuadraticForm, verts::Vector) + factors = ITensor[] + for v in verts + qf_v = ket(qf)[v] + append!(factors, ITensor[qf_v, operator(qf)[v], dag(prime(qf_v))]) + end + return factors +end + +bp_factors(qf::QuadraticForm, v) = bp_factors(qf, [v]) From 68214e1bab60316dd5022517cee07c5db3b4778e Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 27 Oct 2025 16:23:16 -0400 Subject: [PATCH 03/12] Quick Bug Fix --- src/sampling.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/sampling.jl b/src/sampling.jl index 2de77b9..399d04e 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -279,9 +279,7 @@ function certify_sample( return (poverq = p_over_q, bitstring = bitstring) end -certify_sample(ψ, logq_and_bitstring::NamedTuple; kwargs...) = only(certify_sample(ψ, [logq_and_bitstring]; kwargs...)) - function certify_samples(ψ::TensorNetworkState, probs_and_bitstrings::Vector{<:NamedTuple}; alg = "boundarymps", kwargs...) algorithm_check(ψ, "sample", alg) - return [certify_sample(Algorithm(alg), ψ, prob_and_bitstring; kwargs...) for prob_and_bitstring in probs_and_bitstrings] + return [certify_sample(Algorithm(alg), ψ, prob_and_bitstring.bitstring, prob_and_bitstring.logq; kwargs...) for prob_and_bitstring in probs_and_bitstrings] end From 9a53f4c7a9d5e15086f16fbae397117fea603ebe Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 28 Oct 2025 10:15:13 -0400 Subject: [PATCH 04/12] Delete unneeded file --- Project.toml | 2 +- src/Backend/boundarymps.jl | 698 -------------------------------- src/Backend/boundarympscache.jl | 9 +- test/loopcorrections.jl | 2 +- 4 files changed, 6 insertions(+), 705 deletions(-) delete mode 100644 src/Backend/boundarymps.jl diff --git a/Project.toml b/Project.toml index c446aa6..08cc1eb 100644 --- a/Project.toml +++ b/Project.toml @@ -3,7 +3,7 @@ uuid = "4de3b72a-362e-43dd-83ff-3f381eda9f9c" authors = ["JoeyT1994 ", "MSRudolph "] description = "A Julia package for simulating quantum circuits and dynamics with tensor networks of near-arbritrary topology." license = "MIT" -version = "0.1.01" +version = "0.1.02" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/src/Backend/boundarymps.jl b/src/Backend/boundarymps.jl deleted file mode 100644 index 8460be4..0000000 --- a/src/Backend/boundarymps.jl +++ /dev/null @@ -1,698 +0,0 @@ -struct BoundaryMPSCache{V, PV, BPC<:AbstractBeliefPropagationCache{V, PV},PG} <: AbstractBeliefPropagationCache{V, PV} - bp_cache::BPC - partitionedplanargraph::PG - maximum_virtual_dimension::Int64 -end - -const _default_boundarymps_update_alg = "orthogonal" -const _default_boundarymps_update_niters = 40 -const _default_boundarymps_update_tolerance = 1e-12 -const _default_boundarymps_update_cutoff = 1e-12 - -function default_boundarymps_update_kwargs(; cache_is_flat = false, kwargs...) - message_update_alg = Algorithm(ITensorNetworks.default_message_update_alg(cache_is_flat)) - return (; message_update_alg, default_message_update_kwargs(; cache_is_flat, kwargs...)...) -end - -ITensorNetworks.default_message_update_alg(cache_is_flat::Bool = false) = cache_is_flat ? "ITensorMPS" : "orthogonal" - -function default_message_update_kwargs(; cache_is_flat = false, cutoff = _default_boundarymps_update_cutoff, kwargs...) - !cache_is_flat && return return (; niters = _default_boundarymps_update_niters, tolerance = _default_boundarymps_update_tolerance) - return (; cutoff = cutoff, kwargs...) -end - -function is_correct_format(bmpsc::BoundaryMPSCache) - _ppg = ppg(bmpsc) - effective_graph = partitions_graph(_ppg) - if !is_ring_graph(effective_graph) && !is_line_graph(effective_graph) - error("Upon partitioning, graph does not form a line or ring: can't run boundary MPS") - end - for pv in partitionvertices(_ppg) - if !is_line_graph(subgraph(_ppg, pv)) - error("There's a partition that does not form a line: can't run boundary MPS") - end - end - return true -end - -default_cache_update_kwargs(alg::Algorithm"boundarymps") = default_boundarymps_update_kwargs() - -ITensorNetworks.default_update_alg(bmpsc::BoundaryMPSCache) = "bp" -function ITensorNetworks.set_default_kwargs(alg::Algorithm"bp", bmpsc::BoundaryMPSCache) - maxiter = get(alg.kwargs, :maxiter, is_tree(partitions_graph(ppg(bmpsc))) ? 1 : nothing) - edge_sequence = get(alg.kwargs, :edge_sequence, pair.(default_edge_sequence(ppg(bmpsc)))) - verbose = get(alg.kwargs, :verbose, false) - tol = get(alg.kwargs, :tol, nothing) - message_update_alg = ITensorNetworks.set_default_kwargs(get(alg.kwargs, :message_update_alg, Algorithm(ITensorNetworks.default_message_update_alg(is_flat(bmpsc))))) - return Algorithm("bp"; tol, message_update_alg, maxiter, edge_sequence, verbose) -end - -ITensorNetworks.default_normalize(alg::Algorithm"orthogonal") = true -function ITensorNetworks.set_default_kwargs(alg::Algorithm"orthogonal") - normalize = get(alg.kwargs, :normalize, ITensorNetworks.default_normalize(alg)) - tolerance = get(alg.kwargs, :tolerance, _default_boundarymps_update_tolerance) - niters = get(alg.kwargs, :niters, _default_boundarymps_update_niters) - return Algorithm("orthogonal"; tolerance, niters, normalize) -end - -ITensorNetworks.default_normalize(alg::Algorithm"ITensorMPS") = true -function ITensorNetworks.set_default_kwargs(alg::Algorithm"ITensorMPS") - cutoff = get(alg.kwargs, :cutoff, _default_boundarymps_update_cutoff) - normalize = get(alg.kwargs, :normalize, ITensorNetworks.default_normalize(alg)) - return Algorithm("ITensorMPS"; cutoff, normalize) -end - -## Frontend functions - -""" - updatecache(bmpsc::BoundaryMPSCache; alg, message_update_kwargs = (; niters, tolerance)) - -Update the MPS messages inside a boundaryMPS-cache. -""" -function updatecache(bmpsc::BoundaryMPSCache, args...; message_update_alg = ITensorNetworks.default_message_update_alg(is_flat(bmpsc)), - message_update_kwargs = default_message_update_kwargs(; cache_is_flat = is_flat(bmpsc), maxdim = maximum_virtual_dimension(bmpsc)), kwargs...) - return update(bmpsc, args...; message_update_alg, message_update_kwargs..., kwargs...) -end - -""" - build_normsqr_bmps_cache(ψ::AbstractITensorNetwork, message_rank::Int64; cache_construction_kwargs = (;), cache_update_kwargs = default_posdef_boundarymps_update_kwargs()) - -Build the Boundary MPS cache for ψIψ and update it appropriately -""" -function build_normsqr_bmps_cache( - ψ::AbstractITensorNetwork, - message_rank::Int64; - cache_construction_kwargs = (;), - cache_update_kwargs = default_boundarymps_update_kwargs(; cache_is_flat = false, maxdim = message_rank), - update_bp_cache = false, - update_cache = true -) - # build the BP cache - ψIψ = build_normsqr_bp_cache(ψ; update_cache = update_bp_cache) - - # convert BP cache to boundary MPS cache, no further update needed - return build_normsqr_bmps_cache( - ψIψ, - message_rank; - cache_construction_kwargs, - cache_update_kwargs, - update_cache - ) -end - -function build_normsqr_bmps_cache( - ψIψ::AbstractBeliefPropagationCache, - message_rank::Int64; - update_cache = true, - cache_construction_kwargs = (;), - cache_update_kwargs = default_boundarymps_update_kwargs(; cache_is_flat = is_flat(ψIψ), maxdim = message_rank), -) - - ψIψ = BoundaryMPSCache(ψIψ; message_rank, cache_construction_kwargs...) - - if update_cache - ψIψ = updatecache(ψIψ; cache_update_kwargs...) - end - - return ψIψ -end - -is_flat(bmpsc::BoundaryMPSCache) = is_flat(bp_cache(bmpsc)) - -## Backend functions -bp_cache(bmpsc::BoundaryMPSCache) = bmpsc.bp_cache -partitionedplanargraph(bmpsc::BoundaryMPSCache) = bmpsc.partitionedplanargraph -ppg(bmpsc) = partitionedplanargraph(bmpsc) -maximum_virtual_dimension(bmpsc::BoundaryMPSCache) = bmpsc.maximum_virtual_dimension -planargraph(bmpsc::BoundaryMPSCache) = unpartitioned_graph(partitionedplanargraph(bmpsc)) - -function ITensorNetworks.partitioned_tensornetwork(bmpsc::BoundaryMPSCache) - return partitioned_tensornetwork(bp_cache(bmpsc)) -end -ITensorNetworks.messages(bmpsc::BoundaryMPSCache) = messages(bp_cache(bmpsc)) - -function ITensorNetworks.default_bp_maxiter( - alg::Algorithm, - bmpsc::BoundaryMPSCache, -) - return default_bp_maxiter(partitions_graph(ppg(bmpsc))) -end -function ITensorNetworks.default_edge_sequence(alg::Algorithm, bmpsc::BoundaryMPSCache) - return pair.(default_edge_sequence(ppg(bmpsc))) -end - -default_boundarymps_message_rank(tn::AbstractITensorNetwork) = maxlinkdim(tn)^2 -ITensorNetworks.partitions(bmpsc::BoundaryMPSCache) = - parent.(collect(partitionvertices(ppg(bmpsc)))) -NamedGraphs.PartitionedGraphs.partitionedges(bmpsc::BoundaryMPSCache) = pair.(partitionedges(ppg(bmpsc))) - -function ITensorNetworks.cache( - alg::Algorithm"boundarymps", - tn; - bp_cache_construction_kwargs = default_cache_construction_kwargs(Algorithm("bp"), tn), - kwargs..., -) - return BoundaryMPSCache( - BeliefPropagationCache(tn; bp_cache_construction_kwargs...); - kwargs..., - ) -end - -function ITensorNetworks.default_cache_construction_kwargs(alg::Algorithm"boundarymps", tn) - return (; - bp_cache_construction_kwargs = default_cache_construction_kwargs( - Algorithm("bp"), - tn, - ) - ) -end - -function Base.copy(bmpsc::BoundaryMPSCache) - return BoundaryMPSCache( - copy(bp_cache(bmpsc)), - copy(ppg(bmpsc)), - maximum_virtual_dimension(bmpsc), - ) -end - -function ITensorNetworks.default_message( - bmpsc::BoundaryMPSCache, - pe::PartitionEdge; - kwargs..., -) - return default_message(bp_cache(bmpsc), pe::PartitionEdge; kwargs...) -end - -#Get the dimension of the virtual index between the two message tensors on pe1 and pe2 -function virtual_index_dimension( - bmpsc::BoundaryMPSCache, - pe1::PartitionEdge, - pe2::PartitionEdge, -) - pes = planargraph_sorted_partitionedges(bmpsc, planargraph_partitionpair(bmpsc, pe1)) - - if findfirst(x -> x == pe1, pes) > findfirst(x -> x == pe2, pes) - lower_pe, upper_pe = pe2, pe1 - else - lower_pe, upper_pe = pe1, pe2 - end - inds_above = reduce(vcat, linkinds.((bmpsc,), partitionedges_above(bmpsc, lower_pe))) - inds_below = reduce(vcat, linkinds.((bmpsc,), partitionedges_below(bmpsc, upper_pe))) - - return Int(minimum(( - prod(Float64.(dim.(inds_above))), - prod(Float64.(dim.(inds_below))), - Float64(maximum_virtual_dimension(bmpsc)), - ))) -end - -#Vertices of the planargraph -function planargraph_vertices(bmpsc::BoundaryMPSCache, partition) - return vertices(ppg(bmpsc), PartitionVertex(partition)) -end - -#Get partition(s) of vertices of the planargraph -function planargraph_partitions(bmpsc::BoundaryMPSCache, vertices::Vector) - return parent.(partitionvertices(ppg(bmpsc), vertices)) -end - -function planargraph_partition(bmpsc::BoundaryMPSCache, vertex) - return only(planargraph_partitions(bmpsc, [vertex])) -end - -#Get interpartition pairs of partition edges in the underlying partitioned tensornetwork -function planargraph_partitionpair(bmpsc::BoundaryMPSCache, pe::PartitionEdge) - return pair(partitionedge(ppg(bmpsc), parent(pe))) -end - -#Sort (bottom to top) partitoonedges between pair of partitions in the planargraph -function planargraph_sorted_partitionedges(bmpsc::BoundaryMPSCache, partitionpair::Pair) - pg = ppg(bmpsc) - src_vs, dst_vs = vertices(pg, PartitionVertex(first(partitionpair))), - vertices(pg, PartitionVertex(last(partitionpair))) - es = reduce( - vcat, - [ - [src_v => dst_v for dst_v in intersect(neighbors(pg, src_v), dst_vs)] for - src_v in src_vs - ], - ) - es = sort(NamedEdge.(es); by = x -> findfirst(isequal(src(x)), src_vs)) - return PartitionEdge.(es) -end - -#Constructor, inserts missing edge in the planar graph to ensure each partition is connected -#allowing the code to work for arbitrary grids and not just square grids -function BoundaryMPSCache( - bpc::BeliefPropagationCache; - grouping_function::Function = v -> first(v), - group_sorting_function::Function = v -> last(v), - message_rank::Int64 = default_boundarymps_message_rank(tensornetwork(bpc)), - -) - bpc = insert_pseudo_planar_edges(bpc; grouping_function) - planar_graph = partitions_graph(bpc) - vertex_groups = group(grouping_function, collect(vertices(planar_graph))) - vertex_groups = map(x -> sort(x; by = group_sorting_function), vertex_groups) - ppg = PartitionedGraph(planar_graph, vertex_groups) - bmpsc = BoundaryMPSCache(bpc, ppg, message_rank) - @assert is_correct_format(bmpsc) - set_interpartition_messages!(bmpsc) - return bmpsc -end - -function BoundaryMPSCache(tn, args...; kwargs...) - return BoundaryMPSCache(BeliefPropagationCache(tn, args...); kwargs...) -end - -#Functions to get the parellel partitionedges sitting above and below a partitionedge -function partitionedges_above(bmpsc::BoundaryMPSCache, pe::PartitionEdge) - pes = planargraph_sorted_partitionedges(bmpsc, planargraph_partitionpair(bmpsc, pe)) - pe_pos = only(findall(x -> x == pe, pes)) - return PartitionEdge[pes[i] for i = (pe_pos+1):length(pes)] -end - -function partitionedges_below(bmpsc::BoundaryMPSCache, pe::PartitionEdge) - pes = planargraph_sorted_partitionedges(bmpsc, planargraph_partitionpair(bmpsc, pe)) - pe_pos = only(findall(x -> x == pe, pes)) - return PartitionEdge[pes[i] for i = 1:(pe_pos-1)] -end - -function partitionedge_above(bmpsc::BoundaryMPSCache, pe::PartitionEdge) - pes_above = partitionedges_above(bmpsc, pe) - isempty(pes_above) && return nothing - return first(pes_above) -end - -function partitionedge_below(bmpsc::BoundaryMPSCache, pe::PartitionEdge) - pes_below = partitionedges_below(bmpsc, pe) - isempty(pes_below) && return nothing - return last(pes_below) -end - -#Initialise all the interpartition message tensors -function set_interpartition_messages!( - bmpsc::BoundaryMPSCache, - partitionpairs::Vector{<:Pair}, -) - m_keys = keys(messages(bmpsc)) - dtype = datatype(bp_cache(bmpsc)) - for partitionpair in partitionpairs - pes = planargraph_sorted_partitionedges(bmpsc, partitionpair) - for pe in pes - if pe ∉ m_keys - m = dense(delta(linkinds(bmpsc, pe))) - set_message!(bmpsc, pe, ITensor[adapt(dtype)(m)]) - end - end - for i = 1:(length(pes)-1) - virt_dim = virtual_index_dimension(bmpsc, pes[i], pes[i+1]) - ind = Index(virt_dim, "m$(i)$(i+1)") - m1, m2 = only(message(bmpsc, pes[i])), only(message(bmpsc, pes[i+1])) - t = adapt(dtype)(dense(delta(ind))) - set_message!(bmpsc, pes[i], ITensor[m1*t]) - set_message!(bmpsc, pes[i+1], ITensor[m2*t]) - end - end - return bmpsc -end - -function set_interpartition_messages!(bmpsc::BoundaryMPSCache) - partitionpairs = pair.(partitionedges(ppg(bmpsc))) - return set_interpartition_messages!( - bmpsc, - vcat(partitionpairs, reverse.(partitionpairs)), - ) -end - -#Switch the message tensors on partition edges with their reverse (and dagger them) -function switch_message!(bmpsc::BoundaryMPSCache, pe::PartitionEdge) - ms = messages(bmpsc) - me, mer = message(bmpsc, pe), message(bmpsc, reverse(pe)) - set!(ms, pe, dag.(mer)) - set!(ms, reverse(pe), dag.(me)) - return bmpsc -end - -function switch_messages!(bmpsc::BoundaryMPSCache, partitionpair::Pair) - for pe in planargraph_sorted_partitionedges(bmpsc, partitionpair) - switch_message!(bmpsc, pe) - end - return bmpsc -end - -function partition_graph(bmpsc::BoundaryMPSCache, partition) - vs = planargraph_vertices(bmpsc, partition) - return subgraph(unpartitioned_graph(ppg(bmpsc)), vs) -end - -function partition_update!(bmpsc::BoundaryMPSCache, seq::Vector{<:PartitionEdge}) - alg = ITensorNetworks.set_default_kwargs(Algorithm("contract", normalize = false)) - for pe in seq - m = updated_message(alg, bp_cache(bmpsc), pe) - set_message!(bmpsc, pe, m) - end - return bmpsc -end - -#Out-of-place version -function partition_update(bmpsc::BoundaryMPSCache, seq::Vector{<:PartitionEdge}) - bmpsc = copy(bmpsc) - return partition_update!(bmpsc, seq) -end - -#Move the orthogonality centre one step on an interpartition from the message tensor on pe1 to that on pe2 -function gauge_step!( - alg::Algorithm"orthogonal", - bmpsc::BoundaryMPSCache, - pe1::PartitionEdge, - pe2::PartitionEdge; - kwargs..., -) - m1, m2 = only(message(bmpsc, pe1)), only(message(bmpsc, pe2)) - @assert !isempty(commoninds(m1, m2)) - left_inds = uniqueinds(m1, m2) - m1, Y = factorize(m1, left_inds; ortho = "left", kwargs...) - m2 = m2 * Y - set_message!(bmpsc, pe1, ITensor[m1]) - set_message!(bmpsc, pe2, ITensor[m2]) - return bmpsc -end - -#Move the orthogonality centre via a sequence of steps between message tensors -function gauge_walk!( - alg::Algorithm, - bmpsc::BoundaryMPSCache, - seq::Vector; - kwargs..., -) - for (pe1, pe2) in seq - gauge_step!(alg::Algorithm, bmpsc, pe1, pe2; kwargs...) - end - return bmpsc -end - -function inserter!( - alg::Algorithm, - bmpsc::BoundaryMPSCache, - update_pe::PartitionEdge, - m::ITensor; -) - set_message!(bmpsc, reverse(update_pe), ITensor[dag(m)]) - return bmpsc -end - -#Default 1-site extracter -function extracter( - alg::Algorithm"orthogonal", - bmpsc::BoundaryMPSCache, - update_pe::PartitionEdge; -) - message_update_alg = ITensorNetworks.set_default_kwargs(Algorithm("contract"; normalize = false)) - m = only(updated_message(message_update_alg, bmpsc,update_pe)) - return m -end - -function ITensors.commonind(bmpsc::BoundaryMPSCache, pe1::PartitionEdge, pe2::PartitionEdge) - m1, m2 = message(bmpsc, pe1), message(bmpsc, pe2) - return commonind(only(m1), only(m2)) -end - -function merge_internal_tensors(O::Union{MPS, MPO}) - internal_inds = filter(i -> isempty(ITensorMPS.siteinds(O, i)), [i for i in 1:length(O)]) - - while !isempty(internal_inds) - site = first(internal_inds) - tensors = [O[i] for i in setdiff([i for i in 1:length(O)], [site])] - if site != length(O) - tensors[site] = tensors[site] * O[site] - else - tensors[site - 1] = tensors[site - 1] * O[site] - end - - O = typeof(O)(tensors) - - internal_inds = filter(i -> isempty(ITensorMPS.siteinds(O, i)), [i for i in 1:length(O)]) - end - return O -end - -function merge_internal_tensors(O::ITensorNetwork) - O = copy(O) - internal_sites = filter(v -> isempty(siteinds(O, v)), collect(vertices(O))) - - while !isempty(internal_sites) - v = first(internal_sites) - sorted_vs = sort(collect(vertices(O))) - v_pos = findfirst(_v -> _v == v, sorted_vs) - vn = v_pos == length(sorted_vs) ? sorted_vs[v_pos - 1] : sorted_vs[v_pos + 1] - O = ITensorNetworks.contract(O, NamedEdge(v => vn)) - O = ITensorNetworks.combine_linkinds(O) - internal_sites = filter(v -> isempty(siteinds(O, v)), collect(vertices(O))) - end - return O -end - -function ITensorMPS.MPO(bmpsc::BoundaryMPSCache, partition) - sorted_vs = sort(planargraph_vertices(bmpsc, partition)) - ts = [copy(bmpsc[v]) for v in sorted_vs] - O = ITensorMPS.MPO(ts) - return O -end - -function ITensorMPS.MPS(bmpsc::BoundaryMPSCache, partitionpair::Pair) - sorted_pes = planargraph_sorted_partitionedges(bmpsc, partitionpair) - ms = [only(message(bmpsc, pe)) for pe in sorted_pes] - return ITensorMPS.MPS(ms) -end - -function truncate!(bmpsc::BoundaryMPSCache, partitionpair::Pair; truncate_kwargs...) - M = ITensorMPS.MPS(bmpsc, partitionpair) - M = ITensorMPS.truncate(M; truncate_kwargs...) - return set_interpartition_message!(bmpsc, M, partitionpair) -end - -function set_interpartition_message!(bmpsc::BoundaryMPSCache, M::Union{MPS, MPO}, partitionpair::Pair) - sorted_pes = planargraph_sorted_partitionedges(bmpsc, partitionpair) - for i in 1:length(M) - set_message!(bmpsc, sorted_pes[i], ITensor[M[i]]) - end - return bmpsc -end - -function updater!(alg::Algorithm"orthogonal", bmpsc::BoundaryMPSCache, partition_graph, prev_pe, update_pe) - prev_pe == nothing && return bmpsc - - gauge_step!(alg, bmpsc, reverse(prev_pe), reverse(update_pe)) - pupdate_seq = a_star(partition_graph, parent(src(prev_pe)), parent(src(update_pe))) - partition_update!(bmpsc, PartitionEdge.(pupdate_seq)) - return bmpsc -end - -function ITensorNetworks.update_message( - alg::Algorithm"orthogonal", bmpsc::BoundaryMPSCache, partitionpair::Pair) - bmpsc = copy(bmpsc) - delete_partition_messages!(bmpsc, first(partitionpair)) - switch_messages!(bmpsc, partitionpair) - pes = planargraph_sorted_partitionedges(bmpsc, partitionpair) - pg = partition_graph(bmpsc, first(partitionpair)) - update_seq = vcat([pes[i] for i in 1:length(pes)], [pes[i] for i in (length(pes) - 1):-1:2]) - - init_gauge_seq = [(reverse(pes[i]), reverse(pes[i-1])) for i in length(pes):-1:2] - init_update_seq = post_order_dfs_edges(pg, parent(src(first(update_seq)))) - !isempty(init_gauge_seq) && gauge_walk!(alg, bmpsc, init_gauge_seq) - !isempty(init_update_seq) && partition_update!(bmpsc, PartitionEdge.(init_update_seq)) - - prev_cf, prev_pe = 0, nothing - for i = 1:alg.kwargs.niters - cf = 0 - if i == alg.kwargs.niters - update_seq = vcat(update_seq, pes[1]) - end - for update_pe in update_seq - updater!(alg, bmpsc, pg, prev_pe, update_pe) - m = extracter(alg, bmpsc, update_pe) - n = norm(m) - cf += n - if alg.kwargs.normalize - m /= n - end - inserter!(alg, bmpsc, update_pe, m) - prev_pe = update_pe - end - epsilon = abs(cf - prev_cf) / length(update_seq) - !isnothing(alg.kwargs.tolerance) && epsilon < alg.kwargs.tolerance && break - prev_cf = cf - end - delete_partition_messages!(bmpsc, first(partitionpair)) - switch_messages!(bmpsc, partitionpair) - return bmpsc -end - -function prev_partitionpair(bmpsc::BoundaryMPSCache, partitionpair::Pair) - pppg = partitions_graph(ppg(bmpsc)) - vns = neighbors(pppg, first(partitionpair)) - length(vns) == 1 && return nothing - - @assert length(vns) == 2 - v1, v2 = first(vns), last(vns) - last(partitionpair) == v1 && return v2 => first(partitionpair) - last(partitionpair) == v2 && return v1 => first(partitionpair) -end - -function generic_apply(O::MPO, M::MPS; normalize = true, kwargs...) - is_simple_mpo = (length(O) == length(M) && all([length(ITensors.siteinds(O, i)) == 2 for i in 1:length(O)])) - - if is_simple_mpo - out = ITensorMPS.apply(O, M; alg = "naive", kwargs...) - if normalize - out = ITensors.normalize(out) - end - return out - end - - O_tensors = ITensor[] - for i in 1:length(O) - m_ind = filter(j -> !isempty(ITensors.commoninds(O[i], M[j])), [j for j in 1:length(M)]) - if isempty(m_ind) - push!(O_tensors, copy(O[i])) - else - m_ind = only(m_ind) - push!(O_tensors, copy(O[i]) * M[m_ind]) - end - end - O = ITensorNetwork([i for i in 1:length(O_tensors)], O_tensors) - - #Transform away edges that make a loop - loop_edges = filter(e -> abs(src(e) - dst(e)) != 1, edges(O)) - for e in loop_edges - edge_to_split = e - inbetween_vertices = [i for i in (minimum((src(e), dst(e)))+1):(maximum((src(e), dst(e)))-1)] - for v in inbetween_vertices - edge_to_split_ind = only(linkinds(O, edge_to_split)) - O = ITensorNetworks.split_index(O, [edge_to_split]) - d = adapt(datatype(O[v]))(denseblocks(delta(edge_to_split_ind, edge_to_split_ind'))) - O[v] *= d - edge_to_split = NamedEdge(v => maximum((src(e), dst(e)))) - end - end - - O = ITensorNetworks.combine_linkinds(O) - @assert is_tree(O) - O = ITensorMPS.MPS([O[v] for v in vertices(O)]) - O = merge_internal_tensors(O) - - if normalize - O = ITensors.normalize(O) - end - - return truncate(O; kwargs...) -end - -#Update all the message tensors on an interpartition via the ITensorMPS apply function -function ITensorNetworks.update_message( - alg::Algorithm"ITensorMPS", - bmpsc::BoundaryMPSCache, - partitionpair::Pair; - maxdim::Int = maximum_virtual_dimension(bmpsc), -) - bmpsc = copy(bmpsc) - prev_pp = prev_partitionpair(bmpsc, partitionpair) - O = ITensorMPS.MPO(bmpsc, first(partitionpair)) - O = ITensorMPS.truncate(O; alg.kwargs.cutoff, maxdim) - if isnothing(prev_pp) - O = merge_internal_tensors(O) - if alg.kwargs.normalize - O = ITensors.normalize(O) - end - return set_interpartition_message!(bmpsc, O, partitionpair) - end - - M = ITensorMPS.MPS(bmpsc, prev_pp) - M_out = generic_apply(O, M; cutoff = alg.kwargs.cutoff, normalize = alg.kwargs.normalize, maxdim) - return set_interpartition_message!(bmpsc, M_out, partitionpair) -end - -#Environment support, assume all vertices live in the same partition for now -function ITensorNetworks.environment(bmpsc::BoundaryMPSCache, verts::Vector; kwargs...) - vs = parent.((partitionvertices(bp_cache(bmpsc), verts))) - partition = only(planargraph_partitions(bmpsc, parent.(partitionvertices(bmpsc, verts)))) - pg = partition_graph(bmpsc, partition) - update_seq = post_order_dfs_edges(pg,first(vs)) - bmpsc = partition_update(bmpsc, PartitionEdge.(update_seq)) - return environment(bp_cache(bmpsc), verts; kwargs...) -end - -function ITensorNetworks.region_scalar(bmpsc::BoundaryMPSCache, partition) - pg = partition_graph(bmpsc, partition) - v = first(center(pg)) - update_seq = post_order_dfs_edges(pg,v) - bmpsc = partition_update(bmpsc, PartitionEdge.(update_seq)) - return region_scalar(bp_cache(bmpsc), PartitionVertex(v)) -end - -function ITensorNetworks.region_scalar(bmpsc::BoundaryMPSCache, verts::Vector) - partitions = planargraph_partitions(bmpsc, parent.(partitionvertices(bmpsc, verts))) - if length(partitions) == 1 - return region_scalar(bmpsc, only(partitions)) - end - error("Contractions involving more than 1 partition not currently supported") -end - -function ITensorNetworks.region_scalar(bmpsc::BoundaryMPSCache, partitionpair::Pair) - pes = planargraph_sorted_partitionedges(bmpsc, partitionpair) - out = ITensor(one(Bool)) - for pe in pes - out = (out * (only(message(bmpsc, pe)))) * only(message(bmpsc, reverse(pe))) - end - return out[] -end - -function add_partitionedges(pg::PartitionedGraph, pes::Vector{<:PartitionEdge}) - g = partitions_graph(pg) - g = add_edges(g, parent.(pes)) - return PartitionedGraph( - unpartitioned_graph(pg), - g, - partitioned_vertices(pg), - which_partition(pg), - ) -end - -function add_partitionedges(bpc::BeliefPropagationCache, pes::Vector{<:PartitionEdge}) - pg = add_partitionedges(partitioned_tensornetwork(bpc), pes) - return BeliefPropagationCache(pg, messages(bpc)) -end - -#Add partition edges necessary to connect up all vertices in a partition in the planar graph created by the sort function -function insert_pseudo_planar_edges( - bpc::BeliefPropagationCache; - grouping_function = v -> first(v), -) - pg = partitions_graph(bpc) - partitions = unique(grouping_function.(collect(vertices(pg)))) - pseudo_edges = PartitionEdge[] - for p in partitions - vs = sort(filter(v -> grouping_function(v) == p, collect(vertices(pg)))) - for i = 1:(length(vs)-1) - if vs[i] ∉ neighbors(pg, vs[i+1]) - push!(pseudo_edges, PartitionEdge(NamedEdge(vs[i] => vs[i+1]))) - end - end - end - return add_partitionedges(bpc, pseudo_edges) -end - -pair(pe::PartitionEdge) = parent(src(pe)) => parent(dst(pe)) - -function delete_partition_messages!(bmpsc::BoundaryMPSCache, partition) - pg = partition_graph(bmpsc, partition) - pes = PartitionEdge.(edges(pg)) - pes = vcat(pes, reverse.(pes)) - return delete_messages!(bmpsc, filter(pe -> pe ∈ keys(messages(bmpsc)), pes)) -end - -function delete_partitionpair_messages!(bmpsc::BoundaryMPSCache, partitionpair::Pair) - pes = planargraph_sorted_partitionedges(bmpsc, partitionpair) - return delete_messages!(bmpsc, filter(pe -> pe ∈ keys(messages(bmpsc)), pes)) -end diff --git a/src/Backend/boundarympscache.jl b/src/Backend/boundarympscache.jl index f9cd936..13b7f64 100644 --- a/src/Backend/boundarympscache.jl +++ b/src/Backend/boundarympscache.jl @@ -25,16 +25,15 @@ function default_bp_edge_sequence(bmps_cache::BoundaryMPSCache) return PartitionEdge.(forest_cover_edge_sequence(partitions_graph(supergraph(bmps_cache)))) end default_bp_maxiter(bmps_cache::BoundaryMPSCache) = is_tree(partitions_graph(supergraph(bmps_cache))) ? 1 : 5 -function default_message_update_alg(bmps_cache::BoundaryMPSCache) - tn = network(bmps_cache) - if tn isa TensorNetworkState || tn isa BilinearForm || tn isa QuadraticForm +function default_bmps_message_update_alg(tn) + if tn isa TensorNetworkState || tn isa BilinearForm || tn isa QuadraticForm return "orthogonal" elseif tn isa ITensorNetwork return "ITensorMPS" - else - return error("Unrecognized network type inside the cache. Don't know what message update alg to use.") end + return error("Unrecognized network type. Don't know what BMPS message update alg to use.") end +default_message_update_alg(bmps_cache::BoundaryMPSCache) = default_bmps_message_update_alg(network(bmps_cache)) default_normalize(alg::Algorithm"orthogonal") = true default_tolerance(bmps_cache::BoundaryMPSCache) = default_tolerance(ITensors.NDTensors.scalartype(bmps_cache)) diff --git a/test/loopcorrections.jl b/test/loopcorrections.jl index 452cd40..925e9a2 100644 --- a/test/loopcorrections.jl +++ b/test/loopcorrections.jl @@ -18,7 +18,7 @@ Random.seed!(1634) @testset "Loop Corrections" begin nx, ny = 4, 4 - χ = 3 + χ = 2 ITensors.disable_warn_order() gs = [ (named_grid((nx, 1)), "line", 0), From 32c371753e22096d03e745e737152dc902a2c076 Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 28 Oct 2025 10:16:02 -0400 Subject: [PATCH 05/12] Runic formatting --- src/Backend/boundarympscache.jl | 2 +- src/apply.jl | 70 ++++++++++----------- src/gates.jl | 12 ++-- src/quadraticform.jl | 2 +- src/sample.jl | 106 ++++++++++++++++---------------- test/loopcorrections.jl | 6 +- 6 files changed, 100 insertions(+), 98 deletions(-) diff --git a/src/Backend/boundarympscache.jl b/src/Backend/boundarympscache.jl index 13b7f64..53e53a4 100644 --- a/src/Backend/boundarympscache.jl +++ b/src/Backend/boundarympscache.jl @@ -26,7 +26,7 @@ function default_bp_edge_sequence(bmps_cache::BoundaryMPSCache) end default_bp_maxiter(bmps_cache::BoundaryMPSCache) = is_tree(partitions_graph(supergraph(bmps_cache))) ? 1 : 5 function default_bmps_message_update_alg(tn) - if tn isa TensorNetworkState || tn isa BilinearForm || tn isa QuadraticForm + if tn isa TensorNetworkState || tn isa BilinearForm || tn isa QuadraticForm return "orthogonal" elseif tn isa ITensorNetwork return "ITensorMPS" diff --git a/src/apply.jl b/src/apply.jl index 4e43378..2feef6b 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -38,15 +38,15 @@ function adapt_gate(gate::ITensor, ψ_bpc::BeliefPropagationCache) end function apply_gates( - circuit::Vector{<:ITensor}, - ψ_bpc::BeliefPropagationCache; - gate_vertices::Vector = neighbor_vertices.((network(ψ_bpc), ), circuit), - apply_kwargs = (; ), - bp_update_kwargs = default_bp_update_kwargs(ψ_bpc), - update_cache = true, - verbose = false, - transfer_to_gpu = false, -) + circuit::Vector{<:ITensor}, + ψ_bpc::BeliefPropagationCache; + gate_vertices::Vector = neighbor_vertices.((network(ψ_bpc),), circuit), + apply_kwargs = (;), + bp_update_kwargs = default_bp_update_kwargs(ψ_bpc), + update_cache = true, + verbose = false, + transfer_to_gpu = false, + ) ψ_bpc = copy(ψ_bpc) # we keep track of the vertices that have been acted on by 2-qubit gates @@ -92,13 +92,13 @@ end #Apply function for a single gate function apply_gate!( - gate::ITensor, - ψ_bpc::BeliefPropagationCache; - v⃗ = ITensorNetworks.neighbor_vertices(ψ_bpc, gate), - apply_kwargs = _default_apply_kwargs, - transfer_to_gpu = false, -) - envs = length(v⃗) == 1 ? nothing : incoming_messages(ψ_bpc,v⃗) + gate::ITensor, + ψ_bpc::BeliefPropagationCache; + v⃗ = ITensorNetworks.neighbor_vertices(ψ_bpc, gate), + apply_kwargs = _default_apply_kwargs, + transfer_to_gpu = false, + ) + envs = length(v⃗) == 1 ? nothing : incoming_messages(ψ_bpc, v⃗) if !transfer_to_gpu updated_tensors, s_values, err = simple_update(gate, network(ψ_bpc), v⃗; envs, apply_kwargs...) @@ -177,8 +177,8 @@ function simple_update( end function simple_update_cuda( - o::ITensor, ψ, v⃗; envs, normalize_tensors = true, apply_kwargs... - ) + o::ITensor, ψ, v⃗; envs, normalize_tensors = true, apply_kwargs... + ) dtype = datatype(ψ[first(v⃗)]) o = CUDA.cu(o) @@ -193,24 +193,24 @@ function simple_update_cuda( envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) sqrt_envs_v1 = [ - ITensorNetworks.ITensorsExtensions.map_eigvals( - sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian=true - ) for env in envs_v1 + ITensorNetworks.ITensorsExtensions.map_eigvals( + sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian = true + ) for env in envs_v1 ] sqrt_envs_v2 = [ ITensorNetworks.ITensorsExtensions.map_eigvals( - sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian=true - ) for env in envs_v2 + sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian = true + ) for env in envs_v2 ] inv_sqrt_envs_v1 = [ ITensorNetworks.ITensorsExtensions.map_eigvals( - inv ∘ sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian=true - ) for env in envs_v1 + inv ∘ sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian = true + ) for env in envs_v1 ] inv_sqrt_envs_v2 = [ ITensorNetworks.ITensorsExtensions.map_eigvals( - inv ∘ sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian=true - ) for env in envs_v2 + inv ∘ sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian = true + ) for env in envs_v2 ] ψᵥ₁ = contract([ψ1; sqrt_envs_v1]) ψᵥ₂ = contract([ψ2; sqrt_envs_v2]) @@ -224,12 +224,12 @@ function simple_update_cuda( e = v⃗[1] => v⃗[2] singular_values! = Ref(ITensor()) Rᵥ₁, Rᵥ₂, spec = factorize_svd( - oR, - unioninds(rᵥ₁, sᵥ₁); - ortho="none", - tags=edge_tag(e), - singular_values!, - apply_kwargs..., + oR, + unioninds(rᵥ₁, sᵥ₁); + ortho = "none", + tags = edge_tag(e), + singular_values!, + apply_kwargs..., ) err = spec.truncerr s_values = singular_values![] @@ -248,8 +248,8 @@ function simple_update_cuda( return noprime.(updated_tensors), s_values, err end - #Checker for whether the cache needs updating (overlapping gate encountered) -function _cacheupdate_check(affected_indices::Set, gate::ITensor; inds_per_site=1) +#Checker for whether the cache needs updating (overlapping gate encountered) +function _cacheupdate_check(affected_indices::Set, gate::ITensor; inds_per_site = 1) indices = inds(gate) # check if we have a two-site gate and any of the qinds are in the affected_indices. If so update cache diff --git a/src/gates.jl b/src/gates.jl index 61f1ea6..4caae96 100644 --- a/src/gates.jl +++ b/src/gates.jl @@ -160,9 +160,9 @@ end Gate for rotation by Z+ at a given halved angle """ function ITensors.op( - ::OpName"Rz+", ::SiteType"S=1/2"; θ::Number - ) - a = exp( -im * θ * 0.5) + ::OpName"Rz+", ::SiteType"S=1/2"; θ::Number + ) + a = exp(-im * θ * 0.5) mat = zeros(ComplexF64, 2, 2) mat[1, 1] = 1 mat[2, 2] = a @@ -176,9 +176,9 @@ end Gate for rotation by Z+Z+ at a given halved angle """ function ITensors.op( - ::OpName"Rz+z+", ::SiteType"S=1/2"; θ::Number - ) - a = exp( -im * θ * 0.5) + ::OpName"Rz+z+", ::SiteType"S=1/2"; θ::Number + ) + a = exp(-im * θ * 0.5) mat = zeros(ComplexF64, 4, 4) mat[1, 1] = 1 mat[2, 2] = 1 diff --git a/src/quadraticform.jl b/src/quadraticform.jl index f351bf3..1dd9021 100644 --- a/src/quadraticform.jl +++ b/src/quadraticform.jl @@ -26,7 +26,7 @@ for f in [ end #Constructor, bra is taken to be in the vector space of ket so the dual is taken -function QuadraticForm(ket::TensorNetworkState, f::Function = v -> "I") +function QuadraticForm(ket::TensorNetworkState, f::Function = v -> "I") sinds = siteinds(ket) verts = collect(vertices(ket)) dtype = datatype(ket) diff --git a/src/sample.jl b/src/sample.jl index df40e79..3a53caf 100644 --- a/src/sample.jl +++ b/src/sample.jl @@ -1,38 +1,39 @@ using StatsBase function _sample( - ψ::ITensorNetwork, - nsamples::Int64; - projected_message_rank::Int64, - norm_message_rank::Int64, - norm_message_update_kwargs=(; niters = _default_boundarymps_update_niters, tolerance = _default_boundarymps_update_tolerance), - projected_message_update_kwargs = (;cutoff = _default_boundarymps_update_cutoff, maxdim = projected_message_rank), - partition_by = "Row", - kwargs..., -) + ψ::ITensorNetwork, + nsamples::Int64; + projected_message_rank::Int64, + norm_message_rank::Int64, + norm_message_update_kwargs = (; niters = _default_boundarymps_update_niters, tolerance = _default_boundarymps_update_tolerance), + projected_message_update_kwargs = (; cutoff = _default_boundarymps_update_cutoff, maxdim = projected_message_rank), + partition_by = "Row", + kwargs..., + ) grouping_function = partition_by == "Column" ? v -> last(v) : v -> first(v) group_sorting_function = partition_by == "Column" ? v -> first(v) : v -> last(v) ψ, ψψ = symmetric_gauge(ψ) ψ, ψψ = normalize(ψ, ψψ) - norm_MPScache = BoundaryMPSCache(ψψ; message_rank=norm_message_rank, grouping_function, group_sorting_function) + norm_MPScache = BoundaryMPSCache(ψψ; message_rank = norm_message_rank, grouping_function, group_sorting_function) sorted_partitions = sort(ITensorNetworks.partitions(norm_MPScache)) seq = [ - sorted_partitions[i] => sorted_partitions[i-1] for - i = length(sorted_partitions):-1:2 + sorted_partitions[i] => sorted_partitions[i - 1] for + i in length(sorted_partitions):-1:2 ] - norm_message_update_kwargs = (; norm_message_update_kwargs..., normalize=false) + norm_message_update_kwargs = (; norm_message_update_kwargs..., normalize = false) norm_MPScache = ITensorNetworks.update_iteration(Algorithm("bp"; message_update_alg = Algorithm("orthogonal"; norm_message_update_kwargs...)), norm_MPScache, seq) - projected_MPScache = BoundaryMPSCache(ψ; message_rank=projected_message_rank, grouping_function, group_sorting_function) + projected_MPScache = BoundaryMPSCache(ψ; message_rank = projected_message_rank, grouping_function, group_sorting_function) #Generate the bit_strings moving left to right through the network probs_and_bitstrings = NamedTuple[] - for j = 1:nsamples + for j in 1:nsamples p_over_q_approx, logq, bitstring = _get_one_sample( - norm_MPScache, projected_MPScache, sorted_partitions; projected_message_update_kwargs, kwargs...) - push!(probs_and_bitstrings, (poverq=p_over_q_approx, logq=logq, bitstring=bitstring)) + norm_MPScache, projected_MPScache, sorted_partitions; projected_message_update_kwargs, kwargs... + ) + push!(probs_and_bitstrings, (poverq = p_over_q_approx, logq = logq, bitstring = bitstring)) end return probs_and_bitstrings, ψ @@ -73,7 +74,7 @@ end Take nsamples bitstrings from a 2D open boundary tensornetwork by partitioning it and using boundary MPS algorithm with relevant ranks. Returns a vector of (p/q, logq, bitstring) where loqq is log probability of drawing the bitstring and p/q attests to the quality of the bitstring which is accurate only if the projected boundary MPS rank is high enough. """ -function sample_directly_certified(ψ::ITensorNetwork, nsamples::Int64; projected_message_rank=5 * maxlinkdim(ψ), kwargs...) +function sample_directly_certified(ψ::ITensorNetwork, nsamples::Int64; projected_message_rank = 5 * maxlinkdim(ψ), kwargs...) probs_and_bitstrings, _ = _sample(ψ::ITensorNetwork, nsamples::Int64; projected_message_rank, kwargs...) # returns the self-certified p/q, logq and bitstrings return probs_and_bitstrings @@ -95,25 +96,25 @@ Take nsamples bitstrings from a 2D open boundary tensornetwork by partitioning i an independent contraction of to get a measure of p/q. Returns a vector of (p/q, bitstring) where p/q attests to the quality of the bitstring which is accurate only if the certification boundary MPS rank is high enough. """ -function sample_certified(ψ::ITensorNetwork, nsamples::Int64; certification_message_rank=5 * maxlinkdim(ψ), certification_message_update_kwargs = (; cutoff = _default_boundarymps_update_cutoff), kwargs...) +function sample_certified(ψ::ITensorNetwork, nsamples::Int64; certification_message_rank = 5 * maxlinkdim(ψ), certification_message_update_kwargs = (; cutoff = _default_boundarymps_update_cutoff), kwargs...) probs_and_bitstrings, ψ = _sample(ψ::ITensorNetwork, nsamples::Int64; kwargs...) # send the bitstrings and the logq to the certification function - return certify_samples(ψ, probs_and_bitstrings; certification_message_rank, certification_message_update_kwargs, symmetrize_and_normalize=false) + return certify_samples(ψ, probs_and_bitstrings; certification_message_rank, certification_message_update_kwargs, symmetrize_and_normalize = false) end function _get_one_sample( - norm_MPScache::BoundaryMPSCache, - projected_MPScache::BoundaryMPSCache, - sorted_partitions; - projected_message_update_kwargs= (; cutoff = _default_boundarymps_update_cutoff, maxdim = maximum_virtual_dimension(projected_MPScache)), - kwargs..., -) + norm_MPScache::BoundaryMPSCache, + projected_MPScache::BoundaryMPSCache, + sorted_partitions; + projected_message_update_kwargs = (; cutoff = _default_boundarymps_update_cutoff, maxdim = maximum_virtual_dimension(projected_MPScache)), + kwargs..., + ) - projected_message_update_kwargs = (; projected_message_update_kwargs..., normalize=false) + projected_message_update_kwargs = (; projected_message_update_kwargs..., normalize = false) norm_MPScache = copy(norm_MPScache) - bit_string = Dictionary{keytype(vertices(projected_MPScache)),Int}() + bit_string = Dictionary{keytype(vertices(projected_MPScache)), Int}() p_over_q_approx = nothing logq = 0 for (i, partition) in enumerate(sorted_partitions) @@ -130,12 +131,14 @@ function _get_one_sample( if i < length(sorted_partitions) - next_partition = sorted_partitions[i+1] - + next_partition = sorted_partitions[i + 1] + #Alternate fitting procedure here which is faster for small bond dimensions but slower for large - projected_MPScache = ITensorNetworks.update_message(ITensorNetworks.set_default_kwargs(Algorithm("ITensorMPS"; projected_message_update_kwargs...)), + projected_MPScache = ITensorNetworks.update_message( + ITensorNetworks.set_default_kwargs(Algorithm("ITensorMPS"; projected_message_update_kwargs...)), projected_MPScache, - partition => next_partition) + partition => next_partition + ) pes = planargraph_sorted_partitionedges(norm_MPScache, partition => next_partition) @@ -145,25 +148,25 @@ function _get_one_sample( end end - i > 1 && delete_partitionpair_messages!(projected_MPScache, sorted_partitions[i-1] => sorted_partitions[i]) - i > 2 && delete_partitionpair_messages!(norm_MPScache, sorted_partitions[i-2] => sorted_partitions[i-1]) + i > 1 && delete_partitionpair_messages!(projected_MPScache, sorted_partitions[i - 1] => sorted_partitions[i]) + i > 2 && delete_partitionpair_messages!(norm_MPScache, sorted_partitions[i - 2] => sorted_partitions[i - 1]) end return p_over_q_approx, logq, bit_string end function certify_sample( - ψ::ITensorNetwork, bitstring, logq::Number; - certification_message_rank::Int64, - certification_message_update_kwargs = (; maxdim = certification_message_rank, cutoff = _default_boundarymps_update_cutoff), - symmetrize_and_normalize=true, -) + ψ::ITensorNetwork, bitstring, logq::Number; + certification_message_rank::Int64, + certification_message_update_kwargs = (; maxdim = certification_message_rank, cutoff = _default_boundarymps_update_cutoff), + symmetrize_and_normalize = true, + ) if symmetrize_and_normalize ψ, ψψ = symmetric_gauge(ψ) ψ = normalize(ψ, cache! = Ref(ψψ)) end - certification_message_update_kwargs = (; certification_message_update_kwargs..., normalize=false) + certification_message_update_kwargs = (; certification_message_update_kwargs..., normalize = false) ψproj = copy(ψ) s = siteinds(ψ) @@ -173,7 +176,7 @@ function certify_sample( ψproj[v] = ψproj[v] * P * inv(qv) end - bmpsc = BoundaryMPSCache(ψproj; message_rank=certification_message_rank) + bmpsc = BoundaryMPSCache(ψproj; message_rank = certification_message_rank) certification_message_update_kwargs = (; normalize = false, certification_message_update_kwargs...) #This block is two times faster than the two lines below but likely less accurate for smaller maxdims @@ -189,7 +192,7 @@ function certify_sample( p_over_q *= conj(p_over_q) - return (poverq=p_over_q, bitstring=bitstring) + return (poverq = p_over_q, bitstring = bitstring) end certify_sample(ψ, logq_and_bitstring::NamedTuple; kwargs...) = certify_sample(ψ, logq_and_bitstring.bitstring, logq_and_bitstring.logq; kwargs...) @@ -204,25 +207,25 @@ end #Sample along the column/ row specified by pv with the left incoming MPS message input and the right extractable from the cache function sample_partition!( - ψIψ::BoundaryMPSCache, - partition, - bit_string::Dictionary; - kwargs..., -) + ψIψ::BoundaryMPSCache, + partition, + bit_string::Dictionary; + kwargs..., + ) vs = sort(planargraph_vertices(ψIψ, partition)) - seq = PartitionEdge[PartitionEdge(vs[i] => vs[i-1]) for i in length(vs):-1:2] + seq = PartitionEdge[PartitionEdge(vs[i] => vs[i - 1]) for i in length(vs):-1:2] !isempty(seq) && partition_update!(ψIψ, seq) prev_v, traces = nothing, [] logq = 0 for v in vs !isnothing(prev_v) && partition_update!(ψIψ, [PartitionEdge(prev_v => v)]) env = environment(bp_cache(ψIψ), [(v, "operator")]) - seq = contraction_sequence(env; alg="optimal") - ρ = contract(env; sequence=seq) + seq = contraction_sequence(env; alg = "optimal") + ρ = contract(env; sequence = seq) ρ_tr = tr(ρ) push!(traces, ρ_tr) ρ *= inv(ρ_tr) - ρ_diag = collect(real.(diag(ITensors.array(ρ)))) + ρ_diag = collect(real.(diag(ITensors.array(ρ)))) config = StatsBase.sample(1:length(ρ_diag), Weights(ρ_diag)) # config is 1 or 2, but we want 0 or 1 for the sample itself set!(bit_string, v, config - 1) @@ -242,4 +245,3 @@ function sample_partition!( return first(traces), logq, bit_string end - diff --git a/test/loopcorrections.jl b/test/loopcorrections.jl index 925e9a2..ae5f9a7 100644 --- a/test/loopcorrections.jl +++ b/test/loopcorrections.jl @@ -34,7 +34,7 @@ Random.seed!(1634) norm_loopcorrected = norm(ψ; alg = "loopcorrections", max_configuration_size = 2 * (smallest_loop_size) - 1) norm_exact = norm(ψ; alg = "exact") - @test isapprox(norm_bp, norm_exact, atol = 5e-2) - @test isapprox(norm_loopcorrected, norm_exact, atol = 5e-2) + @test isapprox(norm_bp, norm_exact, atol = 5.0e-2) + @test isapprox(norm_loopcorrected, norm_exact, atol = 5.0e-2) end -end \ No newline at end of file +end From 1009cb34501de729562aff6e1247b44e6e9a124b Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 28 Oct 2025 10:18:26 -0400 Subject: [PATCH 06/12] Remove leftover files --- src/apply.jl | 84 ++-------------------------------------------------- src/gates.jl | 34 --------------------- 2 files changed, 3 insertions(+), 115 deletions(-) diff --git a/src/apply.jl b/src/apply.jl index 2feef6b..977fde6 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -45,7 +45,6 @@ function apply_gates( bp_update_kwargs = default_bp_update_kwargs(ψ_bpc), update_cache = true, verbose = false, - transfer_to_gpu = false, ) ψ_bpc = copy(ψ_bpc) @@ -79,7 +78,7 @@ function apply_gates( # actually apply the gate gate = adapt_gate(gate, ψ_bpc) - t = @timed ψ_bpc, truncation_errors[ii] = apply_gate!(gate, ψ_bpc; v⃗ = gate_vertices[ii], apply_kwargs, transfer_to_gpu) + t = @timed ψ_bpc, truncation_errors[ii] = apply_gate!(gate, ψ_bpc; v⃗ = gate_vertices[ii], apply_kwargs) affected_vertices = union(affected_vertices, Set(gate_vertices[ii])) end @@ -95,16 +94,11 @@ function apply_gate!( gate::ITensor, ψ_bpc::BeliefPropagationCache; v⃗ = ITensorNetworks.neighbor_vertices(ψ_bpc, gate), - apply_kwargs = _default_apply_kwargs, - transfer_to_gpu = false, + apply_kwargs = _default_apply_kwargs ) envs = length(v⃗) == 1 ? nothing : incoming_messages(ψ_bpc, v⃗) - if !transfer_to_gpu - updated_tensors, s_values, err = simple_update(gate, network(ψ_bpc), v⃗; envs, apply_kwargs...) - else - updated_tensors, s_values, err = simple_update_cuda(gate, network(ψ_bpc), v⃗; envs, apply_kwargs...) - end + updated_tensors, s_values, err = simple_update(gate, network(ψ_bpc), v⃗; envs, apply_kwargs...) if length(v⃗) == 2 v1, v2 = v⃗ @@ -176,78 +170,6 @@ function simple_update( return noprime.(updated_tensors), s_values, err end -function simple_update_cuda( - o::ITensor, ψ, v⃗; envs, normalize_tensors = true, apply_kwargs... - ) - - dtype = datatype(ψ[first(v⃗)]) - o = CUDA.cu(o) - if length(v⃗) == 1 - ψ1 = CUDA.cu(ψ[first(v⃗)]) - updated_tensors = ITensor[ITensors.apply(o, ψ1)] - s_values, err = nothing, 0 - else - ψ1 = CUDA.cu(ψ[first(v⃗)]) - ψ2 = CUDA.cu(ψ[last(v⃗)]) - cutoff = 10 * eps(real(scalartype(ψ[v⃗[1]]))) - envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs) - envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs) - sqrt_envs_v1 = [ - ITensorNetworks.ITensorsExtensions.map_eigvals( - sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian = true - ) for env in envs_v1 - ] - sqrt_envs_v2 = [ - ITensorNetworks.ITensorsExtensions.map_eigvals( - sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian = true - ) for env in envs_v2 - ] - inv_sqrt_envs_v1 = [ - ITensorNetworks.ITensorsExtensions.map_eigvals( - inv ∘ sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian = true - ) for env in envs_v1 - ] - inv_sqrt_envs_v2 = [ - ITensorNetworks.ITensorsExtensions.map_eigvals( - inv ∘ sqrt, CUDA.cu(env), inds(env)[1], inds(env)[2]; cutoff, ishermitian = true - ) for env in envs_v2 - ] - ψᵥ₁ = contract([ψ1; sqrt_envs_v1]) - ψᵥ₂ = contract([ψ2; sqrt_envs_v2]) - sᵥ₁ = commoninds(ψ1, o) - sᵥ₂ = commoninds(ψ2, o) - Qᵥ₁, Rᵥ₁ = qr(ψᵥ₁, uniqueinds(uniqueinds(ψᵥ₁, ψᵥ₂), sᵥ₁)) - Qᵥ₂, Rᵥ₂ = qr(ψᵥ₂, uniqueinds(uniqueinds(ψᵥ₂, ψᵥ₁), sᵥ₂)) - rᵥ₁ = commoninds(Qᵥ₁, Rᵥ₁) - rᵥ₂ = commoninds(Qᵥ₂, Rᵥ₂) - oR = ITensors.apply(o, Rᵥ₁ * Rᵥ₂) - e = v⃗[1] => v⃗[2] - singular_values! = Ref(ITensor()) - Rᵥ₁, Rᵥ₂, spec = factorize_svd( - oR, - unioninds(rᵥ₁, sᵥ₁); - ortho = "none", - tags = edge_tag(e), - singular_values!, - apply_kwargs..., - ) - err = spec.truncerr - s_values = singular_values![] - s_values = adapt(dtype)(s_values) - Qᵥ₁ = contract([Qᵥ₁; dag.(inv_sqrt_envs_v1)]) - Qᵥ₂ = contract([Qᵥ₂; dag.(inv_sqrt_envs_v2)]) - updated_tensors = [Qᵥ₁ * Rᵥ₁, Qᵥ₂ * Rᵥ₂] - end - - if normalize_tensors - updated_tensors = ITensor[ψᵥ / norm(ψᵥ) for ψᵥ in updated_tensors] - end - - updated_tensors = map(adapt(dtype), updated_tensors) - - return noprime.(updated_tensors), s_values, err -end - #Checker for whether the cache needs updating (overlapping gate encountered) function _cacheupdate_check(affected_indices::Set, gate::ITensor; inds_per_site = 1) indices = inds(gate) diff --git a/src/gates.jl b/src/gates.jl index 4caae96..09009b8 100644 --- a/src/gates.jl +++ b/src/gates.jl @@ -153,39 +153,5 @@ function ITensors.op( return mat end - -""" - ITensors.op(::OpName"Rz+", ::SiteType"S=1/2"; θ::Number) - -Gate for rotation by Z+ at a given halved angle -""" -function ITensors.op( - ::OpName"Rz+", ::SiteType"S=1/2"; θ::Number - ) - a = exp(-im * θ * 0.5) - mat = zeros(ComplexF64, 2, 2) - mat[1, 1] = 1 - mat[2, 2] = a - return mat -end - - -""" - ITensors.op(::OpName"Rz+z+", ::SiteType"S=1/2"; θ::Number) - -Gate for rotation by Z+Z+ at a given halved angle -""" -function ITensors.op( - ::OpName"Rz+z+", ::SiteType"S=1/2"; θ::Number - ) - a = exp(-im * θ * 0.5) - mat = zeros(ComplexF64, 4, 4) - mat[1, 1] = 1 - mat[2, 2] = 1 - mat[3, 3] = 1 - mat[4, 4] = a - return mat -end - ITensors.op(o::OpName"Rxxyy", ::SiteType"Qubit"; θ::Number) = ITensors.op(o, ITensorMPS.SiteType("S=1/2"); θ) ITensors.op(o::OpName"Rxxyyzz", ::SiteType"Qubit"; θ::Number) = ITensors.op(o, ITensorMPS.SiteType("S=1/2"); θ) From 9e36d281050fdd485c995e08c58adeda69be8f16 Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 28 Oct 2025 10:44:27 -0400 Subject: [PATCH 07/12] Formatting, contributors etc --- .github/workflows/CI.yml | 1 - Project.toml | 2 +- docs/make.jl | 22 ++++++++-------- examples/heavyhex_isingdynamics.jl | 2 +- src/constructors.jl | 12 ++++----- src/inner.jl | 2 +- src/norm_sqr.jl | 2 +- src/siteinds.jl | 4 +-- src/tensornetworkstate.jl | 4 +-- test/loopcorrections.jl | 40 ------------------------------ test/runtests.jl | 2 +- test/test_examples.jl | 13 ++++++++++ 12 files changed, 39 insertions(+), 67 deletions(-) delete mode 100644 test/loopcorrections.jl create mode 100644 test/test_examples.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 399c3c9..19fb3ad 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -24,7 +24,6 @@ jobs: matrix: version: - '1.12' - - '1.6' - 'pre' os: - ubuntu-latest diff --git a/Project.toml b/Project.toml index 08cc1eb..aad56cd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "TensorNetworkQuantumSimulator" uuid = "4de3b72a-362e-43dd-83ff-3f381eda9f9c" -authors = ["JoeyT1994 ", "MSRudolph "] +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" diff --git a/docs/make.jl b/docs/make.jl index b563f57..8ab972b 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,23 +1,23 @@ using TensorNetworkQuantumSimulator using Documenter -DocMeta.setdocmeta!(TensorNetworkQuantumSimulator, :DocTestSetup, :(using TensorNetworkQuantumSimulator); recursive=true) +DocMeta.setdocmeta!(TensorNetworkQuantumSimulator, :DocTestSetup, :(using TensorNetworkQuantumSimulator); recursive = true) makedocs(; - modules=[TensorNetworkQuantumSimulator], - authors="Xuanzhao Gao and contributors", - sitename="TensorNetworkQuantumSimulator.jl", - format=Documenter.HTML(; - canonical="https://JoeyT1994.github.io/TensorNetworkQuantumSimulator.jl", - edit_link="main", - assets=String[], + modules = [TensorNetworkQuantumSimulator], + authors = ["JoeyT1994 ", "MSRudolph ", "Xuanzhao Gao ", "and contributors"], + sitename = "TensorNetworkQuantumSimulator.jl", + format = Documenter.HTML(; + canonical = "https://JoeyT1994.github.io/TensorNetworkQuantumSimulator.jl", + edit_link = "main", + assets = String[], ), - pages=[ + pages = [ "Home" => "index.md", ], ) deploydocs(; - repo="github.com/JoeyT1994/TensorNetworkQuantumSimulator.jl", - devbranch="main", + repo = "github.com/JoeyT1994/TensorNetworkQuantumSimulator.jl", + devbranch = "main", ) diff --git a/examples/heavyhex_isingdynamics.jl b/examples/heavyhex_isingdynamics.jl index b6ef82d..e5832ab 100644 --- a/examples/heavyhex_isingdynamics.jl +++ b/examples/heavyhex_isingdynamics.jl @@ -65,7 +65,7 @@ function main() println("Boundary MPS measured magnetisation on central site with MPS rank $(mps_bond_dimension) MPSs is $(only(sz_bmps))") #Sample from q(x) and get p(x) / q(x) for each sample too - nsamples = 250 + nsamples = 50 bitstrings = TN.sample_directly_certified(ψ, nsamples; norm_mps_bond_dimension = mps_bond_dimension) st_dev = Statistics.std(first.(bitstrings)) diff --git a/src/constructors.jl b/src/constructors.jl index af45e5d..1d7d864 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -5,31 +5,31 @@ const stringtostatemap = Dict("I" => [1, 0, 0, 0], "X" => [0, 1, 0, 0], "Y" => [ Tensor network for vacuum state on given graph, i.e all spins up """ -function zerostate(eltype, g::NamedGraph, s::Dictionary = siteinds(g, "S=1/2")) +function zerostate(eltype, g::NamedGraph, s::Dictionary = siteinds("S=1/2", g)) return tensornetworkstate(eltype, v -> "↑", g, s) end -zerostate(g::NamedGraph, s::Dictionary = siteinds(g, "S=1/2")) = zerostate(Float64, g, s) +zerostate(g::NamedGraph, s::Dictionary = siteinds("S=1/2", g)) = zerostate(Float64, g, s) """ topaulitensornetwork(op, g::NamedGraph) Tensor network (in Heisenberg picture). Function should map vertices of the graph to pauli strings. """ -function paulitensornetworkstate(eltype, f::Function, g::NamedGraph, s::Dictionary = siteinds(g, "Pauli")) +function paulitensornetworkstate(eltype, f::Function, g::NamedGraph, s::Dictionary = siteinds("Pauli", g)) h = v -> stringtostatemap[f(v)] return tensornetworkstate(eltype, h, g, s) end -topaulitensornetwork(f::Function, g::NamedGraph, s::Dictionary = siteinds(g, "Pauli")) = topaulitensornetwork(Float64, f, g, s) +topaulitensornetwork(f::Function, g::NamedGraph, s::Dictionary = siteinds("Pauli", g)) = topaulitensornetwork(Float64, f, g, s) """ identitytensornetwork(tninds::IndsNetwork) Tensor network (in Heisenberg picture) for identity matrix on given IndsNetwork """ -function identitytensornetworkstate(eltype, g::NamedGraph, s::Dictionary = siteinds(g, "Pauli")) +function identitytensornetworkstate(eltype, g::NamedGraph, s::Dictionary = siteinds("Pauli", g)) return paulitensornetworkstate(eltype, v -> "I", g, s) end -identitytensornetworkstate(g::NamedGraph, s::Dictionary = siteinds(g, "Pauli")) = identitytensornetworkstate(Float64, g, s) +identitytensornetworkstate(g::NamedGraph, s::Dictionary = siteinds("Pauli", g)) = identitytensornetworkstate(Float64, g, s) diff --git a/src/inner.jl b/src/inner.jl index 71d253e..3170c0d 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -35,7 +35,7 @@ end # Example ```julia - s = siteinds(g, "S=1/2") + s = siteinds("S=1/2", g) ψ = random_tensornetworkstate(ComplexF32, g, s; bond_dimension = 4) ϕ = random_tensornetworkstate(ComplexF32, g, s; bond_dimension = 4) diff --git a/src/norm_sqr.jl b/src/norm_sqr.jl index 3cfc37a..e9848a0 100644 --- a/src/norm_sqr.jl +++ b/src/norm_sqr.jl @@ -28,7 +28,7 @@ end - The computed squared norm as a scalar value. # Example ```julia - s = siteinds(g, "S=1/2") + s = siteinds("S=1/2", g) ψ = random_tensornetworkstate(ComplexF32, g, s; bond_dimension = 4) # Exact norm norm_exact = LinearAlgebra.norm(ψ; alg = "exact") diff --git a/src/siteinds.jl b/src/siteinds.jl index 1cd4368..d37ee03 100644 --- a/src/siteinds.jl +++ b/src/siteinds.jl @@ -1,10 +1,10 @@ using Dictionaries: Dictionary function default_siteinds(g::AbstractGraph) - return siteinds(g, "S=1/2") + return siteinds("S=1/2", g) end -function siteinds(g::AbstractGraph, sitetype::String, sitedimension::Int = site_dimension(sitetype)) +function siteinds(sitetype::String, g::AbstractGraph, sitedimension::Int = site_dimension(sitetype)) vs = collect(vertices(g)) return Dictionary(vs, [Index[Index(sitedimension, site_tag(sitetype))] for v in vs]) end diff --git a/src/tensornetworkstate.jl b/src/tensornetworkstate.jl index 492551a..f9391b5 100644 --- a/src/tensornetworkstate.jl +++ b/src/tensornetworkstate.jl @@ -138,7 +138,7 @@ end - A `TensorNetworkState` object representing the random tensor network state. """ function random_tensornetworkstate(eltype, g::AbstractGraph, sitetype::String, d::Int = site_dimension(sitetype); bond_dimension::Int = 1) - return random_tensornetworkstate(eltype, g, siteinds(g, sitetype, d); bond_dimension) + return random_tensornetworkstate(eltype, g, siteinds(sitetype, g, d); bond_dimension) end """ @@ -192,7 +192,7 @@ end - A `TensorNetworkState` object representing the constructed tensor network state. """ function tensornetworkstate(eltype, f::Function, g::AbstractGraph, sitetype::String, d::Int = site_dimension(sitetype)) - return tensornetworkstate(eltype, f, g, siteinds(g, sitetype, d)) + return tensornetworkstate(eltype, f, g, siteinds(sitetype, g, d)) end function random_tensornetworkstate(g::AbstractGraph, args...; kwargs...) diff --git a/test/loopcorrections.jl b/test/loopcorrections.jl deleted file mode 100644 index ae5f9a7..0000000 --- a/test/loopcorrections.jl +++ /dev/null @@ -1,40 +0,0 @@ -using TensorNetworkQuantumSimulator -import TensorNetworkQuantumSimulator as TN - -using ITensors - -using NamedGraphs -using Graphs -const NG = NamedGraphs -const G = Graphs -using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_graph - -using LinearAlgebra: norm - -using EinExprs: Greedy - -using Random -Random.seed!(1634) - -@testset "Loop Corrections" begin - nx, ny = 4, 4 - χ = 2 - ITensors.disable_warn_order() - gs = [ - (named_grid((nx, 1)), "line", 0), - (named_hexagonal_lattice_graph(nx, ny), "hexagonal", 6), - (named_grid((nx, ny)), "square", 4), - ] - for (g, g_str, smallest_loop_size) in gs - ψ = TN.random_tensornetworkstate(ComplexF32, g, "S=1/2"; bond_dimension = χ) - - ψ = normalize(ψ; alg = "bp") - - norm_bp = norm(ψ; alg = "bp") - norm_loopcorrected = norm(ψ; alg = "loopcorrections", max_configuration_size = 2 * (smallest_loop_size) - 1) - norm_exact = norm(ψ; alg = "exact") - - @test isapprox(norm_bp, norm_exact, atol = 5.0e-2) - @test isapprox(norm_loopcorrected, norm_exact, atol = 5.0e-2) - end -end diff --git a/test/runtests.jl b/test/runtests.jl index c0d01c3..c936e2e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,5 +2,5 @@ using TensorNetworkQuantumSimulator using Test @testset "TensorNetworkQuantumSimulator.jl" begin - include("loopcorrections.jl") + include("test_examples.jl") end diff --git a/test/test_examples.jl b/test/test_examples.jl new file mode 100644 index 0000000..d2e9901 --- /dev/null +++ b/test/test_examples.jl @@ -0,0 +1,13 @@ +@eval module $(gensym()) +using TensorNetworkQuantumSimulator: TensorNetworkQuantumSimulator +using Test: @testset + +@testset "Test examples" begin + example_files = [ + "boundarymps.jl", "heavyhex_isingdynamics.jl", "loopcorrections.jl", "TEBD.jl", "time_evolution.jl", "time_evolution_Heisenberg.jl", + ] + @testset "Test $example_file" for example_file in example_files + include(joinpath(pkgdir(TensorNetworkQuantumSimulator), "examples", example_file)) + end +end +end From f886fb5bc01b8e261b764545a1c880b69bebb178 Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 28 Oct 2025 12:32:26 -0400 Subject: [PATCH 08/12] x64 only. Int64 dropped --- .github/workflows/CI.yml | 1 - src/Backend/loopcorrection.jl | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 19fb3ad..ac2be0f 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -29,7 +29,6 @@ jobs: - ubuntu-latest arch: - x64 - - x86 steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 diff --git a/src/Backend/loopcorrection.jl b/src/Backend/loopcorrection.jl index 41bb874..0fa3a2e 100644 --- a/src/Backend/loopcorrection.jl +++ b/src/Backend/loopcorrection.jl @@ -3,7 +3,7 @@ using ITensorNetworks: underlying_graph function loopcorrected_partitionfunction( bp_cache::BeliefPropagationCache, - max_configuration_size::Int64, + max_configuration_size::Int, ) zbp = partitionfunction(bp_cache) bp_cache = rescale(bp_cache) From 62a613f2b196b6813b594b243352c238b565a0ea Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 28 Oct 2025 12:41:42 -0400 Subject: [PATCH 09/12] Add compat to Docs/Project.toml --- docs/Project.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index fe9f716..d0aac57 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,7 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" TensorNetworkQuantumSimulator = "fb4b0154-8226-4cf3-a39c-47b51461fc74" + +[compat] +julia = "1.12" +Documenter = "1" \ No newline at end of file From 3fadf2ad526ab5576f2f217c3ca6037c7063d979 Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 28 Oct 2025 12:58:15 -0400 Subject: [PATCH 10/12] Int -> Integer --- src/Backend/beliefpropagationcache.jl | 2 +- src/Backend/boundarympscache.jl | 6 +- src/Backend/loopcorrection.jl | 2 +- src/expect.jl | 2 +- src/inner.jl | 4 +- src/norm_sqr.jl | 4 +- src/sample.jl | 247 -------------------------- src/sampling.jl | 38 ++-- src/siteinds.jl | 2 +- src/tensornetworkstate.jl | 20 +-- 10 files changed, 40 insertions(+), 287 deletions(-) delete mode 100644 src/sample.jl diff --git a/src/Backend/beliefpropagationcache.jl b/src/Backend/beliefpropagationcache.jl index 3a60554..e731a0d 100644 --- a/src/Backend/beliefpropagationcache.jl +++ b/src/Backend/beliefpropagationcache.jl @@ -181,7 +181,7 @@ function loop_correlation(bpc::BeliefPropagationCache, loop::Vector{<:NamedEdge} end #Calculate the correlations flowing around each of the primitive loops of the BP cache -function loop_correlations(bpc::BeliefPropagationCache, smallest_loop_size::Int; kwargs...) +function loop_correlations(bpc::BeliefPropagationCache, smallest_loop_size::Integer; kwargs...) g = underlying_graph(bpc) cycles = NamedGraphs.cycle_to_path.(NamedGraphs.unique_simplecycles_limited_length(g, smallest_loop_size)) corrs = [] diff --git a/src/Backend/boundarympscache.jl b/src/Backend/boundarympscache.jl index 53e53a4..29e2a76 100644 --- a/src/Backend/boundarympscache.jl +++ b/src/Backend/boundarympscache.jl @@ -8,7 +8,7 @@ struct BoundaryMPSCache{V, N <: AbstractDataGraph{V}, M <: Union{ITensor, Vector messages::Dictionary{NamedEdge, M} supergraph::PartitionedGraph sorted_edges::Dictionary{PartitionEdge, Vector{NamedEdge}} - mps_bond_dimension::Int + mps_bond_dimension::Integer end default_update_alg(bmps_cache::BoundaryMPSCache) = "bp" @@ -137,7 +137,7 @@ end function BoundaryMPSCache( tn::Union{TensorNetworkState, ITensorNetwork, BilinearForm, QuadraticForm}, - mps_bond_dimension::Int; + mps_bond_dimension::Integer; partition_by = "row", gauge_state = true ) @@ -469,7 +469,7 @@ function update_message!( alg::Algorithm"ITensorMPS", bmps_cache::BoundaryMPSCache, pe::PartitionEdge; - maxdim::Int = mps_bond_dimension(bmps_cache), + maxdim::Integer = mps_bond_dimension(bmps_cache), ) prev_pe = prev_partitionedge(bmps_cache, pe) O = ITensorMPS.MPO(bmps_cache, src(pe)) diff --git a/src/Backend/loopcorrection.jl b/src/Backend/loopcorrection.jl index 0fa3a2e..5ddb744 100644 --- a/src/Backend/loopcorrection.jl +++ b/src/Backend/loopcorrection.jl @@ -3,7 +3,7 @@ using ITensorNetworks: underlying_graph function loopcorrected_partitionfunction( bp_cache::BeliefPropagationCache, - max_configuration_size::Int, + max_configuration_size::Integer, ) zbp = partitionfunction(bp_cache) bp_cache = rescale(bp_cache) diff --git a/src/expect.jl b/src/expect.jl index e701f5a..bf99dd6 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -138,7 +138,7 @@ function expect( observable::Union{Tuple, Vector{<:Tuple}}; cache_update_kwargs = default_bmps_update_kwargs(ψ), partition_by = boundarymps_partitioning(observable), - mps_bond_dimension::Int, + mps_bond_dimension::Integer, kwargs..., ) diff --git a/src/inner.jl b/src/inner.jl index 3170c0d..4f86699 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -23,7 +23,7 @@ end - `"boundarymps"`: Boundary MPS approximation (requires `mps_bond_dimension`). - `"loopcorrections"`: Loop corrections to belief propagation. - Extra kwargs for `alg = "boundarymps"`: - - `mps_bond_dimension::Int`: The bond dimension for the boundary MPS approximation. + - `mps_bond_dimension::Integer`: The bond dimension for the boundary MPS approximation. - `partition_by`: How to partition the graph for boundary MPS (default is `"row"`). - `cache_update_kwargs`: Additional keyword arguments for updating the cache. - Extra kwargs for `alg = "bp"` or `"loopcorrections"`: @@ -84,7 +84,7 @@ function ITensors.inner(alg::Union{Algorithm"bp", Algorithm"loopcorrections"}, return inner(alg, ψϕ_bpc; kwargs...) end -function ITensors.inner(alg::Algorithm"boundarymps", ψ::TensorNetworkState, ϕ::TensorNetworkState; mps_bond_dimension::Int, partition_by = "row", cache_update_kwargs = (;), kwargs...) +function ITensors.inner(alg::Algorithm"boundarymps", ψ::TensorNetworkState, ϕ::TensorNetworkState; mps_bond_dimension::Integer, partition_by = "row", cache_update_kwargs = (;), kwargs...) ψϕ_bmps = BoundaryMPSCache(BilinearForm(ψ, ϕ), mps_bond_dimension; partition_by) maxiter = get(cache_update_kwargs, :maxiter, default_bp_maxiter(ψϕ_bmps)) cache_update_kwargs = (; cache_update_kwargs..., maxiter) diff --git a/src/norm_sqr.jl b/src/norm_sqr.jl index e9848a0..012c00a 100644 --- a/src/norm_sqr.jl +++ b/src/norm_sqr.jl @@ -18,7 +18,7 @@ end - `"loopcorrections"`: Loop corrections to belief propagation (requires `max_configuration_size`). # Keyword Arguments - For `alg = "boundarymps"`: - - `mps_bond_dimension::Int`: The bond dimension for the boundary MPS approximation. + - `mps_bond_dimension::Integer`: The bond dimension for the boundary MPS approximation. - `partition_by`: How to partition the graph for boundary MPS (default is `"row"`). - `cache_update_kwargs`: Additional keyword arguments for updating the cache. - For `alg = "bp"` or `"loopcorrections"`: @@ -75,7 +75,7 @@ function norm_sqr(alg::Union{Algorithm"bp", Algorithm"loopcorrections"}, ψ::Ten return norm_sqr(alg, ψ_bpc; kwargs...) end -function norm_sqr(alg::Algorithm"boundarymps", ψ::TensorNetworkState; mps_bond_dimension::Int, partition_by = "row", cache_update_kwargs = default_bmps_update_kwargs(ψ), kwargs...) +function norm_sqr(alg::Algorithm"boundarymps", ψ::TensorNetworkState; mps_bond_dimension::Integer, partition_by = "row", cache_update_kwargs = default_bmps_update_kwargs(ψ), kwargs...) ψ_bmps = BoundaryMPSCache(ψ, mps_bond_dimension; partition_by) maxiter = get(cache_update_kwargs, :maxiter, default_bp_maxiter(ψ_bmps)) cache_update_kwargs = (; cache_update_kwargs..., maxiter) diff --git a/src/sample.jl b/src/sample.jl deleted file mode 100644 index 3a53caf..0000000 --- a/src/sample.jl +++ /dev/null @@ -1,247 +0,0 @@ -using StatsBase - -function _sample( - ψ::ITensorNetwork, - nsamples::Int64; - projected_message_rank::Int64, - norm_message_rank::Int64, - norm_message_update_kwargs = (; niters = _default_boundarymps_update_niters, tolerance = _default_boundarymps_update_tolerance), - projected_message_update_kwargs = (; cutoff = _default_boundarymps_update_cutoff, maxdim = projected_message_rank), - partition_by = "Row", - kwargs..., - ) - - grouping_function = partition_by == "Column" ? v -> last(v) : v -> first(v) - group_sorting_function = partition_by == "Column" ? v -> first(v) : v -> last(v) - ψ, ψψ = symmetric_gauge(ψ) - ψ, ψψ = normalize(ψ, ψψ) - - norm_MPScache = BoundaryMPSCache(ψψ; message_rank = norm_message_rank, grouping_function, group_sorting_function) - sorted_partitions = sort(ITensorNetworks.partitions(norm_MPScache)) - seq = [ - sorted_partitions[i] => sorted_partitions[i - 1] for - i in length(sorted_partitions):-1:2 - ] - norm_message_update_kwargs = (; norm_message_update_kwargs..., normalize = false) - norm_MPScache = ITensorNetworks.update_iteration(Algorithm("bp"; message_update_alg = Algorithm("orthogonal"; norm_message_update_kwargs...)), norm_MPScache, seq) - - projected_MPScache = BoundaryMPSCache(ψ; message_rank = projected_message_rank, grouping_function, group_sorting_function) - - #Generate the bit_strings moving left to right through the network - probs_and_bitstrings = NamedTuple[] - for j in 1:nsamples - p_over_q_approx, logq, bitstring = _get_one_sample( - norm_MPScache, projected_MPScache, sorted_partitions; projected_message_update_kwargs, kwargs... - ) - push!(probs_and_bitstrings, (poverq = p_over_q_approx, logq = logq, bitstring = bitstring)) - end - - return probs_and_bitstrings, ψ -end - -""" - sample( - ψ::ITensorNetwork, - nsamples::Int64; - projected_message_rank::Int64, - norm_message_rank::Int64, - norm_message_update_kwargs=(; niters = _default_boundarymps_update_niters, tolerance = _default_boundarymps_update_tolerance), - projected_message_update_kwargs = (;cutoff = _default_boundarymps_update_cutoff, maxdim = projected_message_rank), - partition_by = "Row", - kwargs..., - ) - -Take nsamples bitstrings from a 2D open boundary tensornetwork by partitioning it and using boundary MPS algorithm with relevant ranks -""" -function sample(ψ::ITensorNetwork, nsamples::Int64; kwargs...) - probs_and_bitstrings, _ = _sample(ψ::ITensorNetwork, nsamples::Int64; kwargs...) - # returns just the bitstrings - return getindex.(probs_and_bitstrings, :bitstring) -end - -""" - sample_directly_certified( - ψ::ITensorNetwork, - nsamples::Int64; - projected_message_rank::Int64, - norm_message_rank::Int64, - norm_message_update_kwargs=(; niters = _default_boundarymps_update_niters, tolerance = _default_boundarymps_update_tolerance), - projected_message_update_kwargs = (;cutoff = _default_boundarymps_update_cutoff, maxdim = projected_message_rank), - partition_by = "Row", - kwargs..., - ) - -Take nsamples bitstrings from a 2D open boundary tensornetwork by partitioning it and using boundary MPS algorithm with relevant ranks. -Returns a vector of (p/q, logq, bitstring) where loqq is log probability of drawing the bitstring and p/q attests to the quality of the bitstring which is accurate only if the projected boundary MPS rank is high enough. -""" -function sample_directly_certified(ψ::ITensorNetwork, nsamples::Int64; projected_message_rank = 5 * maxlinkdim(ψ), kwargs...) - probs_and_bitstrings, _ = _sample(ψ::ITensorNetwork, nsamples::Int64; projected_message_rank, kwargs...) - # returns the self-certified p/q, logq and bitstrings - return probs_and_bitstrings -end - -""" - sample_certified( - ψ::ITensorNetwork, - nsamples::Int64; - projected_message_rank::Int64, - norm_message_rank::Int64, - norm_message_update_kwargs=(; niters = _default_boundarymps_update_niters, tolerance = _default_boundarymps_update_tolerance), - projected_message_update_kwargs = (;cutoff = _default_boundarymps_update_cutoff, maxdim = projected_message_rank), - partition_by = "Row", - kwargs..., - ) - -Take nsamples bitstrings from a 2D open boundary tensornetwork by partitioning it and using boundary MPS algorithm with relevant ranks. For each sample perform -an independent contraction of to get a measure of p/q. -Returns a vector of (p/q, bitstring) where p/q attests to the quality of the bitstring which is accurate only if the certification boundary MPS rank is high enough. -""" -function sample_certified(ψ::ITensorNetwork, nsamples::Int64; certification_message_rank = 5 * maxlinkdim(ψ), certification_message_update_kwargs = (; cutoff = _default_boundarymps_update_cutoff), kwargs...) - probs_and_bitstrings, ψ = _sample(ψ::ITensorNetwork, nsamples::Int64; kwargs...) - # send the bitstrings and the logq to the certification function - return certify_samples(ψ, probs_and_bitstrings; certification_message_rank, certification_message_update_kwargs, symmetrize_and_normalize = false) -end - -function _get_one_sample( - norm_MPScache::BoundaryMPSCache, - projected_MPScache::BoundaryMPSCache, - sorted_partitions; - projected_message_update_kwargs = (; cutoff = _default_boundarymps_update_cutoff, maxdim = maximum_virtual_dimension(projected_MPScache)), - kwargs..., - ) - - projected_message_update_kwargs = (; projected_message_update_kwargs..., normalize = false) - - norm_MPScache = copy(norm_MPScache) - - bit_string = Dictionary{keytype(vertices(projected_MPScache)), Int}() - p_over_q_approx = nothing - logq = 0 - for (i, partition) in enumerate(sorted_partitions) - - p_over_q_approx, _logq, bit_string, = - sample_partition!(norm_MPScache, partition, bit_string; kwargs...) - vs = planargraph_vertices(norm_MPScache, partition) - logq += _logq - - projected_MPScache = update_factors( - projected_MPScache, - Dict(zip(vs, [copy(norm_MPScache[(v, "ket")]) for v in vs])), - ) - - - if i < length(sorted_partitions) - next_partition = sorted_partitions[i + 1] - - #Alternate fitting procedure here which is faster for small bond dimensions but slower for large - projected_MPScache = ITensorNetworks.update_message( - ITensorNetworks.set_default_kwargs(Algorithm("ITensorMPS"; projected_message_update_kwargs...)), - projected_MPScache, - partition => next_partition - ) - - pes = planargraph_sorted_partitionedges(norm_MPScache, partition => next_partition) - - for pe in pes - m = only(message(projected_MPScache, pe)) - set_message!(norm_MPScache, pe, [m, dag(prime(m))]) - end - end - - i > 1 && delete_partitionpair_messages!(projected_MPScache, sorted_partitions[i - 1] => sorted_partitions[i]) - i > 2 && delete_partitionpair_messages!(norm_MPScache, sorted_partitions[i - 2] => sorted_partitions[i - 1]) - end - - return p_over_q_approx, logq, bit_string -end - -function certify_sample( - ψ::ITensorNetwork, bitstring, logq::Number; - certification_message_rank::Int64, - certification_message_update_kwargs = (; maxdim = certification_message_rank, cutoff = _default_boundarymps_update_cutoff), - symmetrize_and_normalize = true, - ) - if symmetrize_and_normalize - ψ, ψψ = symmetric_gauge(ψ) - ψ = normalize(ψ, cache! = Ref(ψψ)) - end - - certification_message_update_kwargs = (; certification_message_update_kwargs..., normalize = false) - - ψproj = copy(ψ) - s = siteinds(ψ) - qv = sqrt(exp(inv(oftype(logq, length(vertices(ψ)))) * logq)) - for v in vertices(ψ) - P = adapt(datatype(ψproj[v]))(onehot(only(s[v]) => bitstring[v] + 1)) - ψproj[v] = ψproj[v] * P * inv(qv) - end - - bmpsc = BoundaryMPSCache(ψproj; message_rank = certification_message_rank) - certification_message_update_kwargs = (; normalize = false, certification_message_update_kwargs...) - - #This block is two times faster than the two lines below but likely less accurate for smaller maxdims - # pg = partitions_graph(ppg(bmpsc)) - # partition = first(center(pg)) - # seq = [src(e) => dst(e) for e in post_order_dfs_edges(pg, partition)] - #bmpsc = ITensorNetworks.update_iteration(Algorithm("bp"; message_update_alg = Algorithm("ITensorMPS"; certification_message_update_kwargs...)), bmpsc, seq) - #p_over_q = region_scalar(bmpsc, partition) - - bmpsc = ITensorNetworks.update(bmpsc, message_update_alg = Algorithm("ITensorMPS"; certification_message_update_kwargs...)) - p_over_q = scalar(bmpsc) - - - p_over_q *= conj(p_over_q) - - return (poverq = p_over_q, bitstring = bitstring) -end - -certify_sample(ψ, logq_and_bitstring::NamedTuple; kwargs...) = certify_sample(ψ, logq_and_bitstring.bitstring, logq_and_bitstring.logq; kwargs...) - -function certify_samples(ψ::ITensorNetwork, bitstrings, logqs::Vector{<:Number}; kwargs...) - return [certify_sample(ψ, bitstring, logq; kwargs...) for (bitstring, logq) in zip(bitstrings, logqs)] -end - -function certify_samples(ψ::ITensorNetwork, probs_and_bitstrings::Vector{<:NamedTuple}; kwargs...) - return [certify_sample(ψ, prob_and_bitstring; kwargs...) for prob_and_bitstring in probs_and_bitstrings] -end - -#Sample along the column/ row specified by pv with the left incoming MPS message input and the right extractable from the cache -function sample_partition!( - ψIψ::BoundaryMPSCache, - partition, - bit_string::Dictionary; - kwargs..., - ) - vs = sort(planargraph_vertices(ψIψ, partition)) - seq = PartitionEdge[PartitionEdge(vs[i] => vs[i - 1]) for i in length(vs):-1:2] - !isempty(seq) && partition_update!(ψIψ, seq) - prev_v, traces = nothing, [] - logq = 0 - for v in vs - !isnothing(prev_v) && partition_update!(ψIψ, [PartitionEdge(prev_v => v)]) - env = environment(bp_cache(ψIψ), [(v, "operator")]) - seq = contraction_sequence(env; alg = "optimal") - ρ = contract(env; sequence = seq) - ρ_tr = tr(ρ) - push!(traces, ρ_tr) - ρ *= inv(ρ_tr) - ρ_diag = collect(real.(diag(ITensors.array(ρ)))) - config = StatsBase.sample(1:length(ρ_diag), Weights(ρ_diag)) - # config is 1 or 2, but we want 0 or 1 for the sample itself - set!(bit_string, v, config - 1) - s_ind = only(filter(i -> plev(i) == 0, inds(ρ))) - P = adapt(datatype(ρ))(onehot(s_ind => config)) - q = ρ_diag[config] - logq += log(q) - ψv = copy(ψIψ[(v, "ket")]) * inv(sqrt(q)) - ψv = P * ψv - setindex_preserve_graph!(ψIψ, adapt(datatype(ρ))(ITensor(one(Bool))), (v, "operator")) - setindex_preserve_graph!(ψIψ, copy(ψv), (v, "ket")) - setindex_preserve_graph!(ψIψ, dag(prime(copy(ψv))), (v, "bra")) - prev_v = v - end - - delete_partition_messages!(ψIψ, partition) - - return first(traces), logq, bit_string -end diff --git a/src/sampling.jl b/src/sampling.jl index 399d04e..e001ce6 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -3,9 +3,9 @@ using StatsBase function sample( alg::Algorithm"boundarymps", ψ::TensorNetworkState, - nsamples::Int64; - projected_mps_bond_dimension::Int, - norm_mps_bond_dimension::Int, + nsamples::Integer; + projected_mps_bond_dimension::Integer, + norm_mps_bond_dimension::Integer, norm_cache_message_update_kwargs = (;), partition_by = "row", gauge_state = true, @@ -32,9 +32,9 @@ end """ sample( ψ::ITensorNetwork, - nsamples::Int64; - projected_message_rank::Int64, - norm_message_rank::Int64, + nsamples::Integer; + projected_message_rank::Integer, + norm_message_rank::Integer, norm_message_update_kwargs=(; niters = _default_boundarymps_update_niters, tolerance = _default_boundarymps_update_tolerance), projected_message_update_kwargs = (;cutoff = _default_boundarymps_update_cutoff, maxdim = projected_message_rank), partition_by = "Row", @@ -46,7 +46,7 @@ Take nsamples bitstrings, based on the square of the coefficients of the vector Arguments --------- - `ψ::ITensorNetwork`: The tensornetwork state to sample from. -- `nsamples::Int64`: Number of samples to draw. +- `nsamples::Integer`: Number of samples to draw. Keyword Arguments ----------------- @@ -62,7 +62,7 @@ Returns ------- A vector of bitstrings sampled from the probability distribution defined by as a dictionary mapping each vertex to a configuration (0...d). """ -function sample(ψ::TensorNetworkState, nsamples::Int64; alg = "boundarymps", kwargs...) +function sample(ψ::TensorNetworkState, nsamples::Integer; alg = "boundarymps", kwargs...) algorithm_check(ψ, "sample", alg) probs_and_bitstrings, _ = sample(Algorithm(alg), ψ, nsamples; kwargs...) # returns just the bitstrings @@ -72,9 +72,9 @@ end """ sample_directly_certified( ψ::ITensorNetwork, - nsamples::Int64; - projected_message_rank::Int64, - norm_message_rank::Int64, + nsamples::Integer; + projected_message_rank::Integer, + norm_message_rank::Integer, norm_message_update_kwargs=(; niters = _default_boundarymps_update_niters, tolerance = _default_boundarymps_update_tolerance), projected_message_update_kwargs = (;cutoff = _default_boundarymps_update_cutoff, maxdim = projected_message_rank), partition_by = "Row", @@ -87,7 +87,7 @@ The samples are drawn from x~q(x) and for each sample is calculated "on-t Arguments --------- - `ψ::ITensorNetwork`: The tensornetwork state to sample from. -- `nsamples::Int64`: Number of samples to draw. +- `nsamples::Integer`: Number of samples to draw. Keyword Arguments ----------------- @@ -107,7 +107,7 @@ Each NamedTuple contains: - `logq`: Log probability of drawing the bitstring. - `bitstring`: The sampled bitstring as a dictionary mapping each vertex to a configuration (0...d). """ -function sample_directly_certified(ψ::TensorNetworkState, nsamples::Int64; projected_mps_bond_dimension = 5 * maxlinkdim(ψ), alg = "boundarymps", kwargs...) +function sample_directly_certified(ψ::TensorNetworkState, nsamples::Integer; projected_mps_bond_dimension = 5 * maxlinkdim(ψ), alg = "boundarymps", kwargs...) algorithm_check(ψ, "sample", alg) probs_and_bitstrings, _ = sample(Algorithm(alg), ψ, nsamples; projected_mps_bond_dimension, kwargs...) # returns the self-certified p/q, logq and bitstrings @@ -117,9 +117,9 @@ end """ sample_certified( ψ::ITensorNetwork, - nsamples::Int64; - projected_message_rank::Int64, - norm_message_rank::Int64, + nsamples::Integer; + projected_message_rank::Integer, + norm_message_rank::Integer, norm_message_update_kwargs=(; niters = _default_boundarymps_update_niters, tolerance = _default_boundarymps_update_tolerance), projected_message_update_kwargs = (;cutoff = _default_boundarymps_update_cutoff, maxdim = projected_message_rank), partition_by = "Row", @@ -132,7 +132,7 @@ The samples are drawn from x~q(x) and for each sample an independent contraction Arguments --------- - `ψ::ITensorNetwork`: The tensornetwork state to sample from. -- `nsamples::Int64`: Number of samples to draw. +- `nsamples::Integer`: Number of samples to draw. Keyword Arguments ----------------- @@ -162,7 +162,7 @@ end function get_one_sample( norm_bmps_cache::BoundaryMPSCache, seq::Vector{<:PartitionEdge}; - projected_mps_bond_dimension::Int, + projected_mps_bond_dimension::Integer, kwargs..., ) norm_bmps_cache = copy(norm_bmps_cache) @@ -253,7 +253,7 @@ end function certify_sample( alg::Algorithm"boundarymps", ψ::TensorNetworkState, bitstring, logq::Number; - certification_mps_bond_dimension::Int, + certification_mps_bond_dimension::Integer, certification_cache_message_update_kwargs = (;), symmetrize_and_normalize = true, ) diff --git a/src/siteinds.jl b/src/siteinds.jl index d37ee03..8e68c17 100644 --- a/src/siteinds.jl +++ b/src/siteinds.jl @@ -4,7 +4,7 @@ function default_siteinds(g::AbstractGraph) return siteinds("S=1/2", g) end -function siteinds(sitetype::String, g::AbstractGraph, sitedimension::Int = site_dimension(sitetype)) +function siteinds(sitetype::String, g::AbstractGraph, sitedimension::Integer = site_dimension(sitetype)) vs = collect(vertices(g)) return Dictionary(vs, [Index[Index(sitedimension, site_tag(sitetype))] for v in vs]) end diff --git a/src/tensornetworkstate.jl b/src/tensornetworkstate.jl index f9391b5..078f04d 100644 --- a/src/tensornetworkstate.jl +++ b/src/tensornetworkstate.jl @@ -100,19 +100,19 @@ end #TODO: Default to spin 1/2 """ - random_tensornetworkstate(eltype, g::AbstractGraph, siteinds::Dictionary; bond_dimension::Int = 1) + random_tensornetworkstate(eltype, g::AbstractGraph, siteinds::Dictionary; bond_dimension::Integer = 1) Generate a random TensorNetworkState on graph `g` with local state indices given by the dictionary `siteinds`. Arguments: - `eltype`: (Optional) The number type of the tensor elements (e.g. Float64, ComplexF32). Default is Float64. - `g::AbstractGraph`: The underlying graph of the tensor network. - `siteinds::Dictionary`: A dictionary mapping vertices to ITensor indices representing the local states. Defaults to spin 1/2. - - `bond_dimension::Int`: The bond dimension of the virtual indices connecting neighboring tensors (default is 1). + - `bond_dimension::Integer`: The bond dimension of the virtual indices connecting neighboring tensors (default is 1). Returns: - A `TensorNetworkState` object representing the random tensor network state. """ -function random_tensornetworkstate(eltype, g::AbstractGraph, siteinds::Dictionary = default_siteinds(g); bond_dimension::Int = 1) +function random_tensornetworkstate(eltype, g::AbstractGraph, siteinds::Dictionary = default_siteinds(g); bond_dimension::Integer = 1) vs = collect(vertices(g)) l = Dict(e => Index(bond_dimension) for e in edges(g)) l = merge(l, Dict(reverse(e) => l[e] for e in edges(g))) @@ -125,19 +125,19 @@ function random_tensornetworkstate(eltype, g::AbstractGraph, siteinds::Dictionar end """ - random_tensornetworkstate(eltype, g::AbstractGraph, sitetype::String, d::Int = site_dimension(sitetype); bond_dimension::Int = 1) + random_tensornetworkstate(eltype, g::AbstractGraph, sitetype::String, d::Integer = site_dimension(sitetype); bond_dimension::Integer = 1) Generate a random TensorNetworkState on graph `g` with local state indices generated from the `sitetype` string (e.g. "S=1/2", "Pauli") and the local dimension `d` (default is 2 for "S=1/2", 4 for Pauli etc). Arguments: - `eltype`: (Optional) The number type of the tensor elements (e.g. Float64, ComplexF32). Default is Float64. - `g::AbstractGraph`: The underlying graph of the tensor network. - `sitetype::String`: A string representing the type of local site (e.g. "S=1/2", "Pauli"). - - `d::Int`: The local dimension of the site (default is determined by `sitetype`). - - `bond_dimension::Int`: The bond dimension of the virtual indices connecting neighboring tensors (default is 1). + - `d::Integer`: The local dimension of the site (default is determined by `sitetype`). + - `bond_dimension::Integer`: The bond dimension of the virtual indices connecting neighboring tensors (default is 1). Returns: - A `TensorNetworkState` object representing the random tensor network state. """ -function random_tensornetworkstate(eltype, g::AbstractGraph, sitetype::String, d::Int = site_dimension(sitetype); bond_dimension::Int = 1) +function random_tensornetworkstate(eltype, g::AbstractGraph, sitetype::String, d::Integer = site_dimension(sitetype); bond_dimension::Integer = 1) return random_tensornetworkstate(eltype, g, siteinds(sitetype, g, d); bond_dimension) end @@ -177,7 +177,7 @@ function tensornetworkstate(eltype, f::Function, g::AbstractGraph, siteinds::Dic end """ - tensornetworkstate(eltype, f::Function, g::AbstractGraph, sitetype::String, d::Int = site_dimension(sitetype)) + tensornetworkstate(eltype, f::Function, g::AbstractGraph, sitetype::String, d::Integer = site_dimension(sitetype)) Construct a TensorNetworkState on graph `g` where the function `f` maps vertices to local states. The local states can be given as strings (e.g. "↑", "↓", "0", "1", "I", "X", "Y", "Z") or as vectors of numbers (e.g. [1,0], [0,1], [1/sqrt(2), 1/sqrt(2)]). @@ -186,12 +186,12 @@ end - `f::Function`: A function mapping vertices of the graph to local states. - `g::AbstractGraph`: The underlying graph of the tensor network. - `sitetype::String`: A string representing the type of local site (e.g. "S=1/2", "Pauli"). - - `d::Int`: The local dimension of the site (default is determined by `sitetype`). + - `d::Integer`: The local dimension of the site (default is determined by `sitetype`). Returns: - A `TensorNetworkState` object representing the constructed tensor network state. """ -function tensornetworkstate(eltype, f::Function, g::AbstractGraph, sitetype::String, d::Int = site_dimension(sitetype)) +function tensornetworkstate(eltype, f::Function, g::AbstractGraph, sitetype::String, d::Integer = site_dimension(sitetype)) return tensornetworkstate(eltype, f, g, siteinds(sitetype, g, d)) end From 919717266a4e4c6ecf738f207b06802cc1b00723 Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 28 Oct 2025 13:16:21 -0400 Subject: [PATCH 11/12] Try to fix Docs Project.toml --- docs/Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index d0aac57..5fa75da 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,5 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -TensorNetworkQuantumSimulator = "fb4b0154-8226-4cf3-a39c-47b51461fc74" [compat] julia = "1.12" From 1346439acde0374cad042fcd268b9aeaf1d3027c Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 28 Oct 2025 13:26:15 -0400 Subject: [PATCH 12/12] Try to fix Docs Project.toml --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 8ab972b..6cd8d76 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -5,7 +5,7 @@ DocMeta.setdocmeta!(TensorNetworkQuantumSimulator, :DocTestSetup, :(using Tensor makedocs(; modules = [TensorNetworkQuantumSimulator], - authors = ["JoeyT1994 ", "MSRudolph ", "Xuanzhao Gao ", "and contributors"], + authors = "JoeyT1994 , MSRudolph , Xuanzhao Gao , and contributors", sitename = "TensorNetworkQuantumSimulator.jl", format = Documenter.HTML(; canonical = "https://JoeyT1994.github.io/TensorNetworkQuantumSimulator.jl",