From 62ad944f29178ed46136ead051ae80b1df1e7a14 Mon Sep 17 00:00:00 2001 From: Joey Date: Thu, 21 Mar 2024 11:04:07 +0000 Subject: [PATCH 01/35] Refactor inner interface. Add logscalar(tn) functionality --- src/ITensorNetworks.jl | 2 +- src/abstractitensornetwork.jl | 37 -------------------- src/contract.jl | 39 +++++++++++++++++++++ src/inner.jl | 62 +++++++++++++++++++++++++++++++++ test/test_apply.jl | 12 +++---- test/test_belief_propagation.jl | 10 +++--- 6 files changed, 112 insertions(+), 50 deletions(-) create mode 100644 src/inner.jl diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 35f5702e..4bce2297 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -129,7 +129,7 @@ include(joinpath("treetensornetworks", "solvers", "contract.jl")) include(joinpath("treetensornetworks", "solvers", "linsolve.jl")) include(joinpath("treetensornetworks", "solvers", "tree_sweeping.jl")) include("apply.jl") - +include("inner.jl") include("exports.jl") function __init__() diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index 6905694b..c8cdb3e9 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -684,43 +684,6 @@ function norm_network(tn::AbstractITensorNetwork) return tntn end -# TODO: Use or replace with `flatten_networks` -function inner_network( - tn1::AbstractITensorNetwork, - tn2::AbstractITensorNetwork; - map_bra_linkinds=sim, - combine_linkinds=false, - flatten=combine_linkinds, - kwargs..., -) - @assert issetequal(vertices(tn1), vertices(tn2)) - tn1 = map_bra_linkinds(tn1; sites=[]) - inner_net = ⊗(dag(tn1), tn2; kwargs...) - if flatten - for v in vertices(tn1) - inner_net = contract(inner_net, (v, 2) => (v, 1); merged_vertex=v) - end - end - if combine_linkinds - inner_net = ITensorNetworks.combine_linkinds(inner_net) - end - return inner_net -end - -# TODO: Rename `inner`. -function contract_inner( - ϕ::AbstractITensorNetwork, - ψ::AbstractITensorNetwork; - sequence=nothing, - contraction_sequence_kwargs=(;), -) - tn = inner_network(ϕ, ψ; combine_linkinds=true) - if isnothing(sequence) - sequence = contraction_sequence(tn; contraction_sequence_kwargs...) - end - return contract(tn; sequence)[] -end - # TODO: rename `sqnorm` to match https://github.com/JuliaStats/Distances.jl, # or `norm_sqr` to match `LinearAlgebra.norm_sqr` norm_sqr(ψ::AbstractITensorNetwork; sequence) = contract_inner(ψ, ψ; sequence) diff --git a/src/contract.jl b/src/contract.jl index 77c4ee1a..80662cb3 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -44,3 +44,42 @@ function contract_exact( end return ITensor[out] end + +scalar(tn::Union{AbstractITensorNetwork, Vector{ITensor}}; alg = "exact", kwargs...) = scalar(Algorithm(alg), tn; kwargs...) +scalar(alg::Algorithm, tn::Union{AbstractITensorNetwork, Vector{ITensor}}; kwargs...) = contract(alg, tn; kwargs...)[] + +function logscalar(tn::AbstractITensorNetwork; alg::String="exact", kwargs...) + return logscalar(Algorithm(alg), tn; kwargs...) +end + +function logscalar(alg::Algorithm"exact", tn::AbstractITensorNetwork; kwargs...) + return log(scalar(alg, tn; kwargs...)) +end + +function logscalar(alg::Algorithm"bp", tn::AbstractITensorNetwork; (cache!)=nothing, + partitioned_vertices=default_partitioned_vertices(tn), + update_cache=isnothing(cache!), + cache_update_kwargs=default_cache_update_kwargs(cache!),) + + if isnothing(cache!) + cache! = Ref(BeliefPropagationCache(tn, partitioned_vertices)) + end + + if update_cache + cache![] = update(cache![]; cache_update_kwargs...) + end + + pg = partitioned_itensornetwork(bp_cache) + + log_numerator, log_denominator = 0, 0 + for pv in partitionvertices(pg) + incoming_mts = incoming_messages(bp_cache, [pv]) + local_state = factor(bp_cache, pv) + log_numerator += logscalar(vcat(incoming_mts, local_state); alg = "exact") + end + for pe in partitionedges(pg) + log_denominator += logscalar(vcat(message(bp_cache, pe), message(bp_cache, reverse(pe))); alg = "exact") + end + + return log_numerator - log_denominator +end \ No newline at end of file diff --git a/src/inner.jl b/src/inner.jl new file mode 100644 index 00000000..70cca146 --- /dev/null +++ b/src/inner.jl @@ -0,0 +1,62 @@ +default_inner_partitioned_vertices(tn) = group(v -> first(v), vertices(tn)) + +function inner(ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; alg = "exact", kwargs...) + return inner(Algorithm(alg), ϕ, ψ; kwargs...) +end + +function inner(alg::Algorithm"exact", + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + sequence=nothing, + combine_linkinds = true, + flatten = true, + contraction_sequence_kwargs=(;), + ) + tn = flatten_networks(dag(ϕ), ψ; combine_linkinds, flatten) + if isnothing(sequence) + sequence = contraction_sequence(tn; contraction_sequence_kwargs...) + end + return scalar(tn; sequence) +end + +function inner(alg::Algorithm"exact", + A::AbstractITensorNetwork, + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + site_index_map = default_index_map, + combine_linkinds, + flatten, + kwargs..., + ) + tn = flatten_networks(site_index_map(dag(ϕ); links = []), A, ψ; combine_linkinds, flatten) + if isnothing(sequence) + sequence = contraction_sequence(tn; contraction_sequence_kwargs...) + end + return scalar(tn; sequence) +end + +loginner(ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; alg = "exact", kwargs...) = loginner(Algorithm(alg), ϕ, ψ; kwargs...) +loginner(alg::Algorithm"exact", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; kwargs...) = log(inner(alg, ϕ, ψ); kwargs...) + +function loginner(alg::Algorithm"bp", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; + partitioned_vertices = default_inner_partitioned_vertices, + kwargs...) + ϕ_dag = sim(dag(ϕ); sites = []) + tn = disjoint_union("bra" => ϕ_dag, "ket" => ψ) + return logscalar(alg, tn; partitioned_vertices, kwargs...) +end + +function loginner(alg::Algorithm"bp", A::AbstractITensorNetwork, ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; + partitioned_vertices = default_inner_partitioned_vertices, site_index_map = default_index_map, + kwargs...) + ϕ_dag = site_index_map(sim(dag(ϕ); sites = []); links = []) + tn = disjoint_union("operator" => A, "bra" => ϕ_dag, "ket" => ψ) + return logscalar(alg, tn; partitioned_vertices, kwargs...) +end + +inner(alg::Algorithm"bp", ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; kwargs...) = exp(loginner(alg, ϕ, ψ; kwargs...)) + + + diff --git a/test/test_apply.jl b/test/test_apply.jl index 0a815b9a..126aebf0 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -65,14 +65,14 @@ using SplitApplyCombine envisposdef=true, ) fSBP = - contract_inner(ψOSBP, ψOexact) / - sqrt(contract_inner(ψOexact, ψOexact) * contract_inner(ψOSBP, ψOSBP)) + inner(ψOSBP, ψOexact) / + sqrt(inner(ψOexact, ψOexact) * inner(ψOSBP, ψOSBP)) fVidal = - contract_inner(ψOVidal_symm, ψOexact) / - sqrt(contract_inner(ψOexact, ψOexact) * contract_inner(ψOVidal_symm, ψOVidal_symm)) + inner(ψOVidal_symm, ψOexact) / + sqrt(inner(ψOexact, ψOexact) * inner(ψOVidal_symm, ψOVidal_symm)) fGBP = - contract_inner(ψOGBP, ψOexact) / - sqrt(contract_inner(ψOexact, ψOexact) * contract_inner(ψOGBP, ψOGBP)) + inner(ψOGBP, ψOexact) / + sqrt(inner(ψOexact, ψOexact) * inner(ψOGBP, ψOGBP)) @test real(fGBP * conj(fGBP)) >= real(fSBP * conj(fSBP)) diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index a388e758..2085589e 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -36,7 +36,7 @@ ITensors.disable_warn_order() Oψ = copy(ψ) Oψ[v] = apply(op("Sz", s[v]), ψ[v]) - exact_sz = contract_inner(Oψ, ψ) / contract_inner(ψ, ψ) + exact_sz = inner(Oψ, ψ) / inner(ψ, ψ) bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) bpc = update(bpc) @@ -66,7 +66,7 @@ ITensors.disable_warn_order() Oψ = copy(ψ) Oψ[v] = apply(op("Sz", s[v]), ψ[v]) - exact_sz = contract_inner(Oψ, ψ) / contract_inner(ψ, ψ) + exact_sz = inner(Oψ, ψ) / inner(ψ, ψ) bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) bpc = update(bpc) @@ -86,9 +86,7 @@ ITensors.disable_warn_order() ψOψ = ising_network(s, beta; szverts=vs) contract_seq = contraction_sequence(ψψ) - actual_szsz = - ITensors.contract(ψOψ; sequence=contract_seq)[] / - ITensors.contract(ψψ; sequence=contract_seq)[] + actual_szsz = contract(ψOψ; sequence=contract_seq)[] / contract(ψψ; sequence=contract_seq)[] bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) bpc = update(bpc; maxiter=20) @@ -113,7 +111,7 @@ ITensors.disable_warn_order() ψψsplit = split_index(ψψ, NamedEdge.([(v, 1) => (v, 2) for v in vs])) env_tensors = incoming_messages(bpc, [(v, 2) for v in vs]) - rdm = ITensors.contract( + rdm = contract( vcat(env_tensors, ITensor[ψψsplit[vp] for vp in [(v, 2) for v in vs]]) ) From fe27c863e3d4d33e53747a6d9d0d6e08547f6418 Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 22 Mar 2024 15:53:14 +0000 Subject: [PATCH 02/35] Updated tests and examples to work with new inner interface --- .../contraction_sequence.jl | 14 +- src/ITensorNetworks.jl | 2 +- src/abstractitensornetwork.jl | 9 +- src/caches/beliefpropagationcache.jl | 2 +- src/contract.jl | 45 ++++-- src/expect.jl | 14 +- src/inner.jl | 147 ++++++++++++------ test/test_apply.jl | 10 +- test/test_belief_propagation.jl | 9 +- test/test_gauging.jl | 5 +- test/test_tno.jl | 14 +- 11 files changed, 170 insertions(+), 101 deletions(-) diff --git a/examples/contraction_sequence/contraction_sequence.jl b/examples/contraction_sequence/contraction_sequence.jl index 34eb097d..96ac23d4 100644 --- a/examples/contraction_sequence/contraction_sequence.jl +++ b/examples/contraction_sequence/contraction_sequence.jl @@ -14,30 +14,30 @@ s = siteinds("S=1/2", g) χ = 10 ψ = randomITensorNetwork(s; link_space=χ) -tn = norm_sqr_network(ψ) +ψψ = norm_sqr_network(ψ) # Contraction sequence for exactly computing expectation values # contract_edges = map(t -> (1, t...), collect(keys(cartesian_to_linear(system_dims)))) # inner_sequence = reduce((x, y) -> [x, y], contract_edges) println("optimal") -seq_optimal = @time contraction_sequence(tn; alg="optimal") +seq_optimal = @time contraction_sequence(ψψ; alg="optimal") using OMEinsumContractionOrders println("greedy") -seq_greedy = @time contraction_sequence(tn; alg="greedy") -res_greedy = @time contract(tn; sequence=seq_greedy) +seq_greedy = @time contraction_sequence(ψψ; alg="greedy") +res_greedy = @time contract(ψψ; sequence=seq_greedy) println("tree_sa") -seq_tree_sa = @time contraction_sequence(tn; alg="tree_sa") +seq_tree_sa = @time contraction_sequence(ψψ; alg="tree_sa") println("sa_bipartite") -seq_sa_bipartite = @time contraction_sequence(tn; alg="sa_bipartite") +seq_sa_bipartite = @time contraction_sequence(ψψ; alg="sa_bipartite") using KaHyPar println("kahypar_bipartite") seq_kahypar_bipartite = @time contraction_sequence( - tn; alg="kahypar_bipartite", sc_target=200 + ψψ; alg="kahypar_bipartite", sc_target=200 ) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 4bce2297..a1e2c8d9 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -78,7 +78,6 @@ include("opsum.jl") include("sitetype.jl") include("abstractitensornetwork.jl") include("contraction_sequences.jl") -include("expect.jl") include("models.jl") include("tebd.jl") include("itensornetwork.jl") @@ -130,6 +129,7 @@ include(joinpath("treetensornetworks", "solvers", "linsolve.jl")) include(joinpath("treetensornetworks", "solvers", "tree_sweeping.jl")) include("apply.jl") include("inner.jl") +include("expect.jl") include("exports.jl") function __init__() diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index c8cdb3e9..79768243 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -684,10 +684,11 @@ function norm_network(tn::AbstractITensorNetwork) return tntn end -# TODO: rename `sqnorm` to match https://github.com/JuliaStats/Distances.jl, -# or `norm_sqr` to match `LinearAlgebra.norm_sqr` -norm_sqr(ψ::AbstractITensorNetwork; sequence) = contract_inner(ψ, ψ; sequence) - +function inner_network( + ψ::AbstractITensorNetwork, ϕ::AbstractITensorNetwork; flatten=false, kwargs... +) + return flatten_networks(ψ, dag(ϕ); flatten, kwargs...) +end norm_sqr_network(ψ::AbstractITensorNetwork; kwargs...) = inner_network(ψ, ψ; kwargs...) # diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index f5337784..70e4a00b 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -99,7 +99,7 @@ function incoming_messages( ) bpes = boundary_partitionedges(bp_cache, partition_vertices; dir=:in) ms = messages(bp_cache, setdiff(bpes, ignore_edges)) - return reduce(vcat, ms; init=[]) + return reduce(vcat, ms; init=ITensor[]) end function incoming_messages( diff --git a/src/contract.jl b/src/contract.jl index 80662cb3..4291cdf4 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -9,6 +9,8 @@ function contract( return contract(Vector{ITensor}(tn); sequence=sequence_linear_index, kwargs...) end +contract(alg::Algorithm"exact", tn::Vector{ITensor}; kwargs...) = contract(tn; kwargs...) + function contract( alg::Union{Algorithm"density_matrix",Algorithm"ttn_svd"}, tn::AbstractITensorNetwork; @@ -45,22 +47,35 @@ function contract_exact( return ITensor[out] end -scalar(tn::Union{AbstractITensorNetwork, Vector{ITensor}}; alg = "exact", kwargs...) = scalar(Algorithm(alg), tn; kwargs...) -scalar(alg::Algorithm, tn::Union{AbstractITensorNetwork, Vector{ITensor}}; kwargs...) = contract(alg, tn; kwargs...)[] +function scalar(tn::Union{AbstractITensorNetwork,Vector{ITensor}}; alg="exact", kwargs...) + return scalar(Algorithm(alg), tn; kwargs...) +end +function scalar( + alg::Algorithm, tn::Union{AbstractITensorNetwork,Vector{ITensor}}; kwargs... +) + return contract(alg, tn; kwargs...)[] +end -function logscalar(tn::AbstractITensorNetwork; alg::String="exact", kwargs...) +function logscalar( + tn::Union{AbstractITensorNetwork,Vector{ITensor}}; alg::String="exact", kwargs... +) return logscalar(Algorithm(alg), tn; kwargs...) end -function logscalar(alg::Algorithm"exact", tn::AbstractITensorNetwork; kwargs...) - return log(scalar(alg, tn; kwargs...)) +function logscalar( + alg::Algorithm"exact", tn::Union{AbstractITensorNetwork,Vector{ITensor}}; kwargs... +) + return log(complex(scalar(alg, tn; kwargs...))) end -function logscalar(alg::Algorithm"bp", tn::AbstractITensorNetwork; (cache!)=nothing, +function logscalar( + alg::Algorithm"bp", + tn::AbstractITensorNetwork; + (cache!)=nothing, partitioned_vertices=default_partitioned_vertices(tn), update_cache=isnothing(cache!), - cache_update_kwargs=default_cache_update_kwargs(cache!),) - + cache_update_kwargs=default_cache_update_kwargs(cache!), +) if isnothing(cache!) cache! = Ref(BeliefPropagationCache(tn, partitioned_vertices)) end @@ -69,17 +84,19 @@ function logscalar(alg::Algorithm"bp", tn::AbstractITensorNetwork; (cache!)=noth cache![] = update(cache![]; cache_update_kwargs...) end - pg = partitioned_itensornetwork(bp_cache) + pg = partitioned_itensornetwork(cache![]) log_numerator, log_denominator = 0, 0 for pv in partitionvertices(pg) - incoming_mts = incoming_messages(bp_cache, [pv]) - local_state = factor(bp_cache, pv) - log_numerator += logscalar(vcat(incoming_mts, local_state); alg = "exact") + incoming_mts = incoming_messages(cache![], [pv]) + local_state = factor(cache![], pv) + log_numerator += logscalar(vcat(incoming_mts, local_state); alg="exact") end for pe in partitionedges(pg) - log_denominator += logscalar(vcat(message(bp_cache, pe), message(bp_cache, reverse(pe))); alg = "exact") + log_denominator += logscalar( + vcat(message(cache![], pe), message(cache![], reverse(pe))); alg="exact" + ) end return log_numerator - log_denominator -end \ No newline at end of file +end diff --git a/src/expect.jl b/src/expect.jl index 7feb2e11..fd939e42 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -5,6 +5,7 @@ function expect( maxdim=nothing, ortho=false, sequence=nothing, + flatten=true, vertices=vertices(ψ), ) s = siteinds(ψ) @@ -12,13 +13,13 @@ function expect( # ElT = ishermitian(ITensors.op(op, s[vertices[1]])) ? real(ElT) : ElT res = Dictionary(vertices, Vector{ElT}(undef, length(vertices))) if isnothing(sequence) - sequence = contraction_sequence(inner_network(ψ, ψ; flatten=true)) + sequence = contraction_sequence(inner_network(ψ, ψ; flatten)) end - normψ² = norm_sqr(ψ; sequence) + normψ² = norm_sqr(ψ; alg="exact", sequence, flatten) for v in vertices O = ITensor(Op(op, v), s) Oψ = apply(O, ψ; cutoff, maxdim, ortho) - res[v] = contract_inner(ψ, Oψ; sequence) / normψ² + res[v] = inner(ψ, Oψ; alg="exact", sequence, flatten) / normψ² end return res end @@ -30,16 +31,17 @@ function expect( maxdim=nothing, ortho=false, sequence=nothing, + flatten=true, ) s = siteinds(ψ) # h⃗ = Vector{ITensor}(ℋ, s) if isnothing(sequence) - sequence = contraction_sequence(inner_network(ψ, ψ; flatten=true)) + sequence = contraction_sequence(inner_network(ψ, ψ; flatten)) end h⃗ψ = [apply(hᵢ, ψ; cutoff, maxdim, ortho) for hᵢ in ITensors.terms(ℋ)] - ψhᵢψ = [contract_inner(ψ, hᵢψ; sequence) for hᵢψ in h⃗ψ] + ψhᵢψ = [inner(ψ, hᵢψ; alg="exact", flatten, sequence) for hᵢψ in h⃗ψ] ψh⃗ψ = sum(ψhᵢψ) - ψψ = norm_sqr(ψ; sequence) + ψψ = norm_sqr(ψ; alg="exact", sequence, flatten) return ψh⃗ψ / ψψ end diff --git a/src/inner.jl b/src/inner.jl index 70cca146..b1d18d3d 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -1,62 +1,117 @@ default_inner_partitioned_vertices(tn) = group(v -> first(v), vertices(tn)) +default_algorithm(tns::Vector) = all(is_tree.(tns)) ? "bp" : "exact" -function inner(ϕ::AbstractITensorNetwork, - ψ::AbstractITensorNetwork; alg = "exact", kwargs...) - return inner(Algorithm(alg), ϕ, ψ; kwargs...) +function inner( + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + alg=default_algorithm([ϕ, ψ]), + kwargs..., +) + return inner(Algorithm(alg), ϕ, ψ; kwargs...) end -function inner(alg::Algorithm"exact", - ϕ::AbstractITensorNetwork, - ψ::AbstractITensorNetwork; - sequence=nothing, - combine_linkinds = true, - flatten = true, - contraction_sequence_kwargs=(;), - ) - tn = flatten_networks(dag(ϕ), ψ; combine_linkinds, flatten) - if isnothing(sequence) - sequence = contraction_sequence(tn; contraction_sequence_kwargs...) - end - return scalar(tn; sequence) +function inner( + A::AbstractITensorNetwork, + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + alg=default_algorithm([A, ϕ, ψ]), + kwargs..., +) + return inner(Algorithm(alg), A, ϕ, ψ; kwargs...) end -function inner(alg::Algorithm"exact", - A::AbstractITensorNetwork, - ϕ::AbstractITensorNetwork, - ψ::AbstractITensorNetwork; - site_index_map = default_index_map, - combine_linkinds, - flatten, - kwargs..., - ) - tn = flatten_networks(site_index_map(dag(ϕ); links = []), A, ψ; combine_linkinds, flatten) - if isnothing(sequence) - sequence = contraction_sequence(tn; contraction_sequence_kwargs...) - end - return scalar(tn; sequence) +function inner( + alg::Algorithm"exact", + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + sequence=nothing, + contraction_sequence_kwargs=(;), + kwargs..., +) + tn = inner_network(ϕ, ψ; kwargs...) + if isnothing(sequence) + sequence = contraction_sequence(tn; contraction_sequence_kwargs...) + end + return scalar(tn; sequence) end -loginner(ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; alg = "exact", kwargs...) = loginner(Algorithm(alg), ϕ, ψ; kwargs...) -loginner(alg::Algorithm"exact", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; kwargs...) = log(inner(alg, ϕ, ψ); kwargs...) +function inner( + alg::Algorithm"exact", + A::AbstractITensorNetwork, + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + sequence=nothing, + contraction_sequence_kwargs=(;), + site_index_map=prime, + kwargs..., +) + tn = flatten_networks(site_index_map(dag(ϕ); links=[]), A, ψ; kwargs...) + if isnothing(sequence) + sequence = contraction_sequence(tn; contraction_sequence_kwargs...) + end + return scalar(tn; sequence) +end -function loginner(alg::Algorithm"bp", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - partitioned_vertices = default_inner_partitioned_vertices, - kwargs...) - ϕ_dag = sim(dag(ϕ); sites = []) - tn = disjoint_union("bra" => ϕ_dag, "ket" => ψ) - return logscalar(alg, tn; partitioned_vertices, kwargs...) +function loginner( + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + alg=default_algorithm(ϕ, ψ), + kwargs..., +) + return loginner(Algorithm(alg), ϕ, ψ; kwargs...) +end +function loginner( + alg::Algorithm"exact", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; kwargs... +) + return log(inner(alg, ϕ, ψ); kwargs...) end -function loginner(alg::Algorithm"bp", A::AbstractITensorNetwork, ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - partitioned_vertices = default_inner_partitioned_vertices, site_index_map = default_index_map, - kwargs...) - ϕ_dag = site_index_map(sim(dag(ϕ); sites = []); links = []) - tn = disjoint_union("operator" => A, "bra" => ϕ_dag, "ket" => ψ) - return logscalar(alg, tn; partitioned_vertices, kwargs...) +function loginner( + alg::Algorithm"bp", + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + partitioned_verts=default_inner_partitioned_vertices, + link_index_map=sim, + kwargs..., +) + tn = disjoint_union("bra" => link_index_map(dag(ϕ); sites=[]), "ket" => ψ) + return logscalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end -inner(alg::Algorithm"bp", ϕ::AbstractITensorNetwork, - ψ::AbstractITensorNetwork; kwargs...) = exp(loginner(alg, ϕ, ψ; kwargs...)) +function loginner( + alg::Algorithm"bp", + A::AbstractITensorNetwork, + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + partitioned_verts=default_inner_partitioned_vertices, + site_index_map=prime, + link_index_map=sim, + kwargs..., +) + ϕ_dag = link_index_map(site_index_map(dag(ϕ); links=[]); sites=[]) + tn = disjoint_union("operator" => A, "bra" => ϕ_dag, "ket" => ψ) + return logscalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) +end +function inner( + alg::Algorithm"bp", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; kwargs... +) + return exp(loginner(alg, ϕ, ψ; kwargs...)) +end +function inner( + alg::Algorithm"bp", + A::AbstractITensorNetwork, + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + kwargs..., +) + return exp(loginner(alg, A, ϕ, ψ; kwargs...)) +end +# TODO: rename `sqnorm` to match https://github.com/JuliaStats/Distances.jl, +# or `norm_sqr` to match `LinearAlgebra.norm_sqr` +function norm_sqr(ψ::AbstractITensorNetwork; kwargs...) + ψ_mapped = sim(ψ; sites=[]) + return inner(ψ, ψ_mapped; kwargs...) +end diff --git a/test/test_apply.jl b/test/test_apply.jl index 126aebf0..e4dc22f1 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -2,7 +2,7 @@ using ITensorNetworks using ITensorNetworks: incoming_messages, update, - contract_inner, + inner, norm_network, BeliefPropagationCache, VidalITensorNetwork @@ -64,15 +64,11 @@ using SplitApplyCombine print_fidelity_loss=true, envisposdef=true, ) - fSBP = - inner(ψOSBP, ψOexact) / - sqrt(inner(ψOexact, ψOexact) * inner(ψOSBP, ψOSBP)) + fSBP = inner(ψOSBP, ψOexact) / sqrt(inner(ψOexact, ψOexact) * inner(ψOSBP, ψOSBP)) fVidal = inner(ψOVidal_symm, ψOexact) / sqrt(inner(ψOexact, ψOexact) * inner(ψOVidal_symm, ψOVidal_symm)) - fGBP = - inner(ψOGBP, ψOexact) / - sqrt(inner(ψOexact, ψOexact) * inner(ψOGBP, ψOGBP)) + fGBP = inner(ψOGBP, ψOexact) / sqrt(inner(ψOexact, ψOexact) * inner(ψOGBP, ψOGBP)) @test real(fGBP * conj(fGBP)) >= real(fSBP * conj(fSBP)) diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index 2085589e..ea03a38e 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -2,7 +2,7 @@ using ITensorNetworks using ITensorNetworks: ising_network, split_index, - contract_inner, + inner, contract_boundary_mps, BeliefPropagationCache, tensornetwork, @@ -86,7 +86,8 @@ ITensors.disable_warn_order() ψOψ = ising_network(s, beta; szverts=vs) contract_seq = contraction_sequence(ψψ) - actual_szsz = contract(ψOψ; sequence=contract_seq)[] / contract(ψψ; sequence=contract_seq)[] + actual_szsz = + contract(ψOψ; sequence=contract_seq)[] / contract(ψψ; sequence=contract_seq)[] bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) bpc = update(bpc; maxiter=20) @@ -111,9 +112,7 @@ ITensors.disable_warn_order() ψψsplit = split_index(ψψ, NamedEdge.([(v, 1) => (v, 2) for v in vs])) env_tensors = incoming_messages(bpc, [(v, 2) for v in vs]) - rdm = contract( - vcat(env_tensors, ITensor[ψψsplit[vp] for vp in [(v, 2) for v in vs]]) - ) + rdm = contract(vcat(env_tensors, ITensor[ψψsplit[vp] for vp in [(v, 2) for v in vs]])) rdm = array((rdm * combiner(inds(rdm; plev=0)...)) * combiner(inds(rdm; plev=1)...)) rdm /= tr(rdm) diff --git a/test/test_gauging.jl b/test/test_gauging.jl index f0a7d10b..db6083d7 100644 --- a/test/test_gauging.jl +++ b/test/test_gauging.jl @@ -1,7 +1,7 @@ using ITensors using ITensorNetworks using ITensorNetworks: - contract_inner, gauge_error, update, messages, BeliefPropagationCache, VidalITensorNetwork + inner, gauge_error, update, messages, BeliefPropagationCache, VidalITensorNetwork using NamedGraphs using Test using Compat @@ -28,8 +28,7 @@ using SplitApplyCombine bp_cache = cache_ref[] # Test we just did a gauge transform and didn't change the overall network - @test contract_inner(ψ_symm, ψ) / - sqrt(contract_inner(ψ_symm, ψ_symm) * contract_inner(ψ, ψ)) ≈ 1.0 + @test inner(ψ_symm, ψ) / sqrt(inner(ψ_symm, ψ_symm) * inner(ψ, ψ)) ≈ 1.0 #Test all message tensors are approximately diagonal even when we keep running BP bp_cache = update(bp_cache; maxiter=20) diff --git a/test/test_tno.jl b/test/test_tno.jl index 115bf377..d67a0208 100644 --- a/test/test_tno.jl +++ b/test/test_tno.jl @@ -3,7 +3,7 @@ using ITensorNetworks using ITensors using Random -using ITensorNetworks: gate_group_to_tno, get_tnos, group_commuting_itensors, contract_inner +using ITensorNetworks: gate_group_to_tno, get_tnos, group_commuting_itensors, inner @testset "TN operator Basics" begin L = 3 @@ -42,12 +42,12 @@ using ITensorNetworks: gate_group_to_tno, get_tnos, group_commuting_itensors, co noprime!(ψ_tno[v]) end - z1 = contract_inner(ψ_gated, ψ_gated) - z2 = contract_inner(ψ_tnod, ψ_tnod) - z3 = contract_inner(ψ_tno, ψ_tno) - f12 = contract_inner(ψ_tnod, ψ_gated) / sqrt(z1 * z2) - f13 = contract_inner(ψ_tno, ψ_gated) / sqrt(z1 * z3) - f23 = contract_inner(ψ_tno, ψ_tnod) / sqrt(z2 * z3) + z1 = inner(ψ_gated, ψ_gated) + z2 = inner(ψ_tnod, ψ_tnod) + z3 = inner(ψ_tno, ψ_tno) + f12 = inner(ψ_tnod, ψ_gated) / sqrt(z1 * z2) + f13 = inner(ψ_tno, ψ_gated) / sqrt(z1 * z3) + f23 = inner(ψ_tno, ψ_tnod) / sqrt(z2 * z3) @test f12 * conj(f12) ≈ 1.0 @test f13 * conj(f13) ≈ 1.0 @test f23 * conj(f23) ≈ 1.0 From 1af45d85fbab33a052a7305cfc9dc7a84582d845 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 25 Mar 2024 16:44:31 -0400 Subject: [PATCH 03/35] Tree tests now pass to new inner function --- src/abstractitensornetwork.jl | 2 + src/caches/beliefpropagationcache.jl | 34 +++++++++++-- src/contract.jl | 38 +++++++++----- src/imports.jl | 1 + src/inner.jl | 50 +++++++++++++----- .../abstracttreetensornetwork.jl | 30 +++++------ test/test_apply.jl | 7 +-- test/test_belief_propagation.jl | 2 +- test/test_inner.jl | 51 +++++++++++++++++++ test/test_itensornetwork.jl | 14 ++--- 10 files changed, 172 insertions(+), 57 deletions(-) create mode 100644 test/test_inner.jl diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index 79768243..d0dcd6ed 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -903,6 +903,8 @@ end +(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork) = add(tn1, tn2) +ITensors.hasqns(tn::AbstractITensorNetwork) = all([hasqns(tn[v]) for v in vertices(tn)]) + ## # TODO: should this make sure that internal indices ## # don't clash? ## function hvncat( diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index 8b282d6a..a72edf20 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -13,10 +13,12 @@ end default_partitioned_vertices(ψ::AbstractITensorNetwork) = group(v -> v, vertices(ψ)) default_cache_update_kwargs(cache) = (; maxiter=20, tol=1e-5) +#TODO: Define a version of this that works for QN supporting tensors function message_diff(message_a::Vector{ITensor}, message_b::Vector{ITensor}) lhs, rhs = contract(message_a), contract(message_b) - return 0.5 * - norm((denseblocks(lhs) / sum(diag(lhs))) - (denseblocks(rhs) / sum(diag(rhs)))) + tr_lhs = length(inds(lhs)) == 1 ? sum(lhs) : sum(diag(lhs)) + tr_rhs = length(inds(rhs)) == 1 ? sum(rhs) : sum(diag(rhs)) + return 0.5 * norm((denseblocks(lhs) / tr_lhs) - (denseblocks(rhs) / tr_rhs)) end struct BeliefPropagationCache{PTN,MTS,DM} @@ -205,7 +207,7 @@ function update( verbose=false, kwargs..., ) - compute_error = !isnothing(tol) + compute_error = !isnothing(tol) && !hasqns(tensornetwork(bp_cache)) diff = compute_error ? Ref(0.0) : nothing if isnothing(maxiter) error("You need to specify a number of iterations for BP!") @@ -241,3 +243,29 @@ end function update_factor(bp_cache, vertex, factor) return update_factors(bp_cache, [vertex], ITensor[factor]) end + +function vertex_norm(bp_cache::BeliefPropagationCache, pv::PartitionVertex) + incoming_mts = environment(bp_cache, [pv]) + local_state = factor(bp_cache, pv) + return scalar(vcat(incoming_mts, local_state)) +end + +function edge_norm(bp_cache::BeliefPropagationCache, pe::PartitionEdge) + return scalar(vcat(message(bp_cache, pe), message(bp_cache, reverse(pe)))) +end + +function vertex_norms(bp_cache::BeliefPropagationCache, pvs::Vector) + return [vertex_norm(bp_cache, pv) for pv in pvs] +end + +function vertex_norms(bp_cache::BeliefPropagationCache) + return vertex_norms(bp_cache, partitionvertices(partitioned_itensornetwork(bp_cache))) +end + +function edge_norms(bp_cache::BeliefPropagationCache, pes::Vector) + return [edge_norm(bp_cache, pe) for pe in pes] +end + +function edge_norms(bp_cache::BeliefPropagationCache) + return edge_norms(bp_cache, partitionedges(partitioned_itensornetwork(bp_cache))) +end diff --git a/src/contract.jl b/src/contract.jl index 4d2d9b58..db190ab8 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -33,15 +33,16 @@ function contract_density_matrix( return out end -function scalar(tn::Union{AbstractITensorNetwork,Vector{ITensor}}; alg="exact", kwargs...) - return scalar(Algorithm(alg), tn; kwargs...) -end function scalar( alg::Algorithm, tn::Union{AbstractITensorNetwork,Vector{ITensor}}; kwargs... ) return contract(alg, tn; kwargs...)[] end +function scalar(tn::Union{AbstractITensorNetwork,Vector{ITensor}}; alg="exact", kwargs...) + return scalar(Algorithm(alg), tn; kwargs...) +end + function logscalar( tn::Union{AbstractITensorNetwork,Vector{ITensor}}; alg::String="exact", kwargs... ) @@ -70,19 +71,28 @@ function logscalar( cache![] = update(cache![]; cache_update_kwargs...) end - pg = partitioned_itensornetwork(cache![]) + numerator_terms, denominator_terms = vertex_norms(cache![]), edge_norms(cache![]) + + return sum(log.(complex.(numerator_terms))) - sum(log.(complex.((denominator_terms)))) +end - log_numerator, log_denominator = 0, 0 - for pv in partitionvertices(pg) - incoming_mts = incoming_messages(cache![], [pv]) - local_state = factor(cache![], pv) - log_numerator += logscalar(vcat(incoming_mts, local_state); alg="exact") +function scalar( + alg::Algorithm"bp", + tn::AbstractITensorNetwork; + (cache!)=nothing, + partitioned_vertices=default_partitioned_vertices(tn), + update_cache=isnothing(cache!), + cache_update_kwargs=default_cache_update_kwargs(cache!), +) + if isnothing(cache!) + cache! = Ref(BeliefPropagationCache(tn, partitioned_vertices)) end - for pe in partitionedges(pg) - log_denominator += logscalar( - vcat(message(cache![], pe), message(cache![], reverse(pe))); alg="exact" - ) + + if update_cache + cache![] = update(cache![]; cache_update_kwargs...) end - return log_numerator - log_denominator + numerator_terms, denominator_terms = vertex_norms(cache![]), edge_norms(cache![]) + + return prod(numerator_terms) / prod(denominator_terms) end diff --git a/src/imports.jl b/src/imports.jl index 72a07910..42f08594 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -94,6 +94,7 @@ import ITensors: nsite, # promotion and conversion promote_itensor_eltype, + scalar, scalartype, #adding add diff --git a/src/inner.jl b/src/inner.jl index b1d18d3d..28f8d0c3 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -11,8 +11,8 @@ function inner( end function inner( - A::AbstractITensorNetwork, ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; alg=default_algorithm([A, ϕ, ψ]), kwargs..., @@ -37,15 +37,15 @@ end function inner( alg::Algorithm"exact", - A::AbstractITensorNetwork, ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; sequence=nothing, contraction_sequence_kwargs=(;), site_index_map=prime, kwargs..., ) - tn = flatten_networks(site_index_map(dag(ϕ); links=[]), A, ψ; kwargs...) + tn = flatten_networks(dag(ϕ), A, ψ; kwargs...) if isnothing(sequence) sequence = contraction_sequence(tn; contraction_sequence_kwargs...) end @@ -60,6 +60,17 @@ function loginner( ) return loginner(Algorithm(alg), ϕ, ψ; kwargs...) end + +function loginner( + ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + alg=default_algorithm(ϕ, ψ), + kwargs..., +) + return loginner(Algorithm(alg), A, ϕ, ψ; kwargs...) +end + function loginner( alg::Algorithm"exact", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; kwargs... ) @@ -74,29 +85,36 @@ function loginner( link_index_map=sim, kwargs..., ) - tn = disjoint_union("bra" => link_index_map(dag(ϕ); sites=[]), "ket" => ψ) + ϕ_dag = map_inds(link_index_map, dag(ϕ); sites=[]) + tn = disjoint_union("bra" => ϕ_dag, "ket" => ψ) return logscalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end function loginner( alg::Algorithm"bp", - A::AbstractITensorNetwork, ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; partitioned_verts=default_inner_partitioned_vertices, - site_index_map=prime, link_index_map=sim, kwargs..., ) - ϕ_dag = link_index_map(site_index_map(dag(ϕ); links=[]); sites=[]) + ϕ_dag = map_inds(link_index_map, dag(ϕ); sites=[]) tn = disjoint_union("operator" => A, "bra" => ϕ_dag, "ket" => ψ) return logscalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end function inner( - alg::Algorithm"bp", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; kwargs... + alg::Algorithm"bp", + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + partitioned_verts=default_inner_partitioned_vertices, + link_index_map=sim, + kwargs..., ) - return exp(loginner(alg, ϕ, ψ; kwargs...)) + ϕ_dag = map_inds(link_index_map, dag(ϕ); sites=[]) + tn = disjoint_union("bra" => ϕ_dag, "ket" => ψ) + return scalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end function inner( @@ -104,14 +122,22 @@ function inner( A::AbstractITensorNetwork, ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; + partitioned_verts=default_inner_partitioned_vertices, + link_index_map=sim, kwargs..., ) - return exp(loginner(alg, A, ϕ, ψ; kwargs...)) + ϕ_dag = map_inds(link_index_map, dag(ϕ); sites=[]) + tn = disjoint_union("operator" => A, "bra" => ϕ_dag, "ket" => ψ) + return scalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end # TODO: rename `sqnorm` to match https://github.com/JuliaStats/Distances.jl, # or `norm_sqr` to match `LinearAlgebra.norm_sqr` -function norm_sqr(ψ::AbstractITensorNetwork; kwargs...) - ψ_mapped = sim(ψ; sites=[]) +function norm_sqr(ψ::AbstractITensorNetwork; link_index_map=sim, kwargs...) + ψ_mapped = link_index_map(ψ; sites=[]) return inner(ψ, ψ_mapped; kwargs...) end + +function norm(ψ::AbstractITensorNetwork; kwargs...) + return sqrt(abs(real(norm_sqr(ψ; kwargs...)))) +end diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl index 289449ab..e2cd5754 100644 --- a/src/treetensornetworks/abstracttreetensornetwork.jl +++ b/src/treetensornetworks/abstracttreetensornetwork.jl @@ -232,7 +232,7 @@ function logdot(ψ1::TTNT, ψ2::TTNT; kwargs...) where {TTNT<:AbstractTTN} return loginner(ψ1, ψ2; kwargs...) end -# TODO: stick with this traversal or find optimal contraction sequence? +#TODO: stick with this traversal or find optimal contraction sequence? function loginner( ψ1::TTNT, ψ2::TTNT; root_vertex=default_root_vertex(ψ1, ψ2) )::Number where {TTNT<:AbstractTTN} @@ -366,20 +366,20 @@ end # # TODO: implement using multi-graph disjoint union -function inner( - y::AbstractTTN, A::AbstractTTN, x::AbstractTTN; root_vertex=default_root_vertex(x, A, y) -) - traversal_order = reverse(post_order_dfs_vertices(x, root_vertex)) - check_hascommoninds(siteinds, A, x) - check_hascommoninds(siteinds, A, y) - ydag = sim(dag(y); sites=[]) - x = sim(x; sites=[]) - O = ydag[root_vertex] * A[root_vertex] * x[root_vertex] - for v in traversal_order[2:end] - O = O * ydag[v] * A[v] * x[v] - end - return O[] -end +# function inner( +# y::AbstractTTN, A::AbstractTTN, x::AbstractTTN; root_vertex=default_root_vertex(x, A, y) +# ) +# traversal_order = reverse(post_order_dfs_vertices(x, root_vertex)) +# check_hascommoninds(siteinds, A, x) +# check_hascommoninds(siteinds, A, y) +# ydag = sim(dag(y); sites=[]) +# x = sim(x; sites=[]) +# O = ydag[root_vertex] * A[root_vertex] * x[root_vertex] +# for v in traversal_order[2:end] +# O = O * ydag[v] * A[v] * x[v] +# end +# return O[] +# end # TODO: implement using multi-graph disjoint function inner( diff --git a/test/test_apply.jl b/test/test_apply.jl index dbbcab62..c501faa1 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,11 +1,6 @@ using ITensorNetworks using ITensorNetworks: - environment, - update, - inner, - norm_network, - BeliefPropagationCache, - VidalITensorNetwork + environment, update, inner, norm_network, BeliefPropagationCache, VidalITensorNetwork using Test using Compat using ITensors diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index f8b42882..213257e2 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -112,7 +112,7 @@ ITensors.disable_warn_order() bpc = update(bpc; maxiter=20) ψψsplit = split_index(ψψ, NamedEdge.([(v, 1) => (v, 2) for v in vs])) - env_tensors = incoming_messages(bpc, [(v, 2) for v in vs]) + env_tensors = environment(bpc, [(v, 2) for v in vs]) rdm = contract(vcat(env_tensors, ITensor[ψψsplit[vp] for vp in [(v, 2) for v in vs]])) rdm = array((rdm * combiner(inds(rdm; plev=0)...)) * combiner(inds(rdm; plev=1)...)) diff --git a/test/test_inner.jl b/test/test_inner.jl new file mode 100644 index 00000000..8f919a8e --- /dev/null +++ b/test/test_inner.jl @@ -0,0 +1,51 @@ +using Test +using Random +using Graphs +using NamedGraphs +using ITensorNetworks +using SplitApplyCombine + +using ITensorNetworks: logscalar, scalar, inner +using ITensors: siteinds + +@testset "Inner products, BP vs exact comparison" begin + Random.seed!(1234) + L = 12 + χ = 2 + g = NamedGraph(Graphs.SimpleGraph(uniform_tree(L))) + s = siteinds("S=1/2", g) + y = randomITensorNetwork(s; link_space=χ) + x = randomITensorNetwork(s; link_space=χ) + + #First lets do it with the flattened version of the network + xy = inner_network(x, y; flatten=true, combine_linkinds=true) + xy_scalar = scalar(xy) + xy_scalar_bp = scalar(xy; alg="bp", partitioned_vertices=group(v -> v, vertices(xy))) + xy_scalar_logbp = exp( + logscalar(xy; alg="bp", partitioned_vertices=group(v -> v, vertices(xy))) + ) + + @test xy_scalar ≈ xy_scalar_bp + @test xy_scalar_bp ≈ xy_scalar_logbp + @test xy_scalar ≈ xy_scalar_logbp + + #Now lets keep it unflattened and do it via the inner function + xy_scalar = inner(x, y; alg="exact") + xy_scalar_bp = inner(x, y; alg="bp") + xy_scalar_logbp = exp(loginner(x, y; alg="bp")) + + @test xy_scalar ≈ xy_scalar_bp + @test xy_scalar_bp ≈ xy_scalar_logbp + @test xy_scalar ≈ xy_scalar_logbp + + #test contraction of three layers for expectation values + A = TTN(ITensorNetworks.heisenberg(s), s) + xAy_scalar = inner(x', A, y; alg="exact") + xAy_scalar_bp = inner(x', A, y; alg="bp") + xAy_scalar_logbp = exp(loginner(x', A, y; alg="bp")) + + @test xAy_scalar ≈ xAy_scalar_bp + @test xAy_scalar_bp ≈ xAy_scalar_logbp + @test xAy_scalar ≈ xAy_scalar_logbp +end +nothing diff --git a/test/test_itensornetwork.jl b/test/test_itensornetwork.jl index 5f09ecd8..b1e6c097 100644 --- a/test/test_itensornetwork.jl +++ b/test/test_itensornetwork.jl @@ -7,6 +7,8 @@ using NamedGraphs using Random using Test +using ITensorNetworks: norm_sqr, neighbor_itensors + @testset "ITensorNetwork tests" begin @testset "ITensorNetwork Basics" begin Random.seed!(1234) @@ -146,26 +148,26 @@ using Test @testset "orthogonalize" begin tn = randomITensorNetwork(named_grid(4); link_space=2) - Z = contract(inner_network(tn, tn))[] + Z = norm_sqr(tn) tn_ortho = factorize(tn, 4 => 3) # TODO: Error here in arranging the edges. Arrange by hash? - Z̃ = contract(inner_network(tn_ortho, tn_ortho))[] + Z̃ = norm_sqr(tn_ortho) @test nv(tn_ortho) == 5 @test nv(tn) == 4 @test Z ≈ Z̃ tn_ortho = orthogonalize(tn, 4 => 3) - Z̃ = contract(inner_network(tn_ortho, tn_ortho))[] + Z̃ = norm_sqr(tn_ortho) @test nv(tn_ortho) == 4 @test nv(tn) == 4 @test Z ≈ Z̃ tn_ortho = orthogonalize(tn, 1) - Z̃ = contract(inner_network(tn_ortho, tn_ortho))[] + Z̃ = norm_sqr(tn_ortho) @test Z ≈ Z̃ - Z̃ = contract(inner_network(tn_ortho, tn))[] + Z̃ = inner(tn_ortho, tn) @test Z ≈ Z̃ end @@ -202,7 +204,7 @@ using Test s = siteinds("S=1/2", g) ψ = ITensorNetwork(s; link_space=2) - nt = ITensorNetworks.neighbor_itensors(ψ, (1, 1)) + nt = neighbor_itensors(ψ, (1, 1)) @test length(nt) == 2 @test all(map(hascommoninds(ψ[1, 1]), nt)) From eb3fa773fa0a9439c01e1b46930ea0f8af2b895d Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 26 Mar 2024 14:04:47 -0400 Subject: [PATCH 04/35] Bug Fix. --- src/abstractitensornetwork.jl | 59 +++++++++++++------ src/inner.jl | 21 ++++--- test/test_inner.jl | 12 ++-- .../test_solvers/test_dmrg_x.jl | 6 +- 4 files changed, 63 insertions(+), 35 deletions(-) diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index d0dcd6ed..9dbf1075 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -643,10 +643,10 @@ function flatten_networks( combine_linkinds=true, kwargs..., ) - @assert issetequal(vertices(tn1), vertices(tn2)) tn1 = map_bra_linkinds(tn1; sites=[]) flattened_net = ⊗(tn1, tn2; kwargs...) if flatten + @assert issetequal(vertices(tn1), vertices(tn2)) for v in vertices(tn1) flattened_net = contract(flattened_net, (v, 2) => (v, 1); merged_vertex=v) end @@ -667,29 +667,50 @@ function flatten_networks( return flatten_networks(flatten_networks(tn1, tn2; kwargs...), tn3, tn_tail...; kwargs...) end -#Ideally this will dispatch to inner_network but this is a temporary fast version for now -function norm_network(tn::AbstractITensorNetwork) - tnbra = rename_vertices(v -> (v, 1), data_graph(tn)) - tndag = copy(tn) - for v in vertices(tndag) - setindex_preserve_graph!(tndag, dag(tndag[v]), v) - end - tnket = rename_vertices(v -> (v, 2), data_graph(prime(tndag; sites=[]))) - tntn = ITensorNetwork(union(tnbra, tnket)) - for v in vertices(tn) - if !isempty(commoninds(tntn[(v, 1)], tntn[(v, 2)])) - add_edge!(tntn, (v, 1) => (v, 2)) - end - end - return tntn +function inner_network( + x::AbstractITensorNetwork, + y::AbstractITensorNetwork; + flatten=false, + combine_linkinds=false, + map_bra_linkinds=prime, + kwargs..., +) + return flatten_networks(dag(x), y; flatten, combine_linkinds, map_bra_linkinds, kwargs...) end function inner_network( - ψ::AbstractITensorNetwork, ϕ::AbstractITensorNetwork; flatten=false, kwargs... + x::AbstractITensorNetwork, + A::AbstractITensorNetwork, + y::AbstractITensorNetwork; + flatten=false, + combine_linkinds=false, + map_bra_linkinds=prime, + kwargs..., ) - return flatten_networks(ψ, dag(ϕ); flatten, kwargs...) + return flatten_networks( + dag(x), A, y; flatten, combine_linkinds, map_bra_linkinds, kwargs... + ) +end + +inner_network(x::AbstractITensorNetwork; kwargs...) = inner_network(x, x; kwargs...) +const norm_sqr_network = inner_network + +#Ideally this will not be necessary but this is a temporary fast version to avoid the overhead of `disjoint_union` +function norm_sqr_network_fast(ψ::AbstractITensorNetwork) + ψbra = rename_vertices(v -> (v, 1), data_graph(ψ)) + ψdag = copy(ψ) + for v in vertices(ψdag) + setindex_preserve_graph!(ψdag, dag(ψdag[v]), v) + end + ψket = rename_vertices(v -> (v, 2), data_graph(prime(ψdag; sites=[]))) + ψψ = ITensorNetwork(union(ψbra, ψket)) + for v in vertices(ψ) + if !isempty(commoninds(ψψ[(v, 1)], ψψ[(v, 2)])) + add_edge!(ψψ, (v, 1) => (v, 2)) + end + end + return ψψ end -norm_sqr_network(ψ::AbstractITensorNetwork; kwargs...) = inner_network(ψ, ψ; kwargs...) # # Printing diff --git a/src/inner.jl b/src/inner.jl index 28f8d0c3..2a8279fb 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -45,7 +45,7 @@ function inner( site_index_map=prime, kwargs..., ) - tn = flatten_networks(dag(ϕ), A, ψ; kwargs...) + tn = inner_network(ϕ, A, ψ; kwargs...) if isnothing(sequence) sequence = contraction_sequence(tn; contraction_sequence_kwargs...) end @@ -77,6 +77,16 @@ function loginner( return log(inner(alg, ϕ, ψ); kwargs...) end +function loginner( + alg::Algorithm"exact", + ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + kwargs..., +) + return log(inner(alg, ϕ, A, ψ); kwargs...) +end + function loginner( alg::Algorithm"bp", ϕ::AbstractITensorNetwork, @@ -109,7 +119,7 @@ function inner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; partitioned_verts=default_inner_partitioned_vertices, - link_index_map=sim, + link_index_map=prime, kwargs..., ) ϕ_dag = map_inds(link_index_map, dag(ϕ); sites=[]) @@ -123,7 +133,7 @@ function inner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; partitioned_verts=default_inner_partitioned_vertices, - link_index_map=sim, + link_index_map=prime, kwargs..., ) ϕ_dag = map_inds(link_index_map, dag(ϕ); sites=[]) @@ -133,10 +143,7 @@ end # TODO: rename `sqnorm` to match https://github.com/JuliaStats/Distances.jl, # or `norm_sqr` to match `LinearAlgebra.norm_sqr` -function norm_sqr(ψ::AbstractITensorNetwork; link_index_map=sim, kwargs...) - ψ_mapped = link_index_map(ψ; sites=[]) - return inner(ψ, ψ_mapped; kwargs...) -end +norm_sqr(ψ::AbstractITensorNetwork; kwargs...) = inner(ψ, ψ; kwargs...) function norm(ψ::AbstractITensorNetwork; kwargs...) return sqrt(abs(real(norm_sqr(ψ; kwargs...)))) diff --git a/test/test_inner.jl b/test/test_inner.jl index 8f919a8e..bcb8e7e9 100644 --- a/test/test_inner.jl +++ b/test/test_inner.jl @@ -5,8 +5,8 @@ using NamedGraphs using ITensorNetworks using SplitApplyCombine -using ITensorNetworks: logscalar, scalar, inner -using ITensors: siteinds +using ITensorNetworks: logscalar, scalar, inner, loginner +using ITensors: siteinds, dag @testset "Inner products, BP vs exact comparison" begin Random.seed!(1234) @@ -18,7 +18,7 @@ using ITensors: siteinds x = randomITensorNetwork(s; link_space=χ) #First lets do it with the flattened version of the network - xy = inner_network(x, y; flatten=true, combine_linkinds=true) + xy = flatten_networks(dag(x), y; flatten=true, combine_linkinds=true) xy_scalar = scalar(xy) xy_scalar_bp = scalar(xy; alg="bp", partitioned_vertices=group(v -> v, vertices(xy))) xy_scalar_logbp = exp( @@ -29,8 +29,8 @@ using ITensors: siteinds @test xy_scalar_bp ≈ xy_scalar_logbp @test xy_scalar ≈ xy_scalar_logbp - #Now lets keep it unflattened and do it via the inner function - xy_scalar = inner(x, y; alg="exact") + #Now lets do it via the inner function + xy_scalar = inner(x, y; alg="exact", flatten=true, combine_linkinds=true) xy_scalar_bp = inner(x, y; alg="bp") xy_scalar_logbp = exp(loginner(x, y; alg="bp")) @@ -40,7 +40,7 @@ using ITensors: siteinds #test contraction of three layers for expectation values A = TTN(ITensorNetworks.heisenberg(s), s) - xAy_scalar = inner(x', A, y; alg="exact") + xAy_scalar = inner(x', A, y; alg="exact", flatten=true, combine_linkinds=true) xAy_scalar_bp = inner(x', A, y; alg="bp") xAy_scalar_logbp = exp(loginner(x', A, y; alg="bp")) diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl index 18c90539..8b2b1a45 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl @@ -20,9 +20,9 @@ using Test ϕ = dmrg_x(H, ψ; nsites=2, dmrg_x_kwargs...) - @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ', H, ϕ) / inner(ϕ, ϕ) rtol = 1e-1 - @test inner(H, ψ, H, ψ) ≉ inner(ψ', H, ψ)^2 rtol = 1e-7 - @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ', H, ϕ)^2 rtol = 1e-7 + @test inner(H, ψ', ψ) / inner(ψ, ψ) ≈ inner(H, ϕ', ϕ) / inner(ϕ, ϕ) rtol = 1e-1 + @test inner(H, ψ, H, ψ) ≉ inner(H, ψ', ψ)^2 rtol = 1e-7 + @test inner(H, ϕ, H, ϕ) ≈ inner(H, ϕ', ϕ)^2 rtol = 1e-7 ϕ̃ = dmrg_x(H, ϕ; nsites=1, dmrg_x_kwargs...) From e6e8c6f6fce40eefbb9e935f72683cdc5a0bf6ef Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 26 Mar 2024 14:07:49 -0400 Subject: [PATCH 05/35] Argument rearrange DMRGx --- test/test_treetensornetworks/test_solvers/test_dmrg_x.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl index 8b2b1a45..2eba03bb 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl @@ -20,9 +20,9 @@ using Test ϕ = dmrg_x(H, ψ; nsites=2, dmrg_x_kwargs...) - @test inner(H, ψ', ψ) / inner(ψ, ψ) ≈ inner(H, ϕ', ϕ) / inner(ϕ, ϕ) rtol = 1e-1 - @test inner(H, ψ, H, ψ) ≉ inner(H, ψ', ψ)^2 rtol = 1e-7 - @test inner(H, ϕ, H, ϕ) ≈ inner(H, ϕ', ϕ)^2 rtol = 1e-7 + @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ', H, ϕ) / inner(ϕ, ϕ) rtol = 1e-1 + @test inner(H, ψ, H, ψ) ≉ inner(ψ', H, ψ)^2 rtol = 1e-7 + @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ',H, ϕ)^2 rtol = 1e-7 ϕ̃ = dmrg_x(H, ϕ; nsites=1, dmrg_x_kwargs...) From c6f4dac653a6e749745b71fa3df6a654a20fe74e Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 26 Mar 2024 14:10:08 -0400 Subject: [PATCH 06/35] Reinstate Tree Inner for TTN Type --- .../abstracttreetensornetwork.jl | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl index e2cd5754..3bbcf32c 100644 --- a/src/treetensornetworks/abstracttreetensornetwork.jl +++ b/src/treetensornetworks/abstracttreetensornetwork.jl @@ -366,20 +366,20 @@ end # # TODO: implement using multi-graph disjoint union -# function inner( -# y::AbstractTTN, A::AbstractTTN, x::AbstractTTN; root_vertex=default_root_vertex(x, A, y) -# ) -# traversal_order = reverse(post_order_dfs_vertices(x, root_vertex)) -# check_hascommoninds(siteinds, A, x) -# check_hascommoninds(siteinds, A, y) -# ydag = sim(dag(y); sites=[]) -# x = sim(x; sites=[]) -# O = ydag[root_vertex] * A[root_vertex] * x[root_vertex] -# for v in traversal_order[2:end] -# O = O * ydag[v] * A[v] * x[v] -# end -# return O[] -# end +function inner( + y::AbstractTTN, A::AbstractTTN, x::AbstractTTN; root_vertex=default_root_vertex(x, A, y) +) + traversal_order = reverse(post_order_dfs_vertices(x, root_vertex)) + check_hascommoninds(siteinds, A, x) + check_hascommoninds(siteinds, A, y) + ydag = sim(dag(y); sites=[]) + x = sim(x; sites=[]) + O = ydag[root_vertex] * A[root_vertex] * x[root_vertex] + for v in traversal_order[2:end] + O = O * ydag[v] * A[v] * x[v] + end + return O[] +end # TODO: implement using multi-graph disjoint function inner( From 372f9b4e905880bc28c2695ceb8f829fbe495f14 Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 26 Mar 2024 14:11:17 -0400 Subject: [PATCH 07/35] Format --- test/test_treetensornetworks/test_solvers/test_dmrg_x.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl index 2eba03bb..18c90539 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl @@ -22,7 +22,7 @@ using Test @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ', H, ϕ) / inner(ϕ, ϕ) rtol = 1e-1 @test inner(H, ψ, H, ψ) ≉ inner(ψ', H, ψ)^2 rtol = 1e-7 - @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ',H, ϕ)^2 rtol = 1e-7 + @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ', H, ϕ)^2 rtol = 1e-7 ϕ̃ = dmrg_x(H, ϕ; nsites=1, dmrg_x_kwargs...) From fbe906395d76f44950dfa79e7063b14acc229721 Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 26 Mar 2024 14:49:49 -0400 Subject: [PATCH 08/35] Remove import of check_hascommoninds. Fix norm naming --- src/exports.jl | 2 -- src/gauging.jl | 2 +- src/imports.jl | 1 - test/test_apply.jl | 4 ++-- 4 files changed, 3 insertions(+), 6 deletions(-) diff --git a/src/exports.jl b/src/exports.jl index 09dbc3bd..0b363757 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -54,9 +54,7 @@ export AbstractITensorNetwork, default_bp_cache, flatten_networks, inner_network, - norm_network, factorize!, - norm_sqr_network, linkinds_combiners, combine_linkinds, externalinds, diff --git a/src/gauging.jl b/src/gauging.jl index 89a30555..625b1543 100644 --- a/src/gauging.jl +++ b/src/gauging.jl @@ -20,7 +20,7 @@ function copy(ψ::VidalITensorNetwork) end function default_norm_cache(ψ::ITensorNetwork) - ψψ = norm_network(ψ) + ψψ = norm_sqr_network(ψ) return BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) end diff --git a/src/imports.jl b/src/imports.jl index 42f08594..999201d8 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -83,7 +83,6 @@ import ITensors: # permute permute, #commoninds - check_hascommoninds, hascommoninds, # linkdims linkdim, diff --git a/test/test_apply.jl b/test/test_apply.jl index c501faa1..3c0644a0 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,6 +1,6 @@ using ITensorNetworks using ITensorNetworks: - environment, update, inner, norm_network, BeliefPropagationCache, VidalITensorNetwork + environment, update, inner, norm_sqr_network, BeliefPropagationCache, VidalITensorNetwork using Test using Compat using ITensors @@ -19,7 +19,7 @@ using SplitApplyCombine χ = 2 ψ = randomITensorNetwork(s; link_space=χ) v1, v2 = (2, 2), (1, 2) - ψψ = norm_network(ψ) + ψψ = norm_sqr_network(ψ) #Simple Belief Propagation Grouping bp_cache = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) From 19a158c1b5d20fb07f0fc34041473bdf4deafc7f Mon Sep 17 00:00:00 2001 From: Joey Date: Wed, 27 Mar 2024 12:12:57 -0400 Subject: [PATCH 09/35] Refactor flatten_networks --- src/ITensorNetworks.jl | 1 - src/abstractitensornetwork.jl | 85 ++++++++++++++++++----------------- src/inner.jl | 20 ++++----- test/test_inner.jl | 4 +- 4 files changed, 53 insertions(+), 57 deletions(-) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 300ab1d8..6a0b4409 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -43,7 +43,6 @@ using ITensors: AbstractMPS, Algorithm, OneITensor, - check_hascommoninds, commontags, dim, orthocenter, diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index 9dbf1075..cf86c020 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -635,61 +635,62 @@ function split_index( return tn end -function flatten_networks( - tn1::AbstractITensorNetwork, - tn2::AbstractITensorNetwork; - map_bra_linkinds=sim, - flatten=true, - combine_linkinds=true, - kwargs..., +function stack_flatten_combine( + tns::Vector{<:AbstractITensorNetwork}; + map_bra_linkinds=prime, + flatten=false, + combine_linkinds=false, ) - tn1 = map_bra_linkinds(tn1; sites=[]) - flattened_net = ⊗(tn1, tn2; kwargs...) - if flatten - @assert issetequal(vertices(tn1), vertices(tn2)) - for v in vertices(tn1) - flattened_net = contract(flattened_net, (v, 2) => (v, 1); merged_vertex=v) + tns = copy(tns) + stacked_tn = map_bra_linkinds(popfirst!(tns); sites=[]) + current_suffix = 1 + for tn in tns + stacked_tn_vertices = vertices(stacked_tn) + stacked_tn = disjoint_union(current_suffix => stacked_tn, current_suffix + 1 => tn) + + if flatten + @assert issetequal(vertices(tn), stacked_tn_vertices) + for v in vertices(tn) + stacked_tn = contract( + stacked_tn, (v, current_suffix + 1) => (v, current_suffix); merged_vertex=v + ) + end + else + if current_suffix != 1 + #Strip back + stacked_tn_vertices = [(v, current_suffix) for v in stacked_tn_vertices] + stacked_tn = rename_vertices( + v -> v ∈ stacked_tn_vertices ? first(v) : v, stacked_tn + ) + end + current_suffix += 1 + end + if combine_linkinds + stacked_tn = ITensorNetworks.combine_linkinds(stacked_tn) end end - if combine_linkinds - flattened_net = ITensorNetworks.combine_linkinds(flattened_net) - end - return flattened_net + + return stacked_tn +end + +function stack_flatten_combine(tns::AbstractITensorNetwork...; kwargs...) + return stack_flatten_combine([tns...]; kwargs...) end function flatten_networks( - tn1::AbstractITensorNetwork, - tn2::AbstractITensorNetwork, - tn3::AbstractITensorNetwork, - tn_tail::AbstractITensorNetwork...; - kwargs..., + tns::AbstractITensorNetwork...; combine_linkinds=true, map_bra_linkinds=prime ) - return flatten_networks(flatten_networks(tn1, tn2; kwargs...), tn3, tn_tail...; kwargs...) + return stack_flatten_combine(tns...; flatten=true, combine_linkinds, map_bra_linkinds) end -function inner_network( - x::AbstractITensorNetwork, - y::AbstractITensorNetwork; - flatten=false, - combine_linkinds=false, - map_bra_linkinds=prime, - kwargs..., -) - return flatten_networks(dag(x), y; flatten, combine_linkinds, map_bra_linkinds, kwargs...) +function inner_network(x::AbstractITensorNetwork, y::AbstractITensorNetwork; kwargs...) + return stack_flatten_combine(dag(x), y; kwargs...) end function inner_network( - x::AbstractITensorNetwork, - A::AbstractITensorNetwork, - y::AbstractITensorNetwork; - flatten=false, - combine_linkinds=false, - map_bra_linkinds=prime, - kwargs..., + x::AbstractITensorNetwork, A::AbstractITensorNetwork, y::AbstractITensorNetwork; kwargs... ) - return flatten_networks( - dag(x), A, y; flatten, combine_linkinds, map_bra_linkinds, kwargs... - ) + return stack_flatten_combine(dag(x), A, y; kwargs...) end inner_network(x::AbstractITensorNetwork; kwargs...) = inner_network(x, x; kwargs...) diff --git a/src/inner.jl b/src/inner.jl index 2a8279fb..7388299d 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -92,11 +92,10 @@ function loginner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; partitioned_verts=default_inner_partitioned_vertices, - link_index_map=sim, + map_bra_linkinds=sim, kwargs..., ) - ϕ_dag = map_inds(link_index_map, dag(ϕ); sites=[]) - tn = disjoint_union("bra" => ϕ_dag, "ket" => ψ) + tn = inner_network(ϕ, ψ; map_bra_linkinds) return logscalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end @@ -106,11 +105,10 @@ function loginner( A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; partitioned_verts=default_inner_partitioned_vertices, - link_index_map=sim, + map_bra_linkinds=sim, kwargs..., ) - ϕ_dag = map_inds(link_index_map, dag(ϕ); sites=[]) - tn = disjoint_union("operator" => A, "bra" => ϕ_dag, "ket" => ψ) + tn = inner_network(ϕ, A, ψ; map_bra_linkinds) return logscalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end @@ -119,11 +117,10 @@ function inner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; partitioned_verts=default_inner_partitioned_vertices, - link_index_map=prime, + map_bra_linkinds=prime, kwargs..., ) - ϕ_dag = map_inds(link_index_map, dag(ϕ); sites=[]) - tn = disjoint_union("bra" => ϕ_dag, "ket" => ψ) + tn = inner_network(ϕ, ψ; map_bra_linkinds) return scalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end @@ -133,11 +130,10 @@ function inner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; partitioned_verts=default_inner_partitioned_vertices, - link_index_map=prime, + map_bra_linkinds=prime, kwargs..., ) - ϕ_dag = map_inds(link_index_map, dag(ϕ); sites=[]) - tn = disjoint_union("operator" => A, "bra" => ϕ_dag, "ket" => ψ) + tn = inner_network(ϕ, A, ψ; map_bra_linkinds) return scalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end diff --git a/test/test_inner.jl b/test/test_inner.jl index bcb8e7e9..4da47060 100644 --- a/test/test_inner.jl +++ b/test/test_inner.jl @@ -18,7 +18,7 @@ using ITensors: siteinds, dag x = randomITensorNetwork(s; link_space=χ) #First lets do it with the flattened version of the network - xy = flatten_networks(dag(x), y; flatten=true, combine_linkinds=true) + xy = inner_network(x, y; combine_linkinds=true, flatten=true) xy_scalar = scalar(xy) xy_scalar_bp = scalar(xy; alg="bp", partitioned_vertices=group(v -> v, vertices(xy))) xy_scalar_logbp = exp( @@ -39,7 +39,7 @@ using ITensors: siteinds, dag @test xy_scalar ≈ xy_scalar_logbp #test contraction of three layers for expectation values - A = TTN(ITensorNetworks.heisenberg(s), s) + A = ITensorNetwork(TTN(ITensorNetworks.heisenberg(s), s)) xAy_scalar = inner(x', A, y; alg="exact", flatten=true, combine_linkinds=true) xAy_scalar_bp = inner(x', A, y; alg="bp") xAy_scalar_logbp = exp(loginner(x', A, y; alg="bp")) From 9971ff3dd85eeab4a0d3f72789c4217cb224c1de Mon Sep 17 00:00:00 2001 From: Joey Date: Wed, 27 Mar 2024 12:35:39 -0400 Subject: [PATCH 10/35] Remove set nsite and ProjMPS --- src/ITensorNetworks.jl | 4 +--- .../solvers/deprecated/projmpo_apply.jl | 7 +------ .../solvers/deprecated/projmpo_mps2.jl | 10 +--------- src/treetensornetworks/solvers/deprecated/projmps2.jl | 8 +------- 4 files changed, 4 insertions(+), 25 deletions(-) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 6a0b4409..e286ae79 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -45,9 +45,7 @@ using ITensors: OneITensor, commontags, dim, - orthocenter, - ProjMPS, - set_nsite! + orthocenter using KrylovKit: exponentiate, eigsolve, linsolve using NamedGraphs: AbstractNamedGraph, diff --git a/src/treetensornetworks/solvers/deprecated/projmpo_apply.jl b/src/treetensornetworks/solvers/deprecated/projmpo_apply.jl index d9733759..7b73ea13 100644 --- a/src/treetensornetworks/solvers/deprecated/projmpo_apply.jl +++ b/src/treetensornetworks/solvers/deprecated/projmpo_apply.jl @@ -1,4 +1,4 @@ -import ITensors: AbstractProjMPO, makeL!, makeR!, set_nsite! +import ITensors: AbstractProjMPO, makeL!, makeR! import Base: copy """ @@ -36,11 +36,6 @@ function copy(P::ProjMPOApply) return ProjMPOApply(P.lpos, P.rpos, P.nsite, copy(P.psi0), copy(P.H), copy(P.LR)) end -function set_nsite!(P::ProjMPOApply, nsite) - P.nsite = nsite - return P -end - function makeL!(P::ProjMPOApply, psi::MPS, k::Int) # Save the last `L` that is made to help with caching # for DiskProjMPO diff --git a/src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl b/src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl index a24255fe..7a619624 100644 --- a/src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl +++ b/src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl @@ -1,4 +1,4 @@ -import ITensors: AbstractProjMPO, makeL!, makeR!, set_nsite!, contract, nsite +import ITensors: AbstractProjMPO, makeL!, makeR!, contract, nsite import Base: copy mutable struct ProjMPO_MPS2 <: AbstractProjMPO @@ -18,14 +18,6 @@ copy(P::ProjMPO_MPS2) = ProjMPO_MPS2(copy(P.PH), copy(P.Ms)) nsite(P::ProjMPO_MPS2) = nsite(P.PH) -function set_nsite!(P::ProjMPO_MPS2, nsite) - set_nsite!(P.PH, nsite) - for m in P.Ms - set_nsite!(m, nsite) - end - return P -end - function makeL!(P::ProjMPO_MPS2, psi::MPS, k::Int) makeL!(P.PH, psi, k) for m in P.Ms diff --git a/src/treetensornetworks/solvers/deprecated/projmps2.jl b/src/treetensornetworks/solvers/deprecated/projmps2.jl index 9f29fdcc..d908943d 100644 --- a/src/treetensornetworks/solvers/deprecated/projmps2.jl +++ b/src/treetensornetworks/solvers/deprecated/projmps2.jl @@ -1,5 +1,4 @@ -import ITensors: - AbstractProjMPO, makeL!, makeR!, set_nsite!, contract, OneITensor, site_range +import ITensors: AbstractProjMPO, makeL!, makeR!, contract, OneITensor, site_range import Base: copy """ @@ -30,11 +29,6 @@ function copy(P::ProjMPS2) return ProjMPS2(P.lpos, P.rpos, P.nsite, copy(P.M), copy(P.LR)) end -function set_nsite!(P::ProjMPS2, nsite) - P.nsite = nsite - return P -end - function makeL!(P::ProjMPS2, psi::MPS, k::Int) # Save the last `L` that is made to help with caching # for DiskProjMPO From cb4ccc19401fa17bbb9020e4333a48b6af3ac810 Mon Sep 17 00:00:00 2001 From: Joey Date: Wed, 27 Mar 2024 14:14:19 -0400 Subject: [PATCH 11/35] Formatting --- src/abstractitensornetwork.jl | 28 ---------------------------- 1 file changed, 28 deletions(-) diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index cf86c020..834a7e40 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -780,34 +780,6 @@ end # Common index checking # -function hascommoninds( - ::typeof(siteinds), A::AbstractITensorNetwork{V}, B::AbstractITensorNetwork{V} -) where {V} - for v in vertices(A) - !hascommoninds(siteinds(A, v), siteinds(B, v)) && return false - end - return true -end - -function check_hascommoninds( - ::typeof(siteinds), A::AbstractITensorNetwork{V}, B::AbstractITensorNetwork{V} -) where {V} - N = nv(A) - if nv(B) ≠ N - throw( - DimensionMismatch( - "$(typeof(A)) and $(typeof(B)) have mismatched number of vertices $N and $(nv(B))." - ), - ) - end - for v in vertices(A) - !hascommoninds(siteinds(A, v), siteinds(B, v)) && error( - "$(typeof(A)) A and $(typeof(B)) B must share site indices. On vertex $v, A has site indices $(siteinds(A, v)) while B has site indices $(siteinds(B, v)).", - ) - end - return nothing -end - function hassameinds( ::typeof(siteinds), A::AbstractITensorNetwork{V}, B::AbstractITensorNetwork{V} ) where {V} From 8efff49689bd17be8812cf395fd7df09b81bbe34 Mon Sep 17 00:00:00 2001 From: Joey Date: Thu, 28 Mar 2024 10:55:23 -0400 Subject: [PATCH 12/35] Modified inner and forms test for new BiLinearForm code --- src/abstractitensornetwork.jl | 53 ++---------------------- src/caches/beliefpropagationcache.jl | 1 + src/contract.jl | 6 +++ src/formnetworks/bilinearformnetwork.jl | 25 +++++++---- src/formnetworks/quadraticformnetwork.jl | 8 ++-- src/imports.jl | 2 - src/inner.jl | 26 +++++++----- test/test_forms.jl | 9 ++-- test/test_inner.jl | 14 +++---- 9 files changed, 58 insertions(+), 86 deletions(-) diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index 834a7e40..2110e9de 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -635,62 +635,15 @@ function split_index( return tn end -function stack_flatten_combine( - tns::Vector{<:AbstractITensorNetwork}; - map_bra_linkinds=prime, - flatten=false, - combine_linkinds=false, -) - tns = copy(tns) - stacked_tn = map_bra_linkinds(popfirst!(tns); sites=[]) - current_suffix = 1 - for tn in tns - stacked_tn_vertices = vertices(stacked_tn) - stacked_tn = disjoint_union(current_suffix => stacked_tn, current_suffix + 1 => tn) - - if flatten - @assert issetequal(vertices(tn), stacked_tn_vertices) - for v in vertices(tn) - stacked_tn = contract( - stacked_tn, (v, current_suffix + 1) => (v, current_suffix); merged_vertex=v - ) - end - else - if current_suffix != 1 - #Strip back - stacked_tn_vertices = [(v, current_suffix) for v in stacked_tn_vertices] - stacked_tn = rename_vertices( - v -> v ∈ stacked_tn_vertices ? first(v) : v, stacked_tn - ) - end - current_suffix += 1 - end - if combine_linkinds - stacked_tn = ITensorNetworks.combine_linkinds(stacked_tn) - end - end - - return stacked_tn -end - -function stack_flatten_combine(tns::AbstractITensorNetwork...; kwargs...) - return stack_flatten_combine([tns...]; kwargs...) -end - -function flatten_networks( - tns::AbstractITensorNetwork...; combine_linkinds=true, map_bra_linkinds=prime -) - return stack_flatten_combine(tns...; flatten=true, combine_linkinds, map_bra_linkinds) -end - +#Just make this call to form network, rip out flatten function inner_network(x::AbstractITensorNetwork, y::AbstractITensorNetwork; kwargs...) - return stack_flatten_combine(dag(x), y; kwargs...) + return BilinearFormNetwork(x, y; kwargs...) end function inner_network( x::AbstractITensorNetwork, A::AbstractITensorNetwork, y::AbstractITensorNetwork; kwargs... ) - return stack_flatten_combine(dag(x), A, y; kwargs...) + return BilinearFormNetwork(x, A, y; kwargs...) end inner_network(x::AbstractITensorNetwork; kwargs...) = inner_network(x, x; kwargs...) diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index a72edf20..da458080 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -11,6 +11,7 @@ end return default_bp_maxiter(undirected_graph(underlying_graph(g))) end default_partitioned_vertices(ψ::AbstractITensorNetwork) = group(v -> v, vertices(ψ)) +default_partitioned_vertices(f::AbstractFormNetwork) = group(v -> original_state_vertex(f, v), vertices(f)) default_cache_update_kwargs(cache) = (; maxiter=20, tol=1e-5) #TODO: Define a version of this that works for QN supporting tensors diff --git a/src/contract.jl b/src/contract.jl index db190ab8..a756c1fb 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -55,6 +55,12 @@ function logscalar( return log(complex(scalar(alg, tn; kwargs...))) end +#This should just pass to a logscalar(bp_cache::BeliefPropagationCache, ...) +#Catch the 0 case in logscalar not in scalar, then pass to exp(logscalar(...)) +#Check if negative before complex() call +#Break down into scalar factors? +#Make general to all algorithms + function logscalar( alg::Algorithm"bp", tn::AbstractITensorNetwork; diff --git a/src/formnetworks/bilinearformnetwork.jl b/src/formnetworks/bilinearformnetwork.jl index 5519c1e3..35128120 100644 --- a/src/formnetworks/bilinearformnetwork.jl +++ b/src/formnetworks/bilinearformnetwork.jl @@ -1,3 +1,6 @@ +default_dual_site_index_map = prime +default_dual_link_index_map = sim + struct BilinearFormNetwork{ V, TensorNetwork<:AbstractITensorNetwork{V}, @@ -18,9 +21,12 @@ function BilinearFormNetwork( operator_vertex_suffix=default_operator_vertex_suffix(), bra_vertex_suffix=default_bra_vertex_suffix(), ket_vertex_suffix=default_ket_vertex_suffix(), + dual_site_index_map = default_dual_site_index_map, + dual_link_index_map = default_dual_link_index_map ) + bra_mapped = dual_link_index_map(dual_site_index_map(bra; links = []); sites = []) tn = disjoint_union( - operator_vertex_suffix => operator, bra_vertex_suffix => bra, ket_vertex_suffix => ket + operator_vertex_suffix => operator, bra_vertex_suffix => dag(bra_mapped), ket_vertex_suffix => ket ) return BilinearFormNetwork( tn, operator_vertex_suffix, bra_vertex_suffix, ket_vertex_suffix @@ -43,24 +49,29 @@ function copy(blf::BilinearFormNetwork) ) end +#Is the ordering of the indices correct here? CHECK THIS +#Put bra into the vector space!!!! function BilinearFormNetwork( - bra::AbstractITensorNetwork, ket::AbstractITensorNetwork; kwargs... + bra::AbstractITensorNetwork, ket::AbstractITensorNetwork; + dual_site_index_map = default_dual_site_index_map, + kwargs... ) - operator_inds = union_all_inds(siteinds(bra), siteinds(ket)) + @assert issetequal(externalinds(bra), externalinds(ket)) + operator_inds = union_all_inds(siteinds(ket), dual_site_index_map(siteinds(ket))) O = delta_network(operator_inds) - return BilinearFormNetwork(O, bra, ket; kwargs...) + return BilinearFormNetwork(O, bra, ket; dual_site_index_map, kwargs...) end function update( - blf::BilinearFormNetwork, original_state_vertex, bra_state::ITensor, ket_state::ITensor + blf::BilinearFormNetwork, original_bra_state_vertex, original_ket_state_vertex, bra_state::ITensor, ket_state::ITensor ) blf = copy(blf) # TODO: Maybe add a check that it really does preserve the graph. setindex_preserve_graph!( - tensornetwork(blf), bra_state, bra_vertex(blf, original_state_vertex) + tensornetwork(blf), bra_state, bra_vertex(blf, original_bra_state_vertex) ) setindex_preserve_graph!( - tensornetwork(blf), ket_state, ket_vertex(blf, original_state_vertex) + tensornetwork(blf), ket_state, ket_vertex(blf, original_ket_state_vertex) ) return blf end diff --git a/src/formnetworks/quadraticformnetwork.jl b/src/formnetworks/quadraticformnetwork.jl index 8aac841a..52485882 100644 --- a/src/formnetworks/quadraticformnetwork.jl +++ b/src/formnetworks/quadraticformnetwork.jl @@ -41,8 +41,7 @@ function QuadraticFormNetwork( dual_inv_index_map=default_inv_index_map, kwargs..., ) - bra = map_inds(dual_index_map, dag(ket)) - blf = BilinearFormNetwork(operator, bra, ket; kwargs...) + blf = BilinearFormNetwork(operator, ket, ket; dual_site_index_map = dual_index_map, dual_link_index_map = dual_index_map, kwargs...) return QuadraticFormNetwork(blf, dual_index_map, dual_inv_index_map) end @@ -52,14 +51,13 @@ function QuadraticFormNetwork( dual_inv_index_map=default_inv_index_map, kwargs..., ) - bra = map_inds(dual_index_map, dag(ket)) - blf = BilinearFormNetwork(bra, ket; kwargs...) + blf = BilinearFormNetwork(bra, ket; dual_site_index_map = dual_index_map, dual_link_index_map = dual_index_map, kwargs...) return QuadraticFormNetwork(blf, dual_index_map, dual_inv_index_map) end function update(qf::QuadraticFormNetwork, original_state_vertex, ket_state::ITensor) state_inds = inds(ket_state) bra_state = replaceinds(dag(ket_state), state_inds, dual_index_map(qf).(state_inds)) - new_blf = update(bilinear_formnetwork(qf), original_state_vertex, bra_state, ket_state) + new_blf = update(bilinear_formnetwork(qf), original_state_vertex, original_state_vertex, bra_state, ket_state) return QuadraticFormNetwork(new_blf, dual_index_map(qf), dual_index_map(qf)) end diff --git a/src/imports.jl b/src/imports.jl index 999201d8..bb4cfade 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -82,8 +82,6 @@ import ITensors: dag, # permute permute, - #commoninds - hascommoninds, # linkdims linkdim, linkdims, diff --git a/src/inner.jl b/src/inner.jl index 7388299d..b3312d22 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -1,6 +1,11 @@ default_inner_partitioned_vertices(tn) = group(v -> first(v), vertices(tn)) +#Default to BP always?! default_algorithm(tns::Vector) = all(is_tree.(tns)) ? "bp" : "exact" +#Default for map_linkinds should be sim. +#Use form code and just default to identity inbetween x and y +#Have ϕ in the same space as y and then a dual_map kwarg? + function inner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; @@ -10,14 +15,15 @@ function inner( return inner(Algorithm(alg), ϕ, ψ; kwargs...) end +#Make [A, ϕ, ψ] a Tuple function inner( ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - alg=default_algorithm([A, ϕ, ψ]), + alg=default_algorithm([ϕ, A, ψ]), kwargs..., ) - return inner(Algorithm(alg), A, ϕ, ψ; kwargs...) + return inner(Algorithm(alg), ϕ, A, ψ; kwargs...) end function inner( @@ -92,10 +98,10 @@ function loginner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; partitioned_verts=default_inner_partitioned_vertices, - map_bra_linkinds=sim, + dual_link_index_map=sim, kwargs..., ) - tn = inner_network(ϕ, ψ; map_bra_linkinds) + tn = inner_network(ϕ, ψ; dual_link_index_map) return logscalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end @@ -105,10 +111,10 @@ function loginner( A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; partitioned_verts=default_inner_partitioned_vertices, - map_bra_linkinds=sim, + dual_link_index_map=sim, kwargs..., ) - tn = inner_network(ϕ, A, ψ; map_bra_linkinds) + tn = inner_network(ϕ, A, ψ; dual_link_index_map) return logscalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end @@ -117,10 +123,10 @@ function inner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; partitioned_verts=default_inner_partitioned_vertices, - map_bra_linkinds=prime, + dual_link_index_map=prime, kwargs..., ) - tn = inner_network(ϕ, ψ; map_bra_linkinds) + tn = inner_network(ϕ, ψ; dual_link_index_map) return scalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end @@ -130,10 +136,10 @@ function inner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; partitioned_verts=default_inner_partitioned_vertices, - map_bra_linkinds=prime, + dual_link_index_map=prime, kwargs..., ) - tn = inner_network(ϕ, A, ψ; map_bra_linkinds) + tn = inner_network(ϕ, A, ψ; dual_link_index_map) return scalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) end diff --git a/test/test_forms.jl b/test/test_forms.jl index 0bfa2d02..a8f76f05 100644 --- a/test/test_forms.jl +++ b/test/test_forms.jl @@ -20,13 +20,12 @@ using SplitApplyCombine @testset "FormNetworks" begin g = named_grid((1, 4)) - s_ket = siteinds("S=1/2", g) - s_bra = prime(s_ket; links=[]) - s_operator = union_all_inds(s_bra, s_ket) + s = siteinds("S=1/2", g) + s_operator = union_all_inds(s, prime(s)) χ, D = 2, 3 Random.seed!(1234) - ψket = randomITensorNetwork(s_ket; link_space=χ) - ψbra = randomITensorNetwork(s_bra; link_space=χ) + ψket = randomITensorNetwork(s; link_space=χ) + ψbra = randomITensorNetwork(s; link_space=χ) A = randomITensorNetwork(s_operator; link_space=D) blf = BilinearFormNetwork(A, ψbra, ψket) diff --git a/test/test_inner.jl b/test/test_inner.jl index 4da47060..5131d59c 100644 --- a/test/test_inner.jl +++ b/test/test_inner.jl @@ -10,7 +10,7 @@ using ITensors: siteinds, dag @testset "Inner products, BP vs exact comparison" begin Random.seed!(1234) - L = 12 + L = 4 χ = 2 g = NamedGraph(Graphs.SimpleGraph(uniform_tree(L))) s = siteinds("S=1/2", g) @@ -18,11 +18,11 @@ using ITensors: siteinds, dag x = randomITensorNetwork(s; link_space=χ) #First lets do it with the flattened version of the network - xy = inner_network(x, y; combine_linkinds=true, flatten=true) + xy = inner_network(x, y) xy_scalar = scalar(xy) - xy_scalar_bp = scalar(xy; alg="bp", partitioned_vertices=group(v -> v, vertices(xy))) + xy_scalar_bp = scalar(xy; alg="bp") xy_scalar_logbp = exp( - logscalar(xy; alg="bp", partitioned_vertices=group(v -> v, vertices(xy))) + logscalar(xy; alg="bp") ) @test xy_scalar ≈ xy_scalar_bp @@ -30,7 +30,7 @@ using ITensors: siteinds, dag @test xy_scalar ≈ xy_scalar_logbp #Now lets do it via the inner function - xy_scalar = inner(x, y; alg="exact", flatten=true, combine_linkinds=true) + xy_scalar = inner(x, y; alg="exact") xy_scalar_bp = inner(x, y; alg="bp") xy_scalar_logbp = exp(loginner(x, y; alg="bp")) @@ -38,9 +38,9 @@ using ITensors: siteinds, dag @test xy_scalar_bp ≈ xy_scalar_logbp @test xy_scalar ≈ xy_scalar_logbp - #test contraction of three layers for expectation values +# #test contraction of three layers for expectation values A = ITensorNetwork(TTN(ITensorNetworks.heisenberg(s), s)) - xAy_scalar = inner(x', A, y; alg="exact", flatten=true, combine_linkinds=true) + xAy_scalar = inner(x', A, y; alg="exact") xAy_scalar_bp = inner(x', A, y; alg="bp") xAy_scalar_logbp = exp(loginner(x', A, y; alg="bp")) From 74e5b0b9579aef2666af9b0c42274c1df1deb252 Mon Sep 17 00:00:00 2001 From: Joey Date: Thu, 28 Mar 2024 12:52:08 -0400 Subject: [PATCH 13/35] Refactor tests to account for changes --- src/abstractitensornetwork.jl | 36 +++++++++++++++++++++++++++++++-- src/gauging.jl | 8 ++++---- src/inner.jl | 28 ++++++++++--------------- test/test_additensornetworks.jl | 30 ++++++--------------------- test/test_apply.jl | 17 ++++++++-------- test/test_inner.jl | 8 ++++---- 6 files changed, 68 insertions(+), 59 deletions(-) diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index aec11022..272ddfef 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -635,7 +635,39 @@ function split_index( return tn end -#Just make this call to form network, rip out flatten +function flatten_networks( + tn1::AbstractITensorNetwork, + tn2::AbstractITensorNetwork; + map_bra_linkinds=sim, + flatten=true, + combine_linkinds=true, + kwargs..., +) + @assert issetequal(vertices(tn1), vertices(tn2)) + tn1 = map_bra_linkinds(tn1; sites=[]) + flattened_net = ⊗(tn1, tn2; kwargs...) + if flatten + for v in vertices(tn1) + flattened_net = contract(flattened_net, (v, 2) => (v, 1); merged_vertex=v) + end + end + if combine_linkinds + flattened_net = ITensorNetworks.combine_linkinds(flattened_net) + end + return flattened_net +end + +function flatten_networks( + tn1::AbstractITensorNetwork, + tn2::AbstractITensorNetwork, + tn3::AbstractITensorNetwork, + tn_tail::AbstractITensorNetwork...; + kwargs..., +) + return flatten_networks(flatten_networks(tn1, tn2; kwargs...), tn3, tn_tail...; kwargs...) +end + + function inner_network(x::AbstractITensorNetwork, y::AbstractITensorNetwork; kwargs...) return BilinearFormNetwork(x, y; kwargs...) end @@ -643,7 +675,7 @@ end function inner_network( x::AbstractITensorNetwork, A::AbstractITensorNetwork, y::AbstractITensorNetwork; kwargs... ) - return BilinearFormNetwork(x, A, y; kwargs...) + return BilinearFormNetwork(A, x, y; kwargs...) end inner_network(x::AbstractITensorNetwork; kwargs...) = inner_network(x, x; kwargs...) diff --git a/src/gauging.jl b/src/gauging.jl index 625b1543..9efd6ed2 100644 --- a/src/gauging.jl +++ b/src/gauging.jl @@ -20,8 +20,8 @@ function copy(ψ::VidalITensorNetwork) end function default_norm_cache(ψ::ITensorNetwork) - ψψ = norm_sqr_network(ψ) - return BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) + ψψ = disjoint_union("bra" => dag(prime(ψ; sites = [])), "ket" => ψ) + return BeliefPropagationCache(ψψ, group(v -> first(v), vertices(ψψ))) end function ITensorNetwork( @@ -46,7 +46,7 @@ function ITensorNetwork( for e in edges(ψ) vsrc, vdst = src(e), dst(e) - pe = partitionedge(bp_cache, (vsrc, 1) => (vdst, 1)) + pe = partitionedge(bp_cache, (vsrc, "bra") => (vdst, "bra")) set!(mts, pe, copy(ITensor[dense(bond_tensor(ψ_vidal, e))])) set!(mts, reverse(pe), copy(ITensor[dense(bond_tensor(ψ_vidal, e))])) end @@ -75,7 +75,7 @@ function vidalitensornetwork_preserve_cache( vsrc, vdst = src(e), dst(e) ψvsrc, ψvdst = ψ_vidal_site_tensors[vsrc], ψ_vidal_site_tensors[vdst] - pe = partitionedge(cache, (vsrc, 1) => (vdst, 1)) + pe = partitionedge(cache, (vsrc, "bra") => (vdst, "bra")) edge_ind = commoninds(ψvsrc, ψvdst) edge_ind_sim = sim(edge_ind) diff --git a/src/inner.jl b/src/inner.jl index b3312d22..487d623e 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -1,26 +1,21 @@ -default_inner_partitioned_vertices(tn) = group(v -> first(v), vertices(tn)) -#Default to BP always?! -default_algorithm(tns::Vector) = all(is_tree.(tns)) ? "bp" : "exact" +default_algorithm(tns::Tuple) = "bp" #Default for map_linkinds should be sim. -#Use form code and just default to identity inbetween x and y -#Have ϕ in the same space as y and then a dual_map kwarg? function inner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - alg=default_algorithm([ϕ, ψ]), + alg=default_algorithm((ϕ, ψ)), kwargs..., ) return inner(Algorithm(alg), ϕ, ψ; kwargs...) end -#Make [A, ϕ, ψ] a Tuple function inner( ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - alg=default_algorithm([ϕ, A, ψ]), + alg=default_algorithm((ϕ, A, ψ)), kwargs..., ) return inner(Algorithm(alg), ϕ, A, ψ; kwargs...) @@ -48,7 +43,6 @@ function inner( ψ::AbstractITensorNetwork; sequence=nothing, contraction_sequence_kwargs=(;), - site_index_map=prime, kwargs..., ) tn = inner_network(ϕ, A, ψ; kwargs...) @@ -74,7 +68,7 @@ function loginner( alg=default_algorithm(ϕ, ψ), kwargs..., ) - return loginner(Algorithm(alg), A, ϕ, ψ; kwargs...) + return loginner(Algorithm(alg), ϕ, A, ψ; kwargs...) end function loginner( @@ -97,7 +91,7 @@ function loginner( alg::Algorithm"bp", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - partitioned_verts=default_inner_partitioned_vertices, + partitioned_verts=default_partitioned_vertices, dual_link_index_map=sim, kwargs..., ) @@ -110,7 +104,7 @@ function loginner( ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - partitioned_verts=default_inner_partitioned_vertices, + partitioned_verts=default_partitioned_vertices, dual_link_index_map=sim, kwargs..., ) @@ -122,8 +116,8 @@ function inner( alg::Algorithm"bp", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - partitioned_verts=default_inner_partitioned_vertices, - dual_link_index_map=prime, + partitioned_verts=default_partitioned_vertices, + dual_link_index_map=sim, kwargs..., ) tn = inner_network(ϕ, ψ; dual_link_index_map) @@ -132,11 +126,11 @@ end function inner( alg::Algorithm"bp", - A::AbstractITensorNetwork, ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - partitioned_verts=default_inner_partitioned_vertices, - dual_link_index_map=prime, + partitioned_verts=default_partitioned_vertices, + dual_link_index_map=sim, kwargs..., ) tn = inner_network(ϕ, A, ψ; dual_link_index_map) diff --git a/test/test_additensornetworks.jl b/test/test_additensornetworks.jl index 2f6c4246..8a53a003 100644 --- a/test/test_additensornetworks.jl +++ b/test/test_additensornetworks.jl @@ -1,5 +1,4 @@ using ITensorNetworks -using ITensorNetworks: inner_network using Test using Compat using ITensors @@ -14,7 +13,7 @@ using Random @testset "add_itensornetworks" begin Random.seed!(5623) - g = named_grid((2, 3)) + g = named_grid((2, 2)) s = siteinds("S=1/2", g) ψ1 = ITensorNetwork(s, v -> "↑") ψ2 = ITensorNetwork(s, v -> "↓") @@ -28,11 +27,9 @@ using Random ψψ_GHZ = inner_network(ψ_GHZ, ψ_GHZ) ψOψ_GHZ = inner_network(ψ_GHZ, Oψ_GHZ) - @test ITensors.contract(ψOψ_GHZ)[] / ITensors.contract(ψψ_GHZ)[] == 0.0 + @test scalar(ψOψ_GHZ) / scalar(ψψ_GHZ) == 0.0 χ = 3 - g = hexagonal_lattice_graph(1, 2) - s1 = siteinds("S=1/2", g) s2 = copy(s1) rem_edge!(s2, NamedEdge((1, 1) => (1, 2))) @@ -52,26 +49,11 @@ using Random Oψ2 = copy(ψ2) Oψ2[v] = apply(op("Sz", s2[v]), Oψ2[v]) - ψψ_12 = inner_network(ψ12, ψ12) - ψOψ_12 = inner_network(ψ12, Oψ12) - - ψ1ψ2 = inner_network(ψ1, ψ2) - ψ1Oψ2 = inner_network(ψ1, Oψ2) - - ψψ_2 = inner_network(ψ2, ψ2) - ψOψ_2 = inner_network(ψ2, Oψ2) - - ψψ_1 = inner_network(ψ1, ψ1) - ψOψ_1 = inner_network(ψ1, Oψ1) - + alg = "exact" expec_method1 = - ( - ITensors.contract(ψOψ_1)[] + - ITensors.contract(ψOψ_2)[] + - 2 * ITensors.contract(ψ1Oψ2)[] - ) / - (ITensors.contract(ψψ_1)[] + ITensors.contract(ψψ_2)[] + 2 * ITensors.contract(ψ1ψ2)[]) - expec_method2 = ITensors.contract(ψOψ_12)[] / ITensors.contract(ψψ_12)[] + (inner(ψ1, Oψ1; alg) + inner(ψ2, Oψ2; alg) + 2*inner(ψ1, Oψ2; alg)) / + (norm_sqr(ψ1; alg) + norm_sqr(ψ2; alg) + 2 * inner(ψ1, ψ2; alg)) + expec_method2 = inner(ψ12, Oψ12; alg) / norm_sqr(ψ12; alg) @test expec_method1 ≈ expec_method2 end diff --git a/test/test_apply.jl b/test/test_apply.jl index 3c0644a0..d6ca547f 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,6 +1,6 @@ using ITensorNetworks using ITensorNetworks: - environment, update, inner, norm_sqr_network, BeliefPropagationCache, VidalITensorNetwork + environment, update, inner, BeliefPropagationCache, VidalITensorNetwork, norm_sqr_network_fast using Test using Compat using ITensors @@ -19,8 +19,7 @@ using SplitApplyCombine χ = 2 ψ = randomITensorNetwork(s; link_space=χ) v1, v2 = (2, 2), (1, 2) - ψψ = norm_sqr_network(ψ) - + ψψ = disjoint_union("bra" => dag(prime(ψ; sites = [])), "ket" => ψ) #Simple Belief Propagation Grouping bp_cache = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) bp_cache = update(bp_cache; maxiter=20) @@ -31,7 +30,9 @@ using SplitApplyCombine #This grouping will correspond to calculating the environments exactly (each column of the grid is a partition) bp_cache = BeliefPropagationCache(ψψ, group(v -> v[1][1], vertices(ψψ))) bp_cache = update(bp_cache; maxiter=20) - envsGBP = environment(bp_cache, [(v1, 1), (v1, 2), (v2, 1), (v2, 2)]) + envsGBP = environment(bp_cache, [(v1, "bra"), (v1, "ket"), (v2, "bra"), (v2, "ket")]) + + inner_alg = "exact" ngates = 5 @@ -59,11 +60,11 @@ using SplitApplyCombine print_fidelity_loss=true, envisposdef=true, ) - fSBP = inner(ψOSBP, ψOexact) / sqrt(inner(ψOexact, ψOexact) * inner(ψOSBP, ψOSBP)) + fSBP = inner(ψOSBP, ψOexact; alg = inner_alg) / sqrt(inner(ψOexact, ψOexact; alg = inner_alg) * inner(ψOSBP, ψOSBP; alg = inner_alg)) fVidal = - inner(ψOVidal_symm, ψOexact) / - sqrt(inner(ψOexact, ψOexact) * inner(ψOVidal_symm, ψOVidal_symm)) - fGBP = inner(ψOGBP, ψOexact) / sqrt(inner(ψOexact, ψOexact) * inner(ψOGBP, ψOGBP)) + inner(ψOVidal_symm, ψOexact; alg = inner_alg) / + sqrt(inner(ψOexact, ψOexact; alg = inner_alg) * inner(ψOVidal_symm, ψOVidal_symm; alg = inner_alg)) + fGBP = inner(ψOGBP, ψOexact; alg = inner_alg) / sqrt(inner(ψOexact, ψOexact; alg = inner_alg) * inner(ψOGBP, ψOGBP; alg = inner_alg)) @test real(fGBP * conj(fGBP)) >= real(fSBP * conj(fSBP)) diff --git a/test/test_inner.jl b/test/test_inner.jl index 5131d59c..e8ff687e 100644 --- a/test/test_inner.jl +++ b/test/test_inner.jl @@ -38,11 +38,11 @@ using ITensors: siteinds, dag @test xy_scalar_bp ≈ xy_scalar_logbp @test xy_scalar ≈ xy_scalar_logbp -# #test contraction of three layers for expectation values + #test contraction of three layers for expectation values A = ITensorNetwork(TTN(ITensorNetworks.heisenberg(s), s)) - xAy_scalar = inner(x', A, y; alg="exact") - xAy_scalar_bp = inner(x', A, y; alg="bp") - xAy_scalar_logbp = exp(loginner(x', A, y; alg="bp")) + xAy_scalar = inner(x, A, y; alg="exact") + xAy_scalar_bp = inner(x, A, y; alg="bp") + xAy_scalar_logbp = exp(loginner(x, A, y; alg="bp")) @test xAy_scalar ≈ xAy_scalar_bp @test xAy_scalar_bp ≈ xAy_scalar_logbp From 4e8a1ef013c64cb7af0e9c06dcd174aad80384a7 Mon Sep 17 00:00:00 2001 From: Joey Date: Thu, 28 Mar 2024 14:52:07 -0400 Subject: [PATCH 14/35] All tests refactored. logscalar in terms of scalarnorm --- examples/dynamics/2d_ising_imag_tebd.jl | 11 ++--- examples/mincut.jl | 2 - src/abstractitensornetwork.jl | 5 ++- src/caches/beliefpropagationcache.jl | 18 ++++++-- src/contract.jl | 45 ++++++------------- src/expect.jl | 14 +++--- src/formnetworks/bilinearformnetwork.jl | 23 ++++++---- src/formnetworks/quadraticformnetwork.jl | 25 +++++++++-- src/gauging.jl | 2 +- src/inner.jl | 12 ++--- .../solvers/deprecated/projmpo_apply.jl | 1 + .../solvers/deprecated/projmpo_mps2.jl | 1 + .../solvers/deprecated/projmps2.jl | 1 + test/test_additensornetworks.jl | 4 +- test/test_apply.jl | 23 +++++++--- test/test_contraction_sequence.jl | 2 + test/test_inner.jl | 4 +- test/test_itensornetwork.jl | 26 +++++------ test/test_tebd.jl | 15 +++---- 19 files changed, 128 insertions(+), 106 deletions(-) diff --git a/examples/dynamics/2d_ising_imag_tebd.jl b/examples/dynamics/2d_ising_imag_tebd.jl index fe917523..8fa45a19 100644 --- a/examples/dynamics/2d_ising_imag_tebd.jl +++ b/examples/dynamics/2d_ising_imag_tebd.jl @@ -8,7 +8,7 @@ Random.seed!(1234) ITensors.disable_warn_order() -system_dims = (6, 6) +system_dims = (2, 3) n = prod(system_dims) g = named_grid(system_dims) @@ -58,16 +58,13 @@ println("maxdim = $χ") @show β, Δβ @show ortho -# Contraction sequence for exactly computing expectation values -inner_sequence = reduce((x, y) -> [x, y], vec(Tuple.(CartesianIndices(system_dims)))) - println("\nFirst run TEBD without orthogonalization") ψ = @time tebd( group_terms(ℋ, g), ψ_init; β, Δβ, cutoff=1e-8, maxdim=χ, ortho=false, print_frequency=1 ) println("\nMeasure energy expectation value") -E = @time expect(ℋ, ψ; sequence=inner_sequence) +E = @time expect(ℋ, ψ) @show E println("\nThen run TEBD with orthogonalization (more accurate)") @@ -76,11 +73,11 @@ println("\nThen run TEBD with orthogonalization (more accurate)") ) println("\nMeasure energy expectation value") -E = @time expect(ℋ, ψ; sequence=inner_sequence) +E = @time expect(ℋ, ψ) @show E println("\nMeasure magnetization") -Z_dict = @time expect("Z", ψ; sequence=inner_sequence) +Z_dict = @time expect("Z", ψ) Z = [Z_dict[Tuple(I)] for I in CartesianIndices(system_dims)] display(Z) display(heatmap(Z)) diff --git a/examples/mincut.jl b/examples/mincut.jl index 298157cc..ba3c0211 100644 --- a/examples/mincut.jl +++ b/examples/mincut.jl @@ -8,8 +8,6 @@ s = siteinds("S=1/2", g) ψ = ITensorNetwork(s; link_space=10) -# ρ = flatten_networks(dag(ψ), ψ') - # Or: ss = ∪(dag(s), s'; merge_data=union) diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index 272ddfef..41ba3722 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -667,7 +667,6 @@ function flatten_networks( return flatten_networks(flatten_networks(tn1, tn2; kwargs...), tn3, tn_tail...; kwargs...) end - function inner_network(x::AbstractITensorNetwork, y::AbstractITensorNetwork; kwargs...) return BilinearFormNetwork(x, y; kwargs...) end @@ -679,7 +678,9 @@ function inner_network( end inner_network(x::AbstractITensorNetwork; kwargs...) = inner_network(x, x; kwargs...) -const norm_sqr_network = inner_network +function norm_sqr_network(ψ::AbstractITensorNetwork) + return disjoint_union("bra" => dag(prime(ψ; sites=[])), "ket" => ψ) +end #Ideally this will not be necessary but this is a temporary fast version to avoid the overhead of `disjoint_union` function norm_sqr_network_fast(ψ::AbstractITensorNetwork) diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index da458080..fe0d7c27 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -11,7 +11,9 @@ end return default_bp_maxiter(undirected_graph(underlying_graph(g))) end default_partitioned_vertices(ψ::AbstractITensorNetwork) = group(v -> v, vertices(ψ)) -default_partitioned_vertices(f::AbstractFormNetwork) = group(v -> original_state_vertex(f, v), vertices(f)) +function default_partitioned_vertices(f::AbstractFormNetwork) + return group(v -> original_state_vertex(f, v), vertices(f)) +end default_cache_update_kwargs(cache) = (; maxiter=20, tol=1e-5) #TODO: Define a version of this that works for QN supporting tensors @@ -40,8 +42,14 @@ function BeliefPropagationCache(tn, partitioned_vertices; kwargs...) return BeliefPropagationCache(ptn; kwargs...) end -function BeliefPropagationCache(tn; kwargs...) - return BeliefPropagationCache(tn, default_partitioning(tn); kwargs...) +function BeliefPropagationCache( + tn; partitioned_vertices=default_partitioned_vertices(tn), kwargs... +) + return BeliefPropagationCache(tn, partitioned_vertices; kwargs...) +end + +function cache(alg::Algorithm"bp", tn; kwargs...) + return BeliefPropagationCache(tn; kwargs...) end function partitioned_itensornetwork(bp_cache::BeliefPropagationCache) @@ -270,3 +278,7 @@ end function edge_norms(bp_cache::BeliefPropagationCache) return edge_norms(bp_cache, partitionedges(partitioned_itensornetwork(bp_cache))) end + +function scalar_factors(bp_cache::BeliefPropagationCache) + return vertex_norms(bp_cache), edge_norms(bp_cache) +end diff --git a/src/contract.jl b/src/contract.jl index a756c1fb..a7fa4ec3 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -55,50 +55,33 @@ function logscalar( return log(complex(scalar(alg, tn; kwargs...))) end -#This should just pass to a logscalar(bp_cache::BeliefPropagationCache, ...) -#Catch the 0 case in logscalar not in scalar, then pass to exp(logscalar(...)) -#Check if negative before complex() call -#Break down into scalar factors? -#Make general to all algorithms - function logscalar( - alg::Algorithm"bp", + alg::Algorithm, tn::AbstractITensorNetwork; (cache!)=nothing, - partitioned_vertices=default_partitioned_vertices(tn), + cache_construction_kwargs=(;), update_cache=isnothing(cache!), cache_update_kwargs=default_cache_update_kwargs(cache!), ) if isnothing(cache!) - cache! = Ref(BeliefPropagationCache(tn, partitioned_vertices)) + cache! = Ref(cache(alg, tn; cache_construction_kwargs...)) end if update_cache cache![] = update(cache![]; cache_update_kwargs...) end - numerator_terms, denominator_terms = vertex_norms(cache![]), edge_norms(cache![]) - - return sum(log.(complex.(numerator_terms))) - sum(log.(complex.((denominator_terms)))) -end - -function scalar( - alg::Algorithm"bp", - tn::AbstractITensorNetwork; - (cache!)=nothing, - partitioned_vertices=default_partitioned_vertices(tn), - update_cache=isnothing(cache!), - cache_update_kwargs=default_cache_update_kwargs(cache!), -) - if isnothing(cache!) - cache! = Ref(BeliefPropagationCache(tn, partitioned_vertices)) - end - - if update_cache - cache![] = update(cache![]; cache_update_kwargs...) + numerator_terms, denominator_terms = scalar_factors(cache![]) + terms = vcat(numerator_terms, denominator_terms) + if any(==(0), terms) + return -Inf + elseif all(>=(0), terms) + return sum(log.(numerator_terms)) - sum(log.((denominator_terms))) + else + return sum(log.(complex.(numerator_terms))) - sum(log.(complex.((denominator_terms)))) end +end - numerator_terms, denominator_terms = vertex_norms(cache![]), edge_norms(cache![]) - - return prod(numerator_terms) / prod(denominator_terms) +function scalar(alg::Algorithm"bp", tn::AbstractITensorNetwork; kwargs...) + return exp(logscalar(alg, tn; kwargs...)) end diff --git a/src/expect.jl b/src/expect.jl index fd939e42..4060514a 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -5,7 +5,6 @@ function expect( maxdim=nothing, ortho=false, sequence=nothing, - flatten=true, vertices=vertices(ψ), ) s = siteinds(ψ) @@ -13,13 +12,13 @@ function expect( # ElT = ishermitian(ITensors.op(op, s[vertices[1]])) ? real(ElT) : ElT res = Dictionary(vertices, Vector{ElT}(undef, length(vertices))) if isnothing(sequence) - sequence = contraction_sequence(inner_network(ψ, ψ; flatten)) + sequence = contraction_sequence(inner_network(ψ, ψ)) end - normψ² = norm_sqr(ψ; alg="exact", sequence, flatten) + normψ² = norm_sqr(ψ; alg="exact", sequence) for v in vertices O = ITensor(Op(op, v), s) Oψ = apply(O, ψ; cutoff, maxdim, ortho) - res[v] = inner(ψ, Oψ; alg="exact", sequence, flatten) / normψ² + res[v] = inner(ψ, Oψ; alg="exact", sequence) / normψ² end return res end @@ -31,17 +30,16 @@ function expect( maxdim=nothing, ortho=false, sequence=nothing, - flatten=true, ) s = siteinds(ψ) # h⃗ = Vector{ITensor}(ℋ, s) if isnothing(sequence) - sequence = contraction_sequence(inner_network(ψ, ψ; flatten)) + sequence = contraction_sequence(inner_network(ψ, ψ)) end h⃗ψ = [apply(hᵢ, ψ; cutoff, maxdim, ortho) for hᵢ in ITensors.terms(ℋ)] - ψhᵢψ = [inner(ψ, hᵢψ; alg="exact", flatten, sequence) for hᵢψ in h⃗ψ] + ψhᵢψ = [inner(ψ, hᵢψ; alg="exact", sequence) for hᵢψ in h⃗ψ] ψh⃗ψ = sum(ψhᵢψ) - ψψ = norm_sqr(ψ; alg="exact", sequence, flatten) + ψψ = norm_sqr(ψ; alg="exact", sequence) return ψh⃗ψ / ψψ end diff --git a/src/formnetworks/bilinearformnetwork.jl b/src/formnetworks/bilinearformnetwork.jl index 35128120..689b6dbc 100644 --- a/src/formnetworks/bilinearformnetwork.jl +++ b/src/formnetworks/bilinearformnetwork.jl @@ -21,12 +21,14 @@ function BilinearFormNetwork( operator_vertex_suffix=default_operator_vertex_suffix(), bra_vertex_suffix=default_bra_vertex_suffix(), ket_vertex_suffix=default_ket_vertex_suffix(), - dual_site_index_map = default_dual_site_index_map, - dual_link_index_map = default_dual_link_index_map + dual_site_index_map=default_dual_site_index_map, + dual_link_index_map=default_dual_link_index_map, ) - bra_mapped = dual_link_index_map(dual_site_index_map(bra; links = []); sites = []) + bra_mapped = dual_link_index_map(dual_site_index_map(bra; links=[]); sites=[]) tn = disjoint_union( - operator_vertex_suffix => operator, bra_vertex_suffix => dag(bra_mapped), ket_vertex_suffix => ket + operator_vertex_suffix => operator, + bra_vertex_suffix => dag(bra_mapped), + ket_vertex_suffix => ket, ) return BilinearFormNetwork( tn, operator_vertex_suffix, bra_vertex_suffix, ket_vertex_suffix @@ -52,9 +54,10 @@ end #Is the ordering of the indices correct here? CHECK THIS #Put bra into the vector space!!!! function BilinearFormNetwork( - bra::AbstractITensorNetwork, ket::AbstractITensorNetwork; - dual_site_index_map = default_dual_site_index_map, - kwargs... + bra::AbstractITensorNetwork, + ket::AbstractITensorNetwork; + dual_site_index_map=default_dual_site_index_map, + kwargs..., ) @assert issetequal(externalinds(bra), externalinds(ket)) operator_inds = union_all_inds(siteinds(ket), dual_site_index_map(siteinds(ket))) @@ -63,7 +66,11 @@ function BilinearFormNetwork( end function update( - blf::BilinearFormNetwork, original_bra_state_vertex, original_ket_state_vertex, bra_state::ITensor, ket_state::ITensor + blf::BilinearFormNetwork, + original_bra_state_vertex, + original_ket_state_vertex, + bra_state::ITensor, + ket_state::ITensor, ) blf = copy(blf) # TODO: Maybe add a check that it really does preserve the graph. diff --git a/src/formnetworks/quadraticformnetwork.jl b/src/formnetworks/quadraticformnetwork.jl index 52485882..8eaa6246 100644 --- a/src/formnetworks/quadraticformnetwork.jl +++ b/src/formnetworks/quadraticformnetwork.jl @@ -41,7 +41,14 @@ function QuadraticFormNetwork( dual_inv_index_map=default_inv_index_map, kwargs..., ) - blf = BilinearFormNetwork(operator, ket, ket; dual_site_index_map = dual_index_map, dual_link_index_map = dual_index_map, kwargs...) + blf = BilinearFormNetwork( + operator, + ket, + ket; + dual_site_index_map=dual_index_map, + dual_link_index_map=dual_index_map, + kwargs..., + ) return QuadraticFormNetwork(blf, dual_index_map, dual_inv_index_map) end @@ -51,13 +58,25 @@ function QuadraticFormNetwork( dual_inv_index_map=default_inv_index_map, kwargs..., ) - blf = BilinearFormNetwork(bra, ket; dual_site_index_map = dual_index_map, dual_link_index_map = dual_index_map, kwargs...) + blf = BilinearFormNetwork( + bra, + ket; + dual_site_index_map=dual_index_map, + dual_link_index_map=dual_index_map, + kwargs..., + ) return QuadraticFormNetwork(blf, dual_index_map, dual_inv_index_map) end function update(qf::QuadraticFormNetwork, original_state_vertex, ket_state::ITensor) state_inds = inds(ket_state) bra_state = replaceinds(dag(ket_state), state_inds, dual_index_map(qf).(state_inds)) - new_blf = update(bilinear_formnetwork(qf), original_state_vertex, original_state_vertex, bra_state, ket_state) + new_blf = update( + bilinear_formnetwork(qf), + original_state_vertex, + original_state_vertex, + bra_state, + ket_state, + ) return QuadraticFormNetwork(new_blf, dual_index_map(qf), dual_index_map(qf)) end diff --git a/src/gauging.jl b/src/gauging.jl index 9efd6ed2..31cb6731 100644 --- a/src/gauging.jl +++ b/src/gauging.jl @@ -20,7 +20,7 @@ function copy(ψ::VidalITensorNetwork) end function default_norm_cache(ψ::ITensorNetwork) - ψψ = disjoint_union("bra" => dag(prime(ψ; sites = [])), "ket" => ψ) + ψψ = norm_sqr_network(ψ) return BeliefPropagationCache(ψψ, group(v -> first(v), vertices(ψψ))) end diff --git a/src/inner.jl b/src/inner.jl index 487d623e..25e79894 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -91,12 +91,11 @@ function loginner( alg::Algorithm"bp", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - partitioned_verts=default_partitioned_vertices, dual_link_index_map=sim, kwargs..., ) tn = inner_network(ϕ, ψ; dual_link_index_map) - return logscalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) + return logscalar(alg, tn; kwargs...) end function loginner( @@ -104,24 +103,22 @@ function loginner( ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - partitioned_verts=default_partitioned_vertices, dual_link_index_map=sim, kwargs..., ) tn = inner_network(ϕ, A, ψ; dual_link_index_map) - return logscalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) + return logscalar(alg, tn; kwargs...) end function inner( alg::Algorithm"bp", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - partitioned_verts=default_partitioned_vertices, dual_link_index_map=sim, kwargs..., ) tn = inner_network(ϕ, ψ; dual_link_index_map) - return scalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) + return scalar(alg, tn; kwargs...) end function inner( @@ -129,12 +126,11 @@ function inner( ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - partitioned_verts=default_partitioned_vertices, dual_link_index_map=sim, kwargs..., ) tn = inner_network(ϕ, A, ψ; dual_link_index_map) - return scalar(alg, tn; partitioned_vertices=partitioned_verts(tn), kwargs...) + return scalar(alg, tn; kwargs...) end # TODO: rename `sqnorm` to match https://github.com/JuliaStats/Distances.jl, diff --git a/src/treetensornetworks/solvers/deprecated/projmpo_apply.jl b/src/treetensornetworks/solvers/deprecated/projmpo_apply.jl index e69de29b..8b137891 100644 --- a/src/treetensornetworks/solvers/deprecated/projmpo_apply.jl +++ b/src/treetensornetworks/solvers/deprecated/projmpo_apply.jl @@ -0,0 +1 @@ + diff --git a/src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl b/src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl index e69de29b..8b137891 100644 --- a/src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl +++ b/src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl @@ -0,0 +1 @@ + diff --git a/src/treetensornetworks/solvers/deprecated/projmps2.jl b/src/treetensornetworks/solvers/deprecated/projmps2.jl index e69de29b..8b137891 100644 --- a/src/treetensornetworks/solvers/deprecated/projmps2.jl +++ b/src/treetensornetworks/solvers/deprecated/projmps2.jl @@ -0,0 +1 @@ + diff --git a/test/test_additensornetworks.jl b/test/test_additensornetworks.jl index 8a53a003..bf92fbbf 100644 --- a/test/test_additensornetworks.jl +++ b/test/test_additensornetworks.jl @@ -51,8 +51,8 @@ using Random alg = "exact" expec_method1 = - (inner(ψ1, Oψ1; alg) + inner(ψ2, Oψ2; alg) + 2*inner(ψ1, Oψ2; alg)) / - (norm_sqr(ψ1; alg) + norm_sqr(ψ2; alg) + 2 * inner(ψ1, ψ2; alg)) + (inner(ψ1, Oψ1; alg) + inner(ψ2, Oψ2; alg) + 2 * inner(ψ1, Oψ2; alg)) / + (norm_sqr(ψ1; alg) + norm_sqr(ψ2; alg) + 2 * inner(ψ1, ψ2; alg)) expec_method2 = inner(ψ12, Oψ12; alg) / norm_sqr(ψ12; alg) @test expec_method1 ≈ expec_method2 diff --git a/test/test_apply.jl b/test/test_apply.jl index d6ca547f..e247f727 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,6 +1,11 @@ using ITensorNetworks using ITensorNetworks: - environment, update, inner, BeliefPropagationCache, VidalITensorNetwork, norm_sqr_network_fast + environment, + update, + inner, + BeliefPropagationCache, + VidalITensorNetwork, + norm_sqr_network_fast using Test using Compat using ITensors @@ -19,7 +24,7 @@ using SplitApplyCombine χ = 2 ψ = randomITensorNetwork(s; link_space=χ) v1, v2 = (2, 2), (1, 2) - ψψ = disjoint_union("bra" => dag(prime(ψ; sites = [])), "ket" => ψ) + ψψ = disjoint_union("bra" => dag(prime(ψ; sites=[])), "ket" => ψ) #Simple Belief Propagation Grouping bp_cache = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) bp_cache = update(bp_cache; maxiter=20) @@ -60,11 +65,17 @@ using SplitApplyCombine print_fidelity_loss=true, envisposdef=true, ) - fSBP = inner(ψOSBP, ψOexact; alg = inner_alg) / sqrt(inner(ψOexact, ψOexact; alg = inner_alg) * inner(ψOSBP, ψOSBP; alg = inner_alg)) + fSBP = + inner(ψOSBP, ψOexact; alg=inner_alg) / + sqrt(inner(ψOexact, ψOexact; alg=inner_alg) * inner(ψOSBP, ψOSBP; alg=inner_alg)) fVidal = - inner(ψOVidal_symm, ψOexact; alg = inner_alg) / - sqrt(inner(ψOexact, ψOexact; alg = inner_alg) * inner(ψOVidal_symm, ψOVidal_symm; alg = inner_alg)) - fGBP = inner(ψOGBP, ψOexact; alg = inner_alg) / sqrt(inner(ψOexact, ψOexact; alg = inner_alg) * inner(ψOGBP, ψOGBP; alg = inner_alg)) + inner(ψOVidal_symm, ψOexact; alg=inner_alg) / sqrt( + inner(ψOexact, ψOexact; alg=inner_alg) * + inner(ψOVidal_symm, ψOVidal_symm; alg=inner_alg), + ) + fGBP = + inner(ψOGBP, ψOexact; alg=inner_alg) / + sqrt(inner(ψOexact, ψOexact; alg=inner_alg) * inner(ψOGBP, ψOGBP; alg=inner_alg)) @test real(fGBP * conj(fGBP)) >= real(fSBP * conj(fSBP)) diff --git a/test/test_contraction_sequence.jl b/test/test_contraction_sequence.jl index 5470dabe..aae6bd0e 100644 --- a/test/test_contraction_sequence.jl +++ b/test/test_contraction_sequence.jl @@ -5,6 +5,8 @@ using Random using Test using EinExprs: Exhaustive, Greedy, HyPar +using ITensorNetworks: norm_sqr_network + Random.seed!(1234) ITensors.disable_warn_order() diff --git a/test/test_inner.jl b/test/test_inner.jl index e8ff687e..0b90a3cd 100644 --- a/test/test_inner.jl +++ b/test/test_inner.jl @@ -21,9 +21,7 @@ using ITensors: siteinds, dag xy = inner_network(x, y) xy_scalar = scalar(xy) xy_scalar_bp = scalar(xy; alg="bp") - xy_scalar_logbp = exp( - logscalar(xy; alg="bp") - ) + xy_scalar_logbp = exp(logscalar(xy; alg="bp")) @test xy_scalar ≈ xy_scalar_bp @test xy_scalar_bp ≈ xy_scalar_logbp diff --git a/test/test_itensornetwork.jl b/test/test_itensornetwork.jl index b1e6c097..d08ffa7e 100644 --- a/test/test_itensornetwork.jl +++ b/test/test_itensornetwork.jl @@ -88,10 +88,10 @@ using ITensorNetworks: norm_sqr, neighbor_itensors g = named_grid(dims) s = siteinds("S=1/2", g) ψ = ITensorNetwork(s, v -> "↑") - tn = inner_network(ψ, ψ) - tn_2 = contract(tn, ((1, 2), 2) => ((1, 2), 1)) - @test !has_vertex(tn_2, ((1, 2), 2)) - @test tn_2[((1, 2), 1)] ≈ tn[((1, 2), 2)] * tn[((1, 2), 1)] + tn = norm_sqr_network(ψ) + tn_2 = contract(tn, ((1, 2), "ket") => ((1, 2), "bra")) + @test !has_vertex(tn_2, ((1, 2), "ket")) + @test tn_2[((1, 2), "bra")] ≈ tn[((1, 2), "ket")] * tn[((1, 2), "bra")] end @testset "Remove edge (regression test for issue #5)" begin @@ -100,15 +100,15 @@ using ITensorNetworks: norm_sqr, neighbor_itensors s = siteinds("S=1/2", g) ψ = ITensorNetwork(s, v -> "↑") rem_vertex!(ψ, (1, 2)) - tn = inner_network(ψ, ψ) - @test !has_vertex(tn, ((1, 2), 1)) - @test !has_vertex(tn, ((1, 2), 2)) - @test has_vertex(tn, ((1, 1), 1)) - @test has_vertex(tn, ((1, 1), 2)) - @test has_vertex(tn, ((2, 1), 1)) - @test has_vertex(tn, ((2, 1), 2)) - @test has_vertex(tn, ((2, 2), 1)) - @test has_vertex(tn, ((2, 2), 2)) + tn = norm_sqr_network(ψ) + @test !has_vertex(tn, ((1, 2), "bra")) + @test !has_vertex(tn, ((1, 2), "ket")) + @test has_vertex(tn, ((1, 1), "bra")) + @test has_vertex(tn, ((1, 1), "ket")) + @test has_vertex(tn, ((2, 1), "bra")) + @test has_vertex(tn, ((2, 1), "ket")) + @test has_vertex(tn, ((2, 2), "bra")) + @test has_vertex(tn, ((2, 2), "ket")) end @testset "Custom element type" for eltype in (Float32, Float64, ComplexF32, ComplexF64), diff --git a/test/test_tebd.jl b/test/test_tebd.jl index 3a9e0ecc..3600b3ba 100644 --- a/test/test_tebd.jl +++ b/test/test_tebd.jl @@ -5,7 +5,7 @@ using Test ITensors.disable_warn_order() @testset "Ising TEBD" begin - dims = (4, 4) + dims = (2, 3) n = prod(dims) g = named_grid(dims) @@ -33,11 +33,8 @@ ITensors.disable_warn_order() β = 2.0 Δβ = 0.2 - # Sequence for contracting expectation values - inner_sequence = reduce((x, y) -> [x, y], vec(Tuple.(CartesianIndices(dims)))) - ψ_init = ITensorNetwork(s, v -> "↑") - E0 = expect(ℋ, ψ_init; sequence=inner_sequence) + E0 = expect(ℋ, ψ_init) ψ = tebd( group_terms(ℋ, g), ψ_init; @@ -48,7 +45,7 @@ ITensors.disable_warn_order() ortho=false, print_frequency=typemax(Int), ) - E1 = expect(ℋ, ψ; sequence=inner_sequence) + E1 = expect(ℋ, ψ) ψ = tebd( group_terms(ℋ, g), ψ_init; @@ -59,8 +56,8 @@ ITensors.disable_warn_order() ortho=true, print_frequency=typemax(Int), ) - E2 = expect(ℋ, ψ; sequence=inner_sequence) + E2 = expect(ℋ, ψ) @show E0, E1, E2, E_dmrg - @test (((abs((E2 - E1) / E2) < 1e-4) && (E1 < E0)) || (E2 < E1 < E0)) - @test E2 ≈ E_dmrg rtol = 1e-4 + @test (((abs((E2 - E1) / E2) < 1e-3) && (E1 < E0)) || (E2 < E1 < E0)) + @test E2 ≈ E_dmrg rtol = 1e-3 end From ba77b9e6c7546847b8ab84af4d35692bcf1b0ca3 Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 29 Mar 2024 08:48:28 -0400 Subject: [PATCH 15/35] Refactor test bp and test inner for new namespace formatting --- src/boundarymps.jl | 2 +- src/contract.jl | 4 ++++ src/inner.jl | 2 +- test/test_belief_propagation.jl | 32 ++++++++++++++++++++------------ test/test_inner.jl | 12 ++++++------ 5 files changed, 32 insertions(+), 20 deletions(-) diff --git a/src/boundarymps.jl b/src/boundarymps.jl index 8e27868e..7afd5ac7 100644 --- a/src/boundarymps.jl +++ b/src/boundarymps.jl @@ -9,5 +9,5 @@ function contract_boundary_mps(tn::ITensorNetwork; kwargs...) end T = MPO([tn[i1, d2 - 1] for i1 in 1:d1]) vR = MPS([tn[i1, d2] for i1 in 1:d1]) - return inner(dag(vL), T, vR)[] + return ITensors.inner(dag(vL), T, vR)[] end diff --git a/src/contract.jl b/src/contract.jl index 75899200..2b6ee6ae 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -15,6 +15,10 @@ function NDTensors.contract( return contract(Vector{ITensor}(tn); sequence=sequence_linear_index, kwargs...) end +function NDTensors.contract(alg::Algorithm"exact", tensors::Vector{ITensor}; kwargs...) + return contract(tensors; kwargs...) +end + function NDTensors.contract( alg::Union{Algorithm"density_matrix",Algorithm"ttn_svd"}, tn::AbstractITensorNetwork; diff --git a/src/inner.jl b/src/inner.jl index 25e79894..e165551a 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -137,6 +137,6 @@ end # or `norm_sqr` to match `LinearAlgebra.norm_sqr` norm_sqr(ψ::AbstractITensorNetwork; kwargs...) = inner(ψ, ψ; kwargs...) -function norm(ψ::AbstractITensorNetwork; kwargs...) +function LinearAlgebra.norm(ψ::AbstractITensorNetwork; kwargs...) return sqrt(abs(real(norm_sqr(ψ; kwargs...)))) end diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index 213257e2..403ada65 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -1,8 +1,8 @@ +using Test using ITensorNetworks using ITensorNetworks: ising_network, split_index, - inner, contract_boundary_mps, BeliefPropagationCache, tensornetwork, @@ -10,16 +10,24 @@ using ITensorNetworks: update_factor, environment, contract -using Test -using Compat -using ITensors -using LinearAlgebra -using NamedGraphs -using SplitApplyCombine -using Random -using Metis -ITensors.disable_warn_order() +using LinearAlgebra: eigvals, tr +using NamedGraphs: named_grid, named_comb_tree, PartitionVertex, NamedEdge +using Random: seed! +using SplitApplyCombine: group +using ITensors: + siteinds, + apply, + randomITensor, + disable_warn_order, + prime, + op, + ITensor, + inds, + combiner, + array + +disable_warn_order() @testset "belief_propagation" begin @@ -28,7 +36,7 @@ ITensors.disable_warn_order() g = named_grid(g_dims) s = siteinds("S=1/2", g) χ = 4 - Random.seed!(1234) + seed!(1234) ψ = randomITensorNetwork(s; link_space=χ) ψψ = ψ ⊗ prime(dag(ψ); sites=[]) @@ -58,7 +66,7 @@ ITensors.disable_warn_order() g = named_comb_tree((4, 4)) s = siteinds("S=1/2", g) χ = 2 - Random.seed!(1564) + seed!(1564) ψ = randomITensorNetwork(s; link_space=χ) ψψ = ψ ⊗ prime(dag(ψ); sites=[]) diff --git a/test/test_inner.jl b/test/test_inner.jl index 0b90a3cd..a9511a69 100644 --- a/test/test_inner.jl +++ b/test/test_inner.jl @@ -1,18 +1,18 @@ using Test -using Random -using Graphs -using NamedGraphs using ITensorNetworks -using SplitApplyCombine using ITensorNetworks: logscalar, scalar, inner, loginner using ITensors: siteinds, dag +using SplitApplyCombine: group +using Graphs: SimpleGraph, uniform_tree +using NamedGraphs: NamedGraph +using Random: seed! @testset "Inner products, BP vs exact comparison" begin - Random.seed!(1234) + seed!(1234) L = 4 χ = 2 - g = NamedGraph(Graphs.SimpleGraph(uniform_tree(L))) + g = NamedGraph(SimpleGraph(uniform_tree(L))) s = siteinds("S=1/2", g) y = randomITensorNetwork(s; link_space=χ) x = randomITensorNetwork(s; link_space=χ) From 691320e32399ce1adee5d34cfc4a637bd79cd848 Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 29 Mar 2024 16:08:18 -0400 Subject: [PATCH 16/35] Refactor test_apply and test_additensornetworks --- src/contract.jl | 10 ++++++---- src/inner.jl | 31 ++++++++++++++++--------------- test/test_additensornetworks.jl | 18 +++++++----------- test/test_apply.jl | 23 ++++++++--------------- test/test_belief_propagation.jl | 4 +++- 5 files changed, 40 insertions(+), 46 deletions(-) diff --git a/src/contract.jl b/src/contract.jl index 2b6ee6ae..11b0afad 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -1,5 +1,5 @@ using NamedGraphs: vertex_to_parent_vertex -using ITensors: ITensor +using ITensors: ITensor, scalar using ITensors.ContractionSequenceOptimization: deepmap using ITensors.NDTensors: NDTensors, Algorithm, @Algorithm_str, contract using LinearAlgebra: normalize! @@ -41,13 +41,15 @@ function contract_density_matrix( return out end -function scalar( +function ITensors.scalar( alg::Algorithm, tn::Union{AbstractITensorNetwork,Vector{ITensor}}; kwargs... ) return contract(alg, tn; kwargs...)[] end -function scalar(tn::Union{AbstractITensorNetwork,Vector{ITensor}}; alg="exact", kwargs...) +function ITensors.scalar( + tn::Union{AbstractITensorNetwork,Vector{ITensor}}; alg="exact", kwargs... +) return scalar(Algorithm(alg), tn; kwargs...) end @@ -90,6 +92,6 @@ function logscalar( end end -function scalar(alg::Algorithm"bp", tn::AbstractITensorNetwork; kwargs...) +function ITensors.scalar(alg::Algorithm"bp", tn::AbstractITensorNetwork; kwargs...) return exp(logscalar(alg, tn; kwargs...)) end diff --git a/src/inner.jl b/src/inner.jl index e165551a..5b7ba8b1 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -1,8 +1,9 @@ -default_algorithm(tns::Tuple) = "bp" +using ITensors: inner, scalar, loginner +using LinearAlgebra: norm, norm_sqr -#Default for map_linkinds should be sim. +default_algorithm(tns::Tuple) = "bp" -function inner( +function ITensors.inner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; alg=default_algorithm((ϕ, ψ)), @@ -11,7 +12,7 @@ function inner( return inner(Algorithm(alg), ϕ, ψ; kwargs...) end -function inner( +function ITensors.inner( ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; @@ -21,7 +22,7 @@ function inner( return inner(Algorithm(alg), ϕ, A, ψ; kwargs...) end -function inner( +function ITensors.inner( alg::Algorithm"exact", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; @@ -36,7 +37,7 @@ function inner( return scalar(tn; sequence) end -function inner( +function ITensors.inner( alg::Algorithm"exact", ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, @@ -52,7 +53,7 @@ function inner( return scalar(tn; sequence) end -function loginner( +function ITensors.loginner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; alg=default_algorithm(ϕ, ψ), @@ -61,7 +62,7 @@ function loginner( return loginner(Algorithm(alg), ϕ, ψ; kwargs...) end -function loginner( +function ITensors.loginner( ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; @@ -71,13 +72,13 @@ function loginner( return loginner(Algorithm(alg), ϕ, A, ψ; kwargs...) end -function loginner( +function ITensors.loginner( alg::Algorithm"exact", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; kwargs... ) return log(inner(alg, ϕ, ψ); kwargs...) end -function loginner( +function ITensors.loginner( alg::Algorithm"exact", ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, @@ -87,7 +88,7 @@ function loginner( return log(inner(alg, ϕ, A, ψ); kwargs...) end -function loginner( +function ITensors.loginner( alg::Algorithm"bp", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; @@ -98,7 +99,7 @@ function loginner( return logscalar(alg, tn; kwargs...) end -function loginner( +function ITensors.loginner( alg::Algorithm"bp", ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, @@ -110,7 +111,7 @@ function loginner( return logscalar(alg, tn; kwargs...) end -function inner( +function ITensors.inner( alg::Algorithm"bp", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; @@ -121,7 +122,7 @@ function inner( return scalar(alg, tn; kwargs...) end -function inner( +function ITensors.inner( alg::Algorithm"bp", ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, @@ -135,7 +136,7 @@ end # TODO: rename `sqnorm` to match https://github.com/JuliaStats/Distances.jl, # or `norm_sqr` to match `LinearAlgebra.norm_sqr` -norm_sqr(ψ::AbstractITensorNetwork; kwargs...) = inner(ψ, ψ; kwargs...) +LinearAlgebra.norm_sqr(ψ::AbstractITensorNetwork; kwargs...) = inner(ψ, ψ; kwargs...) function LinearAlgebra.norm(ψ::AbstractITensorNetwork; kwargs...) return sqrt(abs(real(norm_sqr(ψ; kwargs...)))) diff --git a/test/test_additensornetworks.jl b/test/test_additensornetworks.jl index bf92fbbf..c49ba0fc 100644 --- a/test/test_additensornetworks.jl +++ b/test/test_additensornetworks.jl @@ -1,18 +1,14 @@ -using ITensorNetworks using Test -using Compat -using ITensors -using Metis -using NamedGraphs -using NamedGraphs: hexagonal_lattice_graph, rem_edge! -using Random -using LinearAlgebra -using SplitApplyCombine +using ITensorNetworks -using Random +using NamedGraphs: hexagonal_lattice_graph, rem_edge!, named_grid, NamedEdge +using SplitApplyCombine: group +using Random: seed! +using ITensors: siteinds, scalar, apply, inner, op +using LinearAlgebra: norm_sqr @testset "add_itensornetworks" begin - Random.seed!(5623) + seed!(5623) g = named_grid((2, 2)) s = siteinds("S=1/2", g) ψ1 = ITensorNetwork(s, v -> "↑") diff --git a/test/test_apply.jl b/test/test_apply.jl index e247f727..177322a8 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,22 +1,15 @@ +using Test using ITensorNetworks using ITensorNetworks: - environment, - update, - inner, - BeliefPropagationCache, - VidalITensorNetwork, - norm_sqr_network_fast -using Test -using Compat -using ITensors -using Metis -using NamedGraphs -using Random -using LinearAlgebra -using SplitApplyCombine + environment, update, BeliefPropagationCache, VidalITensorNetwork, norm_sqr_network_fast + +using ITensors: inner +using SplitApplyCombine: group +using Random: seed! +using NamedGraphs: named_grid, disjoint_union @testset "apply" begin - Random.seed!(5623) + seed!(5623) g_dims = (2, 3) n = prod(g_dims) g = named_grid(g_dims) diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index 403ada65..2e66c908 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -25,7 +25,9 @@ using ITensors: ITensor, inds, combiner, - array + array, + dag, + inner disable_warn_order() From 755cbbd03d896f48f594c9722bcf475158832362 Mon Sep 17 00:00:00 2001 From: Joey Date: Sat, 30 Mar 2024 08:21:31 -0400 Subject: [PATCH 17/35] Formatting, namespace fix on test apply --- src/caches/beliefpropagationcache.jl | 11 +++++++++-- src/contract.jl | 2 +- test/test_apply.jl | 10 +++++----- 3 files changed, 15 insertions(+), 8 deletions(-) diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index 2f066c2e..fdbc8dc6 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -18,6 +18,9 @@ function default_partitioned_vertices(f::AbstractFormNetwork) return group(v -> original_state_vertex(f, v), vertices(f)) end default_cache_update_kwargs(cache) = (; maxiter=20, tol=1e-5) +function default_cache_construction_kwargs(alg::Algorithm"bp", ψ::AbstractITensorNetwork) + return (; partition_vertices=default_partitioned_vertices(ψ)) +end #TODO: Define a version of this that works for QN supporting tensors function message_diff(message_a::Vector{ITensor}, message_b::Vector{ITensor}) @@ -51,8 +54,12 @@ function BeliefPropagationCache( return BeliefPropagationCache(tn, partitioned_vertices; kwargs...) end -function cache(alg::Algorithm"bp", tn; kwargs...) - return BeliefPropagationCache(tn; kwargs...) +function cache( + alg::Algorithm"bp", + tn; + cache_construction_kwargs=default_cache_construction_kwargs(alg, tn), +) + return BeliefPropagationCache(tn; cache_construction_kwargs...) end function partitioned_itensornetwork(bp_cache::BeliefPropagationCache) diff --git a/src/contract.jl b/src/contract.jl index 11b0afad..d53f62c7 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -69,7 +69,7 @@ function logscalar( alg::Algorithm, tn::AbstractITensorNetwork; (cache!)=nothing, - cache_construction_kwargs=(;), + cache_construction_kwargs=default_cache_construction_kwargs(alg, tn), update_cache=isnothing(cache!), cache_update_kwargs=default_cache_update_kwargs(cache!), ) diff --git a/test/test_apply.jl b/test/test_apply.jl index 177322a8..20b1080c 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,12 +1,12 @@ using Test using ITensorNetworks using ITensorNetworks: - environment, update, BeliefPropagationCache, VidalITensorNetwork, norm_sqr_network_fast + environment, update, BeliefPropagationCache, VidalITensorNetwork, norm_sqr_network -using ITensors: inner +using ITensors: inner, siteinds, op, apply using SplitApplyCombine: group using Random: seed! -using NamedGraphs: named_grid, disjoint_union +using NamedGraphs: named_grid, PartitionVertex @testset "apply" begin seed!(5623) @@ -17,7 +17,7 @@ using NamedGraphs: named_grid, disjoint_union χ = 2 ψ = randomITensorNetwork(s; link_space=χ) v1, v2 = (2, 2), (1, 2) - ψψ = disjoint_union("bra" => dag(prime(ψ; sites=[])), "ket" => ψ) + ψψ = norm_sqr_network(ψ) #Simple Belief Propagation Grouping bp_cache = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) bp_cache = update(bp_cache; maxiter=20) @@ -35,7 +35,7 @@ using NamedGraphs: named_grid, disjoint_union ngates = 5 for i in 1:ngates - o = ITensors.op("RandomUnitary", s[v1]..., s[v2]...) + o = op("RandomUnitary", s[v1]..., s[v2]...) ψOexact = apply(o, ψ; cutoff=1e-16) ψOSBP = apply( From a7080eabc2913ce292fe4f72321ebcc6c08318be Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 1 Apr 2024 09:21:28 -0400 Subject: [PATCH 18/35] Better scalar definitions in contract.jl --- src/caches/beliefpropagationcache.jl | 2 +- src/contract.jl | 16 ++- src/formnetworks/bilinearformnetwork.jl | 2 - src/imports.jl | 109 ------------------ .../solvers/deprecated/projmpo_apply.jl | 1 - .../solvers/deprecated/projmpo_mps2.jl | 1 - .../solvers/deprecated/projmps2.jl | 1 - 7 files changed, 13 insertions(+), 119 deletions(-) delete mode 100644 src/imports.jl delete mode 100644 src/treetensornetworks/solvers/deprecated/projmpo_apply.jl delete mode 100644 src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl delete mode 100644 src/treetensornetworks/solvers/deprecated/projmps2.jl diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index fdbc8dc6..3180352b 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -19,7 +19,7 @@ function default_partitioned_vertices(f::AbstractFormNetwork) end default_cache_update_kwargs(cache) = (; maxiter=20, tol=1e-5) function default_cache_construction_kwargs(alg::Algorithm"bp", ψ::AbstractITensorNetwork) - return (; partition_vertices=default_partitioned_vertices(ψ)) + return (; partitioned_vertices=default_partitioned_vertices(ψ)) end #TODO: Define a version of this that works for QN supporting tensors diff --git a/src/contract.jl b/src/contract.jl index d53f62c7..27654c53 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -62,7 +62,15 @@ end function logscalar( alg::Algorithm"exact", tn::Union{AbstractITensorNetwork,Vector{ITensor}}; kwargs... ) - return log(complex(scalar(alg, tn; kwargs...))) + s = scalar(alg, tn; kwargs...) + if s ≈ 0 + tol = 1e-16 + return -Inf + elseif isa(s, AbstractFloat) && s >= 0 + return log(s) + else + return complex(s) + end end function logscalar( @@ -74,7 +82,7 @@ function logscalar( cache_update_kwargs=default_cache_update_kwargs(cache!), ) if isnothing(cache!) - cache! = Ref(cache(alg, tn; cache_construction_kwargs...)) + cache! = Ref(cache(alg, tn; cache_construction_kwargs)) end if update_cache @@ -83,9 +91,9 @@ function logscalar( numerator_terms, denominator_terms = scalar_factors(cache![]) terms = vcat(numerator_terms, denominator_terms) - if any(==(0), terms) + if any(≈(0), terms) return -Inf - elseif all(>=(0), terms) + elseif all(t -> isa(t, AbstractFloat), terms) && all(>=(0), terms) return sum(log.(numerator_terms)) - sum(log.((denominator_terms))) else return sum(log.(complex.(numerator_terms))) - sum(log.(complex.((denominator_terms)))) diff --git a/src/formnetworks/bilinearformnetwork.jl b/src/formnetworks/bilinearformnetwork.jl index aefb1511..b7408e4c 100644 --- a/src/formnetworks/bilinearformnetwork.jl +++ b/src/formnetworks/bilinearformnetwork.jl @@ -51,8 +51,6 @@ function Base.copy(blf::BilinearFormNetwork) ) end -#Is the ordering of the indices correct here? CHECK THIS -#Put bra into the vector space!!!! function BilinearFormNetwork( bra::AbstractITensorNetwork, ket::AbstractITensorNetwork; diff --git a/src/imports.jl b/src/imports.jl deleted file mode 100644 index 999201d8..00000000 --- a/src/imports.jl +++ /dev/null @@ -1,109 +0,0 @@ -import Base: - # types - Vector, - # functions - convert, - copy, - eltype, - getindex, - hvncat, - setindex!, - show, - isapprox, - isassigned, - iterate, - union, - + - -import NamedGraphs: - vertextype, - convert_vertextype, - vertex_to_parent_vertex, - rename_vertices, - disjoint_union, - mincut_partitions, - incident_edges, - boundary_partitionedges - -import .DataGraphs: - underlying_graph, - underlying_graph_type, - vertex_data, - edge_data, - edge_data_type, - reverse_data_direction - -import Graphs: SimpleGraph, is_directed, weights - -import KrylovKit: eigsolve, linsolve - -import LinearAlgebra: factorize, normalize, normalize!, qr, svd - -import Observers: update! - -import ITensors: - # contraction - apply, - contract, - dmrg, - orthogonalize, - isortho, - inner, - loginner, - norm, - lognorm, - expect, - # truncation - truncate, - replacebond!, - replacebond, - # site and link indices - siteind, - siteinds, - linkinds, - # index set functions - uniqueinds, - commoninds, - replaceinds, - hascommoninds, - # priming and tagging - adjoint, - sim, - prime, - setprime, - noprime, - replaceprime, - addtags, - removetags, - replacetags, - settags, - tags, - # dag - dag, - # permute - permute, - #commoninds - hascommoninds, - # linkdims - linkdim, - linkdims, - maxlinkdim, - # projected operators - product, - nsite, - # promotion and conversion - promote_itensor_eltype, - scalar, - scalartype, - #adding - add - -import ITensors.LazyApply: - # extracting terms from a sum - terms -#Algorithm -Algorithm - -using ITensors.ContractionSequenceOptimization: deepmap - -import ITensors.ITensorVisualizationCore: visualize diff --git a/src/treetensornetworks/solvers/deprecated/projmpo_apply.jl b/src/treetensornetworks/solvers/deprecated/projmpo_apply.jl deleted file mode 100644 index 8b137891..00000000 --- a/src/treetensornetworks/solvers/deprecated/projmpo_apply.jl +++ /dev/null @@ -1 +0,0 @@ - diff --git a/src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl b/src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl deleted file mode 100644 index 8b137891..00000000 --- a/src/treetensornetworks/solvers/deprecated/projmpo_mps2.jl +++ /dev/null @@ -1 +0,0 @@ - diff --git a/src/treetensornetworks/solvers/deprecated/projmps2.jl b/src/treetensornetworks/solvers/deprecated/projmps2.jl deleted file mode 100644 index 8b137891..00000000 --- a/src/treetensornetworks/solvers/deprecated/projmps2.jl +++ /dev/null @@ -1 +0,0 @@ - From d94900cfc43fb3c79b9625ae058a58c9bb9e4990 Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 1 Apr 2024 15:15:16 -0400 Subject: [PATCH 19/35] Alphabetize namespaces, polish up --- src/abstractitensornetwork.jl | 43 +------------------ .../binary_tree_partition.jl | 2 +- src/boundarymps.jl | 2 +- src/caches/beliefpropagationcache.jl | 34 ++++++--------- src/contract.jl | 28 ++++-------- src/inner.jl | 10 ++--- test/test_additensornetworks.jl | 10 ++--- test/test_apply.jl | 12 +++--- test/test_belief_propagation.jl | 36 ++++++++-------- test/test_contraction_sequence.jl | 12 +++--- test/test_inner.jl | 8 ++-- 11 files changed, 67 insertions(+), 130 deletions(-) diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index c3cdb3f8..fad6436b 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -724,23 +724,6 @@ function norm_sqr_network(ψ::AbstractITensorNetwork) return disjoint_union("bra" => dag(prime(ψ; sites=[])), "ket" => ψ) end -#Ideally this will not be necessary but this is a temporary fast version to avoid the overhead of `disjoint_union` -function norm_sqr_network_fast(ψ::AbstractITensorNetwork) - ψbra = rename_vertices(v -> (v, 1), data_graph(ψ)) - ψdag = copy(ψ) - for v in vertices(ψdag) - setindex_preserve_graph!(ψdag, dag(ψdag[v]), v) - end - ψket = rename_vertices(v -> (v, 2), data_graph(prime(ψdag; sites=[]))) - ψψ = ITensorNetwork(union(ψbra, ψket)) - for v in vertices(ψ) - if !isempty(commoninds(ψψ[(v, 1)], ψψ[(v, 2)])) - add_edge!(ψψ, (v, 1) => (v, 2)) - end - end - return ψψ -end - # # Printing # @@ -911,28 +894,4 @@ end Base.:+(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork) = add(tn1, tn2) -ITensors.hasqns(tn::AbstractITensorNetwork) = all([hasqns(tn[v]) for v in vertices(tn)]) - -## # TODO: should this make sure that internal indices -## # don't clash? -## function hvncat( -## dim::Int, tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork; new_dim_names=(1, 2) -## ) -## dg = hvncat(dim, data_graph(tn1), data_graph(tn2); new_dim_names) -## -## # Add in missing edges that may be shared -## # across `tn1` and `tn2`. -## vertices1 = vertices(dg)[1:nv(tn1)] -## vertices2 = vertices(dg)[(nv(tn1) + 1):end] -## for v1 in vertices1, v2 in vertices2 -## if hascommoninds(dg[v1], dg[v2]) -## add_edge!(dg, v1 => v2) -## end -## end -## -## # TODO: Allow customization of the output type. -## ## return promote_type(typeof(tn1), typeof(tn2))(dg) -## ## return contract_output(typeof(tn1), typeof(tn2))(dg) -## -## return ITensorNetwork(dg) -## end +ITensors.hasqns(tn::AbstractITensorNetwork) = any(v -> hasqns(tn[v]), vertices(tn)) diff --git a/src/approx_itensornetwork/binary_tree_partition.jl b/src/approx_itensornetwork/binary_tree_partition.jl index c4ee00a8..7eccb4f6 100644 --- a/src/approx_itensornetwork/binary_tree_partition.jl +++ b/src/approx_itensornetwork/binary_tree_partition.jl @@ -129,6 +129,6 @@ function _partition( return rename_vertices(par, name_map) end -function _partition(tn::ITensorNetwork, inds_btree::DataGraph; alg::String) +function _partition(tn::ITensorNetwork, inds_btree::DataGraph; alg) return _partition(Algorithm(alg), tn, inds_btree) end diff --git a/src/boundarymps.jl b/src/boundarymps.jl index 7afd5ac7..8e27868e 100644 --- a/src/boundarymps.jl +++ b/src/boundarymps.jl @@ -9,5 +9,5 @@ function contract_boundary_mps(tn::ITensorNetwork; kwargs...) end T = MPO([tn[i1, d2 - 1] for i1 in 1:d1]) vR = MPS([tn[i1, d2] for i1 in 1:d1]) - return ITensors.inner(dag(vL), T, vR)[] + return inner(dag(vL), T, vR)[] end diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index 3180352b..e5cf67ee 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -54,12 +54,8 @@ function BeliefPropagationCache( return BeliefPropagationCache(tn, partitioned_vertices; kwargs...) end -function cache( - alg::Algorithm"bp", - tn; - cache_construction_kwargs=default_cache_construction_kwargs(alg, tn), -) - return BeliefPropagationCache(tn; cache_construction_kwargs...) +function cache(alg::Algorithm"bp", tn; kwargs...) + return BeliefPropagationCache(tn; kwargs...) end function partitioned_itensornetwork(bp_cache::BeliefPropagationCache) @@ -263,32 +259,26 @@ function update_factor(bp_cache, vertex, factor) return update_factors(bp_cache, [vertex], ITensor[factor]) end -function vertex_norm(bp_cache::BeliefPropagationCache, pv::PartitionVertex) +function region_scalar(bp_cache::BeliefPropagationCache, pv::PartitionVertex) incoming_mts = environment(bp_cache, [pv]) local_state = factor(bp_cache, pv) return scalar(vcat(incoming_mts, local_state)) end -function edge_norm(bp_cache::BeliefPropagationCache, pe::PartitionEdge) +function region_scalar(bp_cache::BeliefPropagationCache, pe::PartitionEdge) return scalar(vcat(message(bp_cache, pe), message(bp_cache, reverse(pe)))) end -function vertex_norms(bp_cache::BeliefPropagationCache, pvs::Vector) - return [vertex_norm(bp_cache, pv) for pv in pvs] -end - -function vertex_norms(bp_cache::BeliefPropagationCache) - return vertex_norms(bp_cache, partitionvertices(partitioned_itensornetwork(bp_cache))) -end - -function edge_norms(bp_cache::BeliefPropagationCache, pes::Vector) - return [edge_norm(bp_cache, pe) for pe in pes] +function vertex_scalars(bp_cache::BeliefPropagationCache, pvs::Vector) + return [region_scalar(bp_cache, pv) for pv in pvs] end -function edge_norms(bp_cache::BeliefPropagationCache) - return edge_norms(bp_cache, partitionedges(partitioned_itensornetwork(bp_cache))) +function edge_scalars(bp_cache::BeliefPropagationCache, pes::Vector) + return [region_scalar(bp_cache, pe) for pe in pes] end -function scalar_factors(bp_cache::BeliefPropagationCache) - return vertex_norms(bp_cache), edge_norms(bp_cache) +function scalar_factors_quotient(bp_cache::BeliefPropagationCache) + pvs = partitionvertices(partitioned_itensornetwork(bp_cache)) + pes = partitionedges(partitioned_itensornetwork(bp_cache)) + return vertex_scalars(bp_cache, pvs), edge_scalars(bp_cache, pes) end diff --git a/src/contract.jl b/src/contract.jl index 27654c53..d3aa25c4 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -4,7 +4,7 @@ using ITensors.ContractionSequenceOptimization: deepmap using ITensors.NDTensors: NDTensors, Algorithm, @Algorithm_str, contract using LinearAlgebra: normalize! -function NDTensors.contract(tn::AbstractITensorNetwork; alg::String="exact", kwargs...) +function NDTensors.contract(tn::AbstractITensorNetwork; alg="exact", kwargs...) return contract(Algorithm(alg), tn; kwargs...) end @@ -54,7 +54,7 @@ function ITensors.scalar( end function logscalar( - tn::Union{AbstractITensorNetwork,Vector{ITensor}}; alg::String="exact", kwargs... + tn::Union{AbstractITensorNetwork,Vector{ITensor}}; alg="exact", kwargs... ) return logscalar(Algorithm(alg), tn; kwargs...) end @@ -63,14 +63,8 @@ function logscalar( alg::Algorithm"exact", tn::Union{AbstractITensorNetwork,Vector{ITensor}}; kwargs... ) s = scalar(alg, tn; kwargs...) - if s ≈ 0 - tol = 1e-16 - return -Inf - elseif isa(s, AbstractFloat) && s >= 0 - return log(s) - else - return complex(s) - end + s = real(s) < 0 ? complex(s) : s + return log(s) end function logscalar( @@ -82,22 +76,18 @@ function logscalar( cache_update_kwargs=default_cache_update_kwargs(cache!), ) if isnothing(cache!) - cache! = Ref(cache(alg, tn; cache_construction_kwargs)) + cache! = Ref(cache(alg, tn; cache_construction_kwargs...)) end if update_cache cache![] = update(cache![]; cache_update_kwargs...) end - numerator_terms, denominator_terms = scalar_factors(cache![]) + numerator_terms, denominator_terms = scalar_factors_quotient(cache![]) terms = vcat(numerator_terms, denominator_terms) - if any(≈(0), terms) - return -Inf - elseif all(t -> isa(t, AbstractFloat), terms) && all(>=(0), terms) - return sum(log.(numerator_terms)) - sum(log.((denominator_terms))) - else - return sum(log.(complex.(numerator_terms))) - sum(log.(complex.((denominator_terms)))) - end + terms .= any(t -> real(t) < 0, terms) ? complex.(terms) : terms + + return sum(log.(numerator_terms)) - sum(log.((denominator_terms))) end function ITensors.scalar(alg::Algorithm"bp", tn::AbstractITensorNetwork; kwargs...) diff --git a/src/inner.jl b/src/inner.jl index 5b7ba8b1..166a2c6f 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -1,12 +1,12 @@ using ITensors: inner, scalar, loginner using LinearAlgebra: norm, norm_sqr -default_algorithm(tns::Tuple) = "bp" +default_contract_alg(tns::Tuple) = "bp" function ITensors.inner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - alg=default_algorithm((ϕ, ψ)), + alg=default_contract_alg((ϕ, ψ)), kwargs..., ) return inner(Algorithm(alg), ϕ, ψ; kwargs...) @@ -16,7 +16,7 @@ function ITensors.inner( ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - alg=default_algorithm((ϕ, A, ψ)), + alg=default_contract_alg((ϕ, A, ψ)), kwargs..., ) return inner(Algorithm(alg), ϕ, A, ψ; kwargs...) @@ -56,7 +56,7 @@ end function ITensors.loginner( ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - alg=default_algorithm(ϕ, ψ), + alg=default_contract_alg((ϕ, ψ)), kwargs..., ) return loginner(Algorithm(alg), ϕ, ψ; kwargs...) @@ -66,7 +66,7 @@ function ITensors.loginner( ϕ::AbstractITensorNetwork, A::AbstractITensorNetwork, ψ::AbstractITensorNetwork; - alg=default_algorithm(ϕ, ψ), + alg=default_contract_alg((ϕ, A, ψ)), kwargs..., ) return loginner(Algorithm(alg), ϕ, A, ψ; kwargs...) diff --git a/test/test_additensornetworks.jl b/test/test_additensornetworks.jl index c49ba0fc..9b52d88b 100644 --- a/test/test_additensornetworks.jl +++ b/test/test_additensornetworks.jl @@ -1,14 +1,14 @@ -using Test using ITensorNetworks -using NamedGraphs: hexagonal_lattice_graph, rem_edge!, named_grid, NamedEdge +using Test: @test, @testset +using NamedGraphs: NamedEdge, hexagonal_lattice_graph, named_grid, rem_edge! using SplitApplyCombine: group -using Random: seed! -using ITensors: siteinds, scalar, apply, inner, op +using Random: Random +using ITensors: apply, inner, op, scalar, siteinds using LinearAlgebra: norm_sqr @testset "add_itensornetworks" begin - seed!(5623) + Random.seed!(5623) g = named_grid((2, 2)) s = siteinds("S=1/2", g) ψ1 = ITensorNetwork(s, v -> "↑") diff --git a/test/test_apply.jl b/test/test_apply.jl index 20b1080c..4eb728b9 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,15 +1,15 @@ -using Test using ITensorNetworks using ITensorNetworks: - environment, update, BeliefPropagationCache, VidalITensorNetwork, norm_sqr_network + BeliefPropagationCache, VidalITensorNetwork, environment, norm_sqr_network, update -using ITensors: inner, siteinds, op, apply +using Test: @test, @testset +using ITensors: apply, op, inner, siteinds using SplitApplyCombine: group -using Random: seed! -using NamedGraphs: named_grid, PartitionVertex +using Random: Random +using NamedGraphs: PartitionVertex, named_grid @testset "apply" begin - seed!(5623) + Random.seed!(5623) g_dims = (2, 3) n = prod(g_dims) g = named_grid(g_dims) diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index 2e66c908..d6608e18 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -1,33 +1,33 @@ -using Test using ITensorNetworks using ITensorNetworks: + BeliefPropagationCache, + contract, + contract_boundary_mps, + environment, ising_network, split_index, - contract_boundary_mps, - BeliefPropagationCache, tensornetwork, update, update_factor, - environment, - contract +using Test: @test, @testset using LinearAlgebra: eigvals, tr -using NamedGraphs: named_grid, named_comb_tree, PartitionVertex, NamedEdge -using Random: seed! +using NamedGraphs: NamedEdge, PartitionVertex, named_comb_tree, named_grid +using Random: Random using SplitApplyCombine: group using ITensors: - siteinds, - apply, - randomITensor, - disable_warn_order, - prime, - op, ITensor, - inds, - combiner, + apply, array, + combiner, dag, - inner + disable_warn_order, + inds, + inner, + op, + prime, + randomITensor, + siteinds, disable_warn_order() @@ -38,7 +38,7 @@ disable_warn_order() g = named_grid(g_dims) s = siteinds("S=1/2", g) χ = 4 - seed!(1234) + Random.seed!(1234) ψ = randomITensorNetwork(s; link_space=χ) ψψ = ψ ⊗ prime(dag(ψ); sites=[]) @@ -68,7 +68,7 @@ disable_warn_order() g = named_comb_tree((4, 4)) s = siteinds("S=1/2", g) χ = 2 - seed!(1564) + Random.seed!(1564) ψ = randomITensorNetwork(s; link_space=χ) ψψ = ψ ⊗ prime(dag(ψ); sites=[]) diff --git a/test/test_contraction_sequence.jl b/test/test_contraction_sequence.jl index aae6bd0e..7474bb75 100644 --- a/test/test_contraction_sequence.jl +++ b/test/test_contraction_sequence.jl @@ -1,17 +1,15 @@ -using ITensors using ITensorNetworks -using OMEinsumContractionOrders -using Random -using Test -using EinExprs: Exhaustive, Greedy, HyPar - using ITensorNetworks: norm_sqr_network -Random.seed!(1234) +using ITensors: contract, siteinds +using EinExprs: Exhaustive, Greedy, HyPar +using Random: Random +using Test: @test, @testset ITensors.disable_warn_order() @testset "contraction_sequence" begin + Random.seed!(1234) dims = (2, 3) g = named_grid(dims) s = siteinds("S=1/2", g) diff --git a/test/test_inner.jl b/test/test_inner.jl index a9511a69..f1a68c78 100644 --- a/test/test_inner.jl +++ b/test/test_inner.jl @@ -1,15 +1,15 @@ using Test using ITensorNetworks -using ITensorNetworks: logscalar, scalar, inner, loginner -using ITensors: siteinds, dag +using ITensorNetworks: inner, loginner, logscalar, scalar +using ITensors: dag, siteinds using SplitApplyCombine: group using Graphs: SimpleGraph, uniform_tree using NamedGraphs: NamedGraph -using Random: seed! +using Random: Random @testset "Inner products, BP vs exact comparison" begin - seed!(1234) + Random.seed!(1234) L = 4 χ = 2 g = NamedGraph(SimpleGraph(uniform_tree(L))) From 9eb55c4bb5f5cee404b6208308a1ea04a275329d Mon Sep 17 00:00:00 2001 From: Joey Date: Mon, 1 Apr 2024 15:22:00 -0400 Subject: [PATCH 20/35] Formatting --- test/test_belief_propagation.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index d6608e18..3cf961fb 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -8,7 +8,7 @@ using ITensorNetworks: split_index, tensornetwork, update, - update_factor, + update_factor using Test: @test, @testset using LinearAlgebra: eigvals, tr @@ -17,19 +17,18 @@ using Random: Random using SplitApplyCombine: group using ITensors: ITensor, - apply, + apply, array, combiner, dag, disable_warn_order, - inds, + inds, inner, op, prime, randomITensor, siteinds, - -disable_warn_order() + disable_warn_order() @testset "belief_propagation" begin From 9b33331f00813f349ed69e4d3fc23de15d23e037 Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 2 Apr 2024 08:50:08 -0400 Subject: [PATCH 21/35] Bug Fix --- test/test_belief_propagation.jl | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index 3cf961fb..8f625246 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -16,21 +16,10 @@ using NamedGraphs: NamedEdge, PartitionVertex, named_comb_tree, named_grid using Random: Random using SplitApplyCombine: group using ITensors: - ITensor, - apply, - array, - combiner, - dag, - disable_warn_order, - inds, - inner, - op, - prime, - randomITensor, - siteinds, - disable_warn_order() + ITensor, apply, array, combiner, dag, inds, inner, op, prime, randomITensor, siteinds @testset "belief_propagation" begin + ITensors.disable_warn_order() #First test on an MPS, should be exact g_dims = (1, 6) From be39bc4bdd172c4053df90fc59594c091d0b593c Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 2 Apr 2024 10:07:27 -0400 Subject: [PATCH 22/35] Bug Fix --- test/test_belief_propagation.jl | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index 8f625246..ba0c3d46 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -16,7 +16,18 @@ using NamedGraphs: NamedEdge, PartitionVertex, named_comb_tree, named_grid using Random: Random using SplitApplyCombine: group using ITensors: - ITensor, apply, array, combiner, dag, inds, inner, op, prime, randomITensor, siteinds + ITensors, + ITensor, + apply, + array, + combiner, + dag, + inds, + inner, + op, + prime, + randomITensor, + siteinds @testset "belief_propagation" begin ITensors.disable_warn_order() From 84da9d709c3525960982c1164fd1196bf88519c2 Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 2 Apr 2024 14:42:47 -0400 Subject: [PATCH 23/35] Fixed logscalar bug --- src/caches/beliefpropagationcache.jl | 2 +- src/contract.jl | 9 +++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index e5cf67ee..909169de 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -1,5 +1,5 @@ using ITensors.ITensorMPS: ITensorMPS -using NamedGraphs: boundary_partitionedges +using NamedGraphs: boundary_partitionedges, partitionvertices, partitionedges default_message(inds_e) = ITensor[denseblocks(delta(inds_e))] default_messages(ptn::PartitionedGraph) = Dictionary() diff --git a/src/contract.jl b/src/contract.jl index d3aa25c4..8919bba8 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -84,8 +84,13 @@ function logscalar( end numerator_terms, denominator_terms = scalar_factors_quotient(cache![]) - terms = vcat(numerator_terms, denominator_terms) - terms .= any(t -> real(t) < 0, terms) ? complex.(terms) : terms + numerator_terms = + any(t -> real(t) < 0, numerator_terms) ? complex.(numerator_terms) : numerator_terms + denominator_terms = if any(t -> real(t) < 0, denominator_terms) + complex.(denominator_terms) + else + denominator_terms + end return sum(log.(numerator_terms)) - sum(log.((denominator_terms))) end From 4bf3bab53cf4f035f140adbf6c5f35514f98ad3c Mon Sep 17 00:00:00 2001 From: Joey Date: Tue, 2 Apr 2024 16:17:50 -0400 Subject: [PATCH 24/35] Quick Bug Fix, add_itensornetworks --- src/abstractitensornetwork.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index fad6436b..a5cdb2d4 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -869,7 +869,7 @@ function ITensorMPS.add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork #Create vertices of tn12 as direct sum of tn1[v] and tn2[v]. Work out the matching indices by matching edges. Make index tags those of tn1[v] for v in vertices(tn1) - @assert siteinds(tn1, v) == siteinds(tn2, v) + @assert issetequal(siteinds(tn1, v), siteinds(tn2, v)) e1_v = filter(x -> src(x) == v || dst(x) == v, edges_tn1) e2_v = filter(x -> src(x) == v || dst(x) == v, edges_tn2) From c750a080978ab99a4ec65f840128bc8fd4de6c54 Mon Sep 17 00:00:00 2001 From: Joey Date: Wed, 3 Apr 2024 10:51:43 -0400 Subject: [PATCH 25/35] Merged upstream changes --- examples/dynamics/2d_ising_imag_tebd.jl | 84 ------------------------- 1 file changed, 84 deletions(-) delete mode 100644 examples/dynamics/2d_ising_imag_tebd.jl diff --git a/examples/dynamics/2d_ising_imag_tebd.jl b/examples/dynamics/2d_ising_imag_tebd.jl deleted file mode 100644 index 57011b2b..00000000 --- a/examples/dynamics/2d_ising_imag_tebd.jl +++ /dev/null @@ -1,84 +0,0 @@ -using ITensors -using ITensorNetworks -using ITensorUnicodePlots -using UnicodePlots -using Random - -Random.seed!(1234) - -ITensors.disable_warn_order() - -system_dims = (2, 3) -n = prod(system_dims) -g = named_grid(system_dims) - -h = 2.0 - -@show h -@show system_dims - -s = siteinds("S=1/2", g) - -# -# DMRG comparison -# - -g_dmrg = rename_vertices(g, cartesian_to_linear(system_dims)) -ℋ_dmrg = ising(g_dmrg; h) -s_dmrg = [only(s[v]) for v in vertices(s)] -H_dmrg = MPO(ℋ_dmrg, s_dmrg) -ψ_dmrg_init = MPS(s_dmrg, j -> "↑") -@show inner(ψ_dmrg_init', H_dmrg, ψ_dmrg_init) -E_dmrg, ψ_dmrg = dmrg( - H_dmrg, ψ_dmrg_init; nsweeps=20, maxdim=[fill(10, 10); 20], cutoff=1e-8 -) -@show E_dmrg -Z_dmrg = reshape(expect(ψ_dmrg, "Z"), system_dims) - -display(Z_dmrg) -display(heatmap(Z_dmrg)) - -# -# PEPS TEBD optimization -# - -ℋ = ising(g; h) - -χ = 2 - -# Enable orthogonalizing the PEPS using a local gauge transformation -ortho = true - -ψ_init = ITensorNetwork(s, v -> "↑") - -β = 1.0 -Δβ = 0.1 - -println("maxdim = $χ") -@show β, Δβ -@show ortho - -println("\nFirst run TEBD without orthogonalization") -ψ = @time tebd( - group_terms(ℋ, g), ψ_init; β, Δβ, cutoff=1e-8, maxdim=χ, ortho=false, print_frequency=1 -) - -println("\nMeasure energy expectation value") -E = @time expect(ℋ, ψ) -@show E - -println("\nThen run TEBD with orthogonalization (more accurate)") -ψ = @time tebd( - group_terms(ℋ, g), ψ_init; β, Δβ, cutoff=1e-8, maxdim=χ, ortho, print_frequency=1 -) - -println("\nMeasure energy expectation value") -E = @time expect(ℋ, ψ) -@show E - -println("\nMeasure magnetization") -Z_dict = @time expect("Z", ψ) -Z = [Z_dict[Tuple(I)] for I in CartesianIndices(system_dims)] -display(Z) -display(heatmap(Z)) - From 4227726a35219b666ac1bb51cf81d14b58c96864 Mon Sep 17 00:00:00 2001 From: Joey Date: Wed, 3 Apr 2024 12:55:44 -0400 Subject: [PATCH 26/35] Fix tests to be compatable with upstream changes --- test/test_additensornetworks.jl | 3 ++- test/test_apply.jl | 5 ++--- test/test_belief_propagation.jl | 3 +-- test/test_gauging.jl | 3 +-- 4 files changed, 6 insertions(+), 8 deletions(-) diff --git a/test/test_additensornetworks.jl b/test/test_additensornetworks.jl index c123ebcd..9689a38c 100644 --- a/test/test_additensornetworks.jl +++ b/test/test_additensornetworks.jl @@ -2,7 +2,8 @@ using Graphs: rem_edge!, vertices using NamedGraphs: NamedEdge, hexagonal_lattice_graph, named_grid using ITensorNetworks: ITensorNetwork, inner_network, randomITensorNetwork, siteinds -using ITensors: ITensors, apply, op +using ITensors: ITensors, apply, op, scalar, inner +using LinearAlgebra: norm_sqr using Random: Random using Test: @test, @testset diff --git a/test/test_apply.jl b/test/test_apply.jl index 698d710b..7251a604 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -6,13 +6,12 @@ using ITensorNetworks: ITensorNetwork, VidalITensorNetwork, apply, - contract_inner, environment, - norm_network, + norm_sqr_network, randomITensorNetwork, siteinds, update -using ITensors: ITensors +using ITensors: ITensors, inner, op using NamedGraphs: PartitionVertex, named_grid using Random: Random using SplitApplyCombine: group diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index 56f4f6e4..a971b9f4 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -9,7 +9,6 @@ using ITensorNetworks: apply, combine_linkinds, contract, - contract_inner, contract_boundary_mps, contraction_sequence, environment, @@ -22,7 +21,7 @@ using ITensorNetworks: tensornetwork, update, update_factor -using ITensors: ITensors, ITensor, combiner, dag, inds, op, prime, randomITensor +using ITensors: ITensors, ITensor, combiner, dag, inner, inds, op, prime, randomITensor using ITensors.NDTensors: array using LinearAlgebra: eigvals, tr using NamedGraphs: NamedEdge, PartitionVertex, named_comb_tree, named_grid diff --git a/test/test_gauging.jl b/test/test_gauging.jl index a9b437ae..33b4cd8b 100644 --- a/test/test_gauging.jl +++ b/test/test_gauging.jl @@ -4,13 +4,12 @@ using ITensorNetworks: BeliefPropagationCache, ITensorNetwork, VidalITensorNetwork, - contract_inner, gauge_error, messages, randomITensorNetwork, siteinds, update -using ITensors: diagITensor, inds +using ITensors: diagITensor, inds, inner using ITensors.NDTensors: vector using LinearAlgebra: diag using NamedGraphs: named_grid From 0664f247e82f07de26eb4ad17f5c5ee09e1c0630 Mon Sep 17 00:00:00 2001 From: Joey Date: Wed, 3 Apr 2024 13:40:11 -0400 Subject: [PATCH 27/35] Fix tests to be compatable with upstream changes --- test/test_belief_propagation.jl | 2 +- test/test_inner.jl | 3 ++- test/test_itensornetwork.jl | 2 ++ test/test_tno.jl | 3 +-- 4 files changed, 6 insertions(+), 4 deletions(-) diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index a971b9f4..3fc22cd3 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -21,7 +21,7 @@ using ITensorNetworks: tensornetwork, update, update_factor -using ITensors: ITensors, ITensor, combiner, dag, inner, inds, op, prime, randomITensor +using ITensors: ITensors, ITensor, combiner, dag, inds, inner, op, prime, randomITensor using ITensors.NDTensors: array using LinearAlgebra: eigvals, tr using NamedGraphs: NamedEdge, PartitionVertex, named_comb_tree, named_grid diff --git a/test/test_inner.jl b/test/test_inner.jl index f1a68c78..976f5d3d 100644 --- a/test/test_inner.jl +++ b/test/test_inner.jl @@ -1,7 +1,8 @@ using Test using ITensorNetworks -using ITensorNetworks: inner, loginner, logscalar, scalar +using ITensorNetworks: + TTN, inner, inner_network, loginner, logscalar, randomITensorNetwork, scalar using ITensors: dag, siteinds using SplitApplyCombine: group using Graphs: SimpleGraph, uniform_tree diff --git a/test/test_itensornetwork.jl b/test/test_itensornetwork.jl index e8f2ebc3..da771763 100644 --- a/test/test_itensornetwork.jl +++ b/test/test_itensornetwork.jl @@ -25,6 +25,7 @@ using ITensors: hascommoninds, hasinds, inds, + inner, itensor, order, sim, @@ -40,6 +41,7 @@ using ITensorNetworks: inner_network, internalinds, linkinds, + norm_sqr_network, orthogonalize, randomITensorNetwork, siteinds diff --git a/test/test_tno.jl b/test/test_tno.jl index 727ba80c..f535c343 100644 --- a/test/test_tno.jl +++ b/test/test_tno.jl @@ -2,7 +2,6 @@ using Graphs: vertices using ITensorNetworks: apply, - contract_inner, flatten_networks, group_commuting_itensors, gate_group_to_tno, @@ -10,7 +9,7 @@ using ITensorNetworks: ising, randomITensorNetwork, siteinds -using ITensors: ITensor, noprime +using ITensors: ITensor, inner, noprime using NamedGraphs: named_grid using Test: @test, @testset From 388c15d3913f6a48e1e12a92bbb69a9d730b957a Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 5 Apr 2024 10:33:38 -0400 Subject: [PATCH 28/35] Refactor due to upstream changes --- README.md | 12 ++++----- examples/contract_bp.jl | 40 ++++++++++++++++++++++++++++ src/caches/beliefpropagationcache.jl | 4 +-- src/utility.jl | 1 + test/test_additensornetworks.jl | 2 +- test/test_forms.jl | 4 +-- test/test_inner.jl | 16 ++++++++--- 7 files changed, 64 insertions(+), 15 deletions(-) create mode 100644 examples/contract_bp.jl diff --git a/README.md b/README.md index 8a92276d..fe3fd65c 100644 --- a/README.md +++ b/README.md @@ -39,7 +39,7 @@ julia> using ITensorNetworks: ITensorNetwork, siteinds julia> using NamedGraphs: named_grid, subgraph julia> tn = ITensorNetwork(named_grid(4); link_space=2) -ITensorNetworks.ITensorNetwork{Int64} with 4 vertices: +ITensorNetwork{Int64} with 4 vertices: 4-element Vector{Int64}: 1 2 @@ -90,7 +90,7 @@ and here is a similar example for making a tensor network on a grid (a tensor pr ```julia julia> tn = ITensorNetwork(named_grid((2, 2)); link_space=2) -ITensorNetworks.ITensorNetwork{Tuple{Int64, Int64}} with 4 vertices: +ITensorNetwork{Tuple{Int64, Int64}} with 4 vertices: 4-element Vector{Tuple{Int64, Int64}}: (1, 1) (2, 1) @@ -125,7 +125,7 @@ julia> neighbors(tn, (1, 2)) (2, 2) julia> tn_1 = subgraph(v -> v[1] == 1, tn) -ITensorNetworks.ITensorNetwork{Tuple{Int64, Int64}} with 2 vertices: +ITensorNetwork{Tuple{Int64, Int64}} with 2 vertices: 2-element Vector{Tuple{Int64, Int64}}: (1, 1) (1, 2) @@ -139,7 +139,7 @@ with vertex data: (1, 2) │ ((dim=2|id=598|"1×1,1×2"), (dim=2|id=683|"1×2,2×2")) julia> tn_2 = subgraph(v -> v[1] == 2, tn) -ITensorNetworks.ITensorNetwork{Tuple{Int64, Int64}} with 2 vertices: +ITensorNetwork{Tuple{Int64, Int64}} with 2 vertices: 2-element Vector{Tuple{Int64, Int64}}: (2, 1) (2, 2) @@ -184,7 +184,7 @@ and edge data: 0-element Dictionaries.Dictionary{NamedGraphs.NamedEdge{Int64}, Vector{ITensors.Index}} julia> tn1 = ITensorNetwork(s; link_space=2) -ITensorNetworks.ITensorNetwork{Int64} with 3 vertices: +ITensorNetwork{Int64} with 3 vertices: 3-element Vector{Int64}: 1 2 @@ -201,7 +201,7 @@ with vertex data: 3 │ ((dim=2|id=558|"S=1/2,Site,n=3"), (dim=2|id=430|"2,3")) julia> tn2 = ITensorNetwork(s; link_space=2) -ITensorNetworks.ITensorNetwork{Int64} with 3 vertices: +ITensorNetwork{Int64} with 3 vertices: 3-element Vector{Int64}: 1 2 diff --git a/examples/contract_bp.jl b/examples/contract_bp.jl new file mode 100644 index 00000000..17eea650 --- /dev/null +++ b/examples/contract_bp.jl @@ -0,0 +1,40 @@ +# https://docs.julialang.org/en/v1/base/base/#Base.names +function imported_names(M::Module) + return setdiff(names(M; imported=true), names(M; imported=false)) +end + +function imported_functions(M::Module) + return map(f -> getfield(M, f), imported_names(M)) +end + +using FileUtils: replace_in_files + +function prepend_imports(M::Module) + fs = imported_functions(M) + # foreach(f -> @show((parentmodule(f), f)), fs) + src_path = joinpath(pkgdir(M), "src") + @show src_path + for f in fs + function_name = last(split(string(f), ".")) + module_name = last(split(string(parentmodule(f)), ".")) + annotated_function_name = module_name * "." * function_name + @show function_name, annotated_function_name + replace_in_files( + src_path, + "function $(function_name)(" => "function $(annotated_function_name)("; + recursive=true, + ignore_dirs=[".git"], + showdiffs=false, + ) + replace_in_files( + src_path, + "\n$(function_name)(" => "\n$(annotated_function_name)("; + recursive=true, + ignore_dirs=[".git"], + showdiffs=false, + ) + end +end + +using ITensorNetworks: ITensorNetworks +prepend_imports(ITensorNetworks) diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index 60014959..fbab6c67 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -285,7 +285,7 @@ function edge_scalars(bp_cache::BeliefPropagationCache, pes::Vector) end function scalar_factors_quotient(bp_cache::BeliefPropagationCache) - pvs = partitionvertices(partitioned_itensornetwork(bp_cache)) - pes = partitionedges(partitioned_itensornetwork(bp_cache)) + pvs = partitionvertices(partitioned_tensornetwork(bp_cache)) + pes = partitionedges(partitioned_tensornetwork(bp_cache)) return vertex_scalars(bp_cache, pvs), edge_scalars(bp_cache, pes) end diff --git a/src/utility.jl b/src/utility.jl index 5d155f64..0a4150a9 100644 --- a/src/utility.jl +++ b/src/utility.jl @@ -1,3 +1,4 @@ +using ITensors: OpSum """ Relabel sites in OpSum according to given site map """ diff --git a/test/test_additensornetworks.jl b/test/test_additensornetworks.jl index f127df77..8a9965f7 100644 --- a/test/test_additensornetworks.jl +++ b/test/test_additensornetworks.jl @@ -1,7 +1,7 @@ @eval module $(gensym()) using Graphs: rem_edge!, vertices using NamedGraphs: NamedEdge, hexagonal_lattice_graph, named_grid -using ITensorNetworks: ITensorNetwork, inner_network, randomITensorNetwork, siteinds +using ITensorNetworks: ITensorNetwork, inner_network, random_tensornetwork, siteinds using ITensors: ITensors, apply, op, scalar, inner using LinearAlgebra: norm_sqr using Random: Random diff --git a/test/test_forms.jl b/test/test_forms.jl index 3cdb041b..4ff3e2d0 100644 --- a/test/test_forms.jl +++ b/test/test_forms.jl @@ -31,8 +31,8 @@ using Random: Random s_operator = union_all_inds(s, prime(s)) χ, D = 2, 3 Random.seed!(1234) - ψket = random_tensornetwork(s_ket; link_space=χ) - ψbra = random_tensornetwork(s_bra; link_space=χ) + ψket = random_tensornetwork(s; link_space=χ) + ψbra = random_tensornetwork(s; link_space=χ) A = random_tensornetwork(s_operator; link_space=D) blf = BilinearFormNetwork(A, ψbra, ψket) diff --git a/test/test_inner.jl b/test/test_inner.jl index 976f5d3d..9570fb06 100644 --- a/test/test_inner.jl +++ b/test/test_inner.jl @@ -2,7 +2,15 @@ using Test using ITensorNetworks using ITensorNetworks: - TTN, inner, inner_network, loginner, logscalar, randomITensorNetwork, scalar + inner, + inner_network, + loginner, + logscalar, + random_tensornetwork, + scalar, + ttn, + underlying_graph +using ITensorNetworks.ModelHamiltonians: heisenberg using ITensors: dag, siteinds using SplitApplyCombine: group using Graphs: SimpleGraph, uniform_tree @@ -15,8 +23,8 @@ using Random: Random χ = 2 g = NamedGraph(SimpleGraph(uniform_tree(L))) s = siteinds("S=1/2", g) - y = randomITensorNetwork(s; link_space=χ) - x = randomITensorNetwork(s; link_space=χ) + y = random_tensornetwork(s; link_space=χ) + x = random_tensornetwork(s; link_space=χ) #First lets do it with the flattened version of the network xy = inner_network(x, y) @@ -38,7 +46,7 @@ using Random: Random @test xy_scalar ≈ xy_scalar_logbp #test contraction of three layers for expectation values - A = ITensorNetwork(TTN(ITensorNetworks.heisenberg(s), s)) + A = ITensorNetwork(ttn(heisenberg(g), s)) xAy_scalar = inner(x, A, y; alg="exact") xAy_scalar_bp = inner(x, A, y; alg="bp") xAy_scalar_logbp = exp(loginner(x, A, y; alg="bp")) From ad21df985861851f083c184c00223352933fb732 Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 5 Apr 2024 11:24:42 -0400 Subject: [PATCH 29/35] Remove contract_bp example. remove inner(x) --- examples/contract_bp.jl | 40 ----------------------------------- src/abstractitensornetwork.jl | 1 - 2 files changed, 41 deletions(-) delete mode 100644 examples/contract_bp.jl diff --git a/examples/contract_bp.jl b/examples/contract_bp.jl deleted file mode 100644 index 17eea650..00000000 --- a/examples/contract_bp.jl +++ /dev/null @@ -1,40 +0,0 @@ -# https://docs.julialang.org/en/v1/base/base/#Base.names -function imported_names(M::Module) - return setdiff(names(M; imported=true), names(M; imported=false)) -end - -function imported_functions(M::Module) - return map(f -> getfield(M, f), imported_names(M)) -end - -using FileUtils: replace_in_files - -function prepend_imports(M::Module) - fs = imported_functions(M) - # foreach(f -> @show((parentmodule(f), f)), fs) - src_path = joinpath(pkgdir(M), "src") - @show src_path - for f in fs - function_name = last(split(string(f), ".")) - module_name = last(split(string(parentmodule(f)), ".")) - annotated_function_name = module_name * "." * function_name - @show function_name, annotated_function_name - replace_in_files( - src_path, - "function $(function_name)(" => "function $(annotated_function_name)("; - recursive=true, - ignore_dirs=[".git"], - showdiffs=false, - ) - replace_in_files( - src_path, - "\n$(function_name)(" => "\n$(annotated_function_name)("; - recursive=true, - ignore_dirs=[".git"], - showdiffs=false, - ) - end -end - -using ITensorNetworks: ITensorNetworks -prepend_imports(ITensorNetworks) diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index b4a4e061..c5847eae 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -738,7 +738,6 @@ function inner_network( return BilinearFormNetwork(A, x, y; kwargs...) end -inner_network(x::AbstractITensorNetwork; kwargs...) = inner_network(x, x; kwargs...) function norm_sqr_network(ψ::AbstractITensorNetwork) return disjoint_union("bra" => dag(prime(ψ; sites=[])), "ket" => ψ) end From 488109536825be397adb5de9e0be009f23e9e9d1 Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 5 Apr 2024 11:29:34 -0400 Subject: [PATCH 30/35] Added ToDO for norm_sqr_network --- src/abstractitensornetwork.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index c5847eae..77fb73b8 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -738,6 +738,9 @@ function inner_network( return BilinearFormNetwork(A, x, y; kwargs...) end +# TODO: We should make this pass to inner_network and then to BiLinearForm. +# Parts of the code (tests relying on norm_sqr being two layer and the gauging code +# which relies on specific message tensors) currently would break in that case so we need to resolve function norm_sqr_network(ψ::AbstractITensorNetwork) return disjoint_union("bra" => dag(prime(ψ; sites=[])), "ket" => ψ) end From 825d32423fea62af6b02d86f565d65d0efc6e048 Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 5 Apr 2024 11:35:13 -0400 Subject: [PATCH 31/35] Optional positional argument for scalars(bp_cache) --- src/caches/beliefpropagationcache.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index fbab6c67..e8b07c42 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -276,16 +276,20 @@ function region_scalar(bp_cache::BeliefPropagationCache, pe::PartitionEdge) return scalar(vcat(message(bp_cache, pe), message(bp_cache, reverse(pe)))) end -function vertex_scalars(bp_cache::BeliefPropagationCache, pvs::Vector) +function vertex_scalars( + bp_cache::BeliefPropagationCache, + pvs::Vector=partitionvertices(partitioned_tensornetwork(bp_cache)), +) return [region_scalar(bp_cache, pv) for pv in pvs] end -function edge_scalars(bp_cache::BeliefPropagationCache, pes::Vector) +function edge_scalars( + bp_cache::BeliefPropagationCache, + pes::Vector=partitionedges(partitioned_tensornetwork(bp_cache)), +) return [region_scalar(bp_cache, pe) for pe in pes] end function scalar_factors_quotient(bp_cache::BeliefPropagationCache) - pvs = partitionvertices(partitioned_tensornetwork(bp_cache)) - pes = partitionedges(partitioned_tensornetwork(bp_cache)) - return vertex_scalars(bp_cache, pvs), edge_scalars(bp_cache, pes) + return vertex_scalars(bp_cache), edge_scalars(bp_cache) end From 6c91629c286d6b78b32b0c0fd45b7f4f036f02bb Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 5 Apr 2024 12:56:21 -0400 Subject: [PATCH 32/35] Improved, more general message_diff function --- src/caches/beliefpropagationcache.jl | 15 +- test/test_belief_propagation.jl | 201 ++++++++++++++------------- test/test_gauging.jl | 8 +- 3 files changed, 114 insertions(+), 110 deletions(-) diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index e8b07c42..0f06effd 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -10,6 +10,7 @@ using NamedGraphs: boundary_partitionedges, partitionvertices, partitionedges default_message(inds_e) = ITensor[denseblocks(delta(inds_e))] default_messages(ptn::PartitionedGraph) = Dictionary() +default_message_norm(m::ITensor) = norm(m) function default_message_update(contract_list::Vector{ITensor}; kwargs...) sequence = optimal_contraction_sequence(contract_list) updated_messages = contract(contract_list; sequence, kwargs...) @@ -29,12 +30,12 @@ function default_cache_construction_kwargs(alg::Algorithm"bp", ψ::AbstractITens return (; partitioned_vertices=default_partitioned_vertices(ψ)) end -#TODO: Define a version of this that works for QN supporting tensors -function message_diff(message_a::Vector{ITensor}, message_b::Vector{ITensor}) +function message_diff( + message_a::Vector{ITensor}, message_b::Vector{ITensor}; message_norm=default_message_norm +) lhs, rhs = contract(message_a), contract(message_b) - tr_lhs = length(inds(lhs)) == 1 ? sum(lhs) : sum(diag(lhs)) - tr_rhs = length(inds(rhs)) == 1 ? sum(rhs) : sum(diag(rhs)) - return 0.5 * norm((denseblocks(lhs) / tr_lhs) - (denseblocks(rhs) / tr_rhs)) + norm_lhs, norm_rhs = message_norm(lhs), message_norm(rhs) + return 0.5 * norm((denseblocks(lhs) / norm_lhs) - (denseblocks(rhs) / norm_rhs)) end struct BeliefPropagationCache{PTN,MTS,DM} @@ -229,12 +230,12 @@ function update( verbose=false, kwargs..., ) - compute_error = !isnothing(tol) && !hasqns(tensornetwork(bp_cache)) - diff = compute_error ? Ref(0.0) : nothing + compute_error = !isnothing(tol) if isnothing(maxiter) error("You need to specify a number of iterations for BP!") end for i in 1:maxiter + diff = compute_error ? Ref(0.0) : nothing bp_cache = update(bp_cache, edges; (update_diff!)=diff, kwargs...) if compute_error && (diff.x / length(edges)) <= tol if verbose diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index b186bb40..f26f3f2e 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -32,136 +32,137 @@ using Test: @test, @testset @testset "belief_propagation" begin ITensors.disable_warn_order() - #First test on an MPS, should be exact - g_dims = (1, 6) - g = named_grid(g_dims) - s = siteinds("S=1/2", g) - χ = 4 - Random.seed!(1234) - ψ = random_tensornetwork(s; link_space=χ) + # #First test on an MPS, should be exact + # g_dims = (1, 6) + # g = named_grid(g_dims) + # s = siteinds("S=1/2", g) + # χ = 4 + # Random.seed!(1234) + # ψ = random_tensornetwork(s; link_space=χ) - ψψ = ψ ⊗ prime(dag(ψ); sites=[]) + # ψψ = ψ ⊗ prime(dag(ψ); sites=[]) - v = (1, 3) + # v = (1, 3) - Oψ = copy(ψ) - Oψ[v] = apply(op("Sz", s[v]), ψ[v]) - exact_sz = inner(Oψ, ψ) / inner(ψ, ψ) + # Oψ = copy(ψ) + # Oψ[v] = apply(op("Sz", s[v]), ψ[v]) + # exact_sz = inner(Oψ, ψ) / inner(ψ, ψ) - bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) - bpc = update(bpc) - env_tensors = environment(bpc, [PartitionVertex(v)]) - numerator = contract(vcat(env_tensors, ITensor[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))]))[] - denominator = contract(vcat(env_tensors, ITensor[ψ[v], op("I", s[v]), dag(prime(ψ[v]))]))[] + # bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) + # bpc = update(bpc) + # env_tensors = environment(bpc, [PartitionVertex(v)]) + # numerator = contract(vcat(env_tensors, ITensor[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))]))[] + # denominator = contract(vcat(env_tensors, ITensor[ψ[v], op("I", s[v]), dag(prime(ψ[v]))]))[] - @test abs.((numerator / denominator) - exact_sz) <= 1e-14 + # @test abs.((numerator / denominator) - exact_sz) <= 1e-14 - #Test updating the underlying tensornetwork in the cache - v = first(vertices(ψψ)) - new_tensor = randomITensor(inds(ψψ[v])) - bpc = update_factor(bpc, v, new_tensor) - ψψ_updated = tensornetwork(bpc) - @test ψψ_updated[v] == new_tensor + # #Test updating the underlying tensornetwork in the cache + # v = first(vertices(ψψ)) + # new_tensor = randomITensor(inds(ψψ[v])) + # bpc = update_factor(bpc, v, new_tensor) + # ψψ_updated = tensornetwork(bpc) + # @test ψψ_updated[v] == new_tensor - #Now test on a tree, should also be exact - g = named_comb_tree((4, 4)) - s = siteinds("S=1/2", g) - χ = 2 - Random.seed!(1564) - ψ = random_tensornetwork(s; link_space=χ) + # #Now test on a tree, should also be exact + # g = named_comb_tree((4, 4)) + # s = siteinds("S=1/2", g) + # χ = 2 + # Random.seed!(1564) + # ψ = random_tensornetwork(s; link_space=χ) - ψψ = ψ ⊗ prime(dag(ψ); sites=[]) + # ψψ = ψ ⊗ prime(dag(ψ); sites=[]) - v = (1, 3) + # v = (1, 3) - Oψ = copy(ψ) - Oψ[v] = apply(op("Sz", s[v]), ψ[v]) - exact_sz = inner(Oψ, ψ) / inner(ψ, ψ) + # Oψ = copy(ψ) + # Oψ[v] = apply(op("Sz", s[v]), ψ[v]) + # exact_sz = inner(Oψ, ψ) / inner(ψ, ψ) - bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) - bpc = update(bpc) - env_tensors = environment(bpc, [PartitionVertex(v)]) - numerator = contract(vcat(env_tensors, ITensor[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))]))[] - denominator = contract(vcat(env_tensors, ITensor[ψ[v], op("I", s[v]), dag(prime(ψ[v]))]))[] + # bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) + # bpc = update(bpc) + # env_tensors = environment(bpc, [PartitionVertex(v)]) + # numerator = contract(vcat(env_tensors, ITensor[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))]))[] + # denominator = contract(vcat(env_tensors, ITensor[ψ[v], op("I", s[v]), dag(prime(ψ[v]))]))[] - @test abs.((numerator / denominator) - exact_sz) <= 1e-14 + # @test abs.((numerator / denominator) - exact_sz) <= 1e-14 # # #Now test two-site expec taking on the partition function of the Ising model. Not exact, but close g_dims = (3, 4) g = named_grid(g_dims) s = IndsNetwork(g; link_space=2) - beta = 0.2 + beta, h = 0.3, 0.5 vs = [(2, 3), (3, 3)] - ψψ = ModelNetworks.ising_network(s, beta) - ψOψ = ModelNetworks.ising_network(s, beta; szverts=vs) + ψψ = ModelNetworks.ising_network(s, beta; h) + ψOψ = ModelNetworks.ising_network(s, beta; h, szverts=vs) contract_seq = contraction_sequence(ψψ) actual_szsz = contract(ψOψ; sequence=contract_seq)[] / contract(ψψ; sequence=contract_seq)[] - bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) - bpc = update(bpc; maxiter=20) + bpc = BeliefPropagationCache(ψψ, group(v -> v, vertices(ψψ))) + bpc = update(bpc; maxiter=20, verbose=true, tol=1e-5) env_tensors = environment(bpc, vs) numerator = contract(vcat(env_tensors, ITensor[ψOψ[v] for v in vs]))[] denominator = contract(vcat(env_tensors, ITensor[ψψ[v] for v in vs]))[] + @show abs.((numerator / denominator) - actual_szsz) @test abs.((numerator / denominator) - actual_szsz) <= 0.05 # # #Test forming a two-site RDM. Check it has the correct size, trace 1 and is PSD - g_dims = (3, 3) - g = named_grid(g_dims) - s = siteinds("S=1/2", g) - vs = [(2, 2), (2, 3)] - χ = 3 - ψ = random_tensornetwork(s; link_space=χ) - ψψ = ψ ⊗ prime(dag(ψ); sites=[]) - - bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) - bpc = update(bpc; maxiter=20) - - ψψsplit = split_index(ψψ, NamedEdge.([(v, 1) => (v, 2) for v in vs])) - env_tensors = environment(bpc, [(v, 2) for v in vs]) - rdm = contract(vcat(env_tensors, ITensor[ψψsplit[vp] for vp in [(v, 2) for v in vs]])) - - rdm = array((rdm * combiner(inds(rdm; plev=0)...)) * combiner(inds(rdm; plev=1)...)) - rdm /= tr(rdm) - - eigs = eigvals(rdm) - @test size(rdm) == (2^length(vs), 2^length(vs)) - @test all(>=(0), real(eigs)) && all(==(0), imag(eigs)) - - # # #Test more advanced block BP with MPS message tensors on a grid - g_dims = (4, 3) - g = named_grid(g_dims) - s = siteinds("S=1/2", g) - χ = 2 - ψ = random_tensornetwork(s; link_space=χ) - v = (2, 2) - - ψψ = flatten_networks(ψ, dag(ψ); combine_linkinds=false, map_bra_linkinds=prime) - Oψ = copy(ψ) - Oψ[v] = apply(op("Sz", s[v]), ψ[v]) - ψOψ = flatten_networks(ψ, dag(Oψ); combine_linkinds=false, map_bra_linkinds=prime) - - combiners = linkinds_combiners(ψψ) - ψψ = combine_linkinds(ψψ, combiners) - ψOψ = combine_linkinds(ψOψ, combiners) - - bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) - bpc = update( - bpc; - message_update=ITensorNetworks.contract_density_matrix, - message_update_kwargs=(; cutoff=1e-6, maxdim=4), - ) - - env_tensors = environment(bpc, [v]) - numerator = contract(vcat(env_tensors, ITensor[ψOψ[v]]))[] - denominator = contract(vcat(env_tensors, ITensor[ψψ[v]]))[] - - exact_sz = - contract_boundary_mps(ψOψ; cutoff=1e-16) / contract_boundary_mps(ψψ; cutoff=1e-16) - - @test abs.((numerator / denominator) - exact_sz) <= 1e-5 + # g_dims = (3, 3) + # g = named_grid(g_dims) + # s = siteinds("S=1/2", g) + # vs = [(2, 2), (2, 3)] + # χ = 3 + # ψ = random_tensornetwork(s; link_space=χ) + # ψψ = ψ ⊗ prime(dag(ψ); sites=[]) + + # bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) + # bpc = update(bpc; maxiter=20) + + # ψψsplit = split_index(ψψ, NamedEdge.([(v, 1) => (v, 2) for v in vs])) + # env_tensors = environment(bpc, [(v, 2) for v in vs]) + # rdm = contract(vcat(env_tensors, ITensor[ψψsplit[vp] for vp in [(v, 2) for v in vs]])) + + # rdm = array((rdm * combiner(inds(rdm; plev=0)...)) * combiner(inds(rdm; plev=1)...)) + # rdm /= tr(rdm) + + # eigs = eigvals(rdm) + # @test size(rdm) == (2^length(vs), 2^length(vs)) + # @test all(>=(0), real(eigs)) && all(==(0), imag(eigs)) + + # # # #Test more advanced block BP with MPS message tensors on a grid + # g_dims = (4, 3) + # g = named_grid(g_dims) + # s = siteinds("S=1/2", g) + # χ = 2 + # ψ = random_tensornetwork(s; link_space=χ) + # v = (2, 2) + + # ψψ = flatten_networks(ψ, dag(ψ); combine_linkinds=false, map_bra_linkinds=prime) + # Oψ = copy(ψ) + # Oψ[v] = apply(op("Sz", s[v]), ψ[v]) + # ψOψ = flatten_networks(ψ, dag(Oψ); combine_linkinds=false, map_bra_linkinds=prime) + + # combiners = linkinds_combiners(ψψ) + # ψψ = combine_linkinds(ψψ, combiners) + # ψOψ = combine_linkinds(ψOψ, combiners) + + # bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) + # bpc = update( + # bpc; + # message_update=ITensorNetworks.contract_density_matrix, + # message_update_kwargs=(; cutoff=1e-6, maxdim=4), + # ) + + # env_tensors = environment(bpc, [v]) + # numerator = contract(vcat(env_tensors, ITensor[ψOψ[v]]))[] + # denominator = contract(vcat(env_tensors, ITensor[ψψ[v]]))[] + + # exact_sz = + # contract_boundary_mps(ψOψ; cutoff=1e-16) / contract_boundary_mps(ψψ; cutoff=1e-16) + + # @test abs.((numerator / denominator) - exact_sz) <= 1e-5 end end diff --git a/test/test_gauging.jl b/test/test_gauging.jl index b907a710..1c7bff7d 100644 --- a/test/test_gauging.jl +++ b/test/test_gauging.jl @@ -27,8 +27,10 @@ using Test: @test, @testset ψ = random_tensornetwork(s; link_space=χ) # Move directly to vidal gauge - ψ_vidal = VidalITensorNetwork(ψ) - @test gauge_error(ψ_vidal) < 1e-5 + ψ_vidal = VidalITensorNetwork( + ψ; cache_update_kwargs=(; maxiter=20, tol=1e-12, verbose=true) + ) + @test gauge_error(ψ_vidal) < 1e-8 # Move to symmetric gauge cache_ref = Ref{BeliefPropagationCache}() @@ -39,7 +41,7 @@ using Test: @test, @testset @test inner(ψ_symm, ψ) / sqrt(inner(ψ_symm, ψ_symm) * inner(ψ, ψ)) ≈ 1.0 #Test all message tensors are approximately diagonal even when we keep running BP - bp_cache = update(bp_cache; maxiter=20) + bp_cache = update(bp_cache; maxiter=10) for m_e in values(messages(bp_cache)) @test diagITensor(vector(diag(only(m_e))), inds(only(m_e))) ≈ only(m_e) atol = 1e-8 end From fa0f8a77a2b5ece1caf1cba527cdf6e0ba0e3c54 Mon Sep 17 00:00:00 2001 From: Joey Date: Wed, 10 Apr 2024 09:44:10 -0400 Subject: [PATCH 33/35] Removed support for Vector{ITensor} in contract.jl functions for now --- src/ITensorsExtensions/ITensorsExtensions.jl | 82 ++++++++ src/abstractitensornetwork.jl | 3 +- src/caches/beliefpropagationcache.jl | 4 +- src/contract.jl | 33 +--- test/test_belief_propagation.jl | 197 ++++++++++--------- test/test_itensornetwork.jl | 4 +- 6 files changed, 191 insertions(+), 132 deletions(-) create mode 100644 src/ITensorsExtensions/ITensorsExtensions.jl diff --git a/src/ITensorsExtensions/ITensorsExtensions.jl b/src/ITensorsExtensions/ITensorsExtensions.jl new file mode 100644 index 00000000..66350c8f --- /dev/null +++ b/src/ITensorsExtensions/ITensorsExtensions.jl @@ -0,0 +1,82 @@ +module ITensorsExtensions +using LinearAlgebra: LinearAlgebra, eigen, pinv +using ITensors: + ITensor, + Index, + commonind, + dag, + hasqns, + inds, + isdiag, + itensor, + map_diag, + noncommonind, + noprime, + replaceinds, + space, + sqrt_decomp +using ITensors.NDTensors: + NDTensors, + Block, + Tensor, + blockdim, + blockoffsets, + denseblocks, + diaglength, + getdiagindex, + nzblocks, + setdiagindex!, + svd, + tensor, + DiagBlockSparseTensor, + DenseTensor, + BlockOffsets +using Observers: update!, insert_function! + +function NDTensors.blockoffsets(dense::DenseTensor) + return BlockOffsets{ndims(dense)}([Block(ntuple(Returns(1), ndims(dense)))], [0]) +end +function NDTensors.nzblocks(dense::DenseTensor) + return nzblocks(blockoffsets(dense)) +end +NDTensors.blockdim(ind::Int, ::Block{1}) = ind +NDTensors.blockdim(i::Index{Int}, b::Integer) = blockdim(i, Block(b)) +NDTensors.blockdim(i::Index{Int}, b::Block) = blockdim(space(i), b) + +LinearAlgebra.isdiag(it::ITensor) = isdiag(tensor(it)) + +# Convenience functions +sqrt_diag(it::ITensor) = map_diag(sqrt, it) +inv_diag(it::ITensor) = map_diag(inv, it) +invsqrt_diag(it::ITensor) = map_diag(inv ∘ sqrt, it) +pinv_diag(it::ITensor) = map_diag(pinv, it) +pinvsqrt_diag(it::ITensor) = map_diag(pinv ∘ sqrt, it) + +function map_itensor( + f::Function, A::ITensor, lind=first(inds(A)); regularization=nothing, kwargs... +) + USV = svd(A, lind; kwargs...) + U, S, V, spec, u, v = USV + S = map_diag(s -> f(s + regularization), S) + sqrtDL, δᵤᵥ, sqrtDR = sqrt_decomp(S, u, v) + sqrtDR = denseblocks(sqrtDR) * denseblocks(δᵤᵥ) + L, R = U * sqrtDL, V * sqrtDR + return L * R +end + +# Analagous to `denseblocks`. +# Extract the diagonal entries into a diagonal tensor. +function diagblocks(D::Tensor) + nzblocksD = nzblocks(D) + T = DiagBlockSparseTensor(eltype(D), nzblocksD, inds(D)) + for b in nzblocksD + for n in 1:diaglength(D) + setdiagindex!(T, getdiagindex(D, n), n) + end + end + return T +end + +diagblocks(it::ITensor) = itensor(diagblocks(tensor(it))) + +end diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index 77fb73b8..bfcb9397 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -738,9 +738,10 @@ function inner_network( return BilinearFormNetwork(A, x, y; kwargs...) end -# TODO: We should make this pass to inner_network and then to BiLinearForm. +# TODO: We should make this use the QuadraticFormNetwork constructor here. # Parts of the code (tests relying on norm_sqr being two layer and the gauging code # which relies on specific message tensors) currently would break in that case so we need to resolve +# We could have the option in the Form constructors to pre-contract the operator into the bra or ket function norm_sqr_network(ψ::AbstractITensorNetwork) return disjoint_union("bra" => dag(prime(ψ; sites=[])), "ket" => ψ) end diff --git a/src/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index 0f06effd..35769012 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -270,11 +270,11 @@ end function region_scalar(bp_cache::BeliefPropagationCache, pv::PartitionVertex) incoming_mts = environment(bp_cache, [pv]) local_state = factor(bp_cache, pv) - return scalar(vcat(incoming_mts, local_state)) + return contract(vcat(incoming_mts, local_state))[] end function region_scalar(bp_cache::BeliefPropagationCache, pe::PartitionEdge) - return scalar(vcat(message(bp_cache, pe), message(bp_cache, reverse(pe)))) + return contract(vcat(message(bp_cache, pe), message(bp_cache, reverse(pe))))[] end function vertex_scalars( diff --git a/src/contract.jl b/src/contract.jl index c24432c4..a5f3fdd7 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -15,10 +15,6 @@ function NDTensors.contract( return contract(Vector{ITensor}(tn); sequence=sequence_linear_index, kwargs...) end -function NDTensors.contract(alg::Algorithm"exact", tensors::Vector{ITensor}; kwargs...) - return contract(tensors; kwargs...) -end - function NDTensors.contract( alg::Union{Algorithm"density_matrix",Algorithm"ttn_svd"}, tn::AbstractITensorNetwork; @@ -28,40 +24,19 @@ function NDTensors.contract( return approx_tensornetwork(alg, tn, output_structure; kwargs...) end -function contract_density_matrix( - contract_list::Vector{ITensor}; normalize=true, contractor_kwargs... -) - tn, _ = contract( - ITensorNetwork(contract_list); alg="density_matrix", contractor_kwargs... - ) - out = Vector{ITensor}(tn) - if normalize - out .= normalize!.(copy.(out)) - end - return out -end - -function ITensors.scalar( - alg::Algorithm, tn::Union{AbstractITensorNetwork,Vector{ITensor}}; kwargs... -) +function ITensors.scalar(alg::Algorithm, tn::AbstractITensorNetwork; kwargs...) return contract(alg, tn; kwargs...)[] end -function ITensors.scalar( - tn::Union{AbstractITensorNetwork,Vector{ITensor}}; alg="exact", kwargs... -) +function ITensors.scalar(tn::AbstractITensorNetwork; alg="exact", kwargs...) return scalar(Algorithm(alg), tn; kwargs...) end -function logscalar( - tn::Union{AbstractITensorNetwork,Vector{ITensor}}; alg="exact", kwargs... -) +function logscalar(tn::AbstractITensorNetwork; alg="exact", kwargs...) return logscalar(Algorithm(alg), tn; kwargs...) end -function logscalar( - alg::Algorithm"exact", tn::Union{AbstractITensorNetwork,Vector{ITensor}}; kwargs... -) +function logscalar(alg::Algorithm"exact", tn::AbstractITensorNetwork; kwargs...) s = scalar(alg, tn; kwargs...) s = real(s) < 0 ? complex(s) : s return log(s) diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index f26f3f2e..6c380852 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -5,11 +5,13 @@ using ITensorNetworks: ITensorNetworks, BeliefPropagationCache, IndsNetwork, + ITensorNetwork, ⊗, apply, combine_linkinds, contract, contract_boundary_mps, + contract_density_matrix, contraction_sequence, environment, flatten_networks, @@ -32,61 +34,61 @@ using Test: @test, @testset @testset "belief_propagation" begin ITensors.disable_warn_order() - # #First test on an MPS, should be exact - # g_dims = (1, 6) - # g = named_grid(g_dims) - # s = siteinds("S=1/2", g) - # χ = 4 - # Random.seed!(1234) - # ψ = random_tensornetwork(s; link_space=χ) + #First test on an MPS, should be exact + g_dims = (1, 6) + g = named_grid(g_dims) + s = siteinds("S=1/2", g) + χ = 4 + Random.seed!(1234) + ψ = random_tensornetwork(s; link_space=χ) - # ψψ = ψ ⊗ prime(dag(ψ); sites=[]) + ψψ = ψ ⊗ prime(dag(ψ); sites=[]) - # v = (1, 3) + v = (1, 3) - # Oψ = copy(ψ) - # Oψ[v] = apply(op("Sz", s[v]), ψ[v]) - # exact_sz = inner(Oψ, ψ) / inner(ψ, ψ) + Oψ = copy(ψ) + Oψ[v] = apply(op("Sz", s[v]), ψ[v]) + exact_sz = inner(Oψ, ψ) / inner(ψ, ψ) - # bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) - # bpc = update(bpc) - # env_tensors = environment(bpc, [PartitionVertex(v)]) - # numerator = contract(vcat(env_tensors, ITensor[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))]))[] - # denominator = contract(vcat(env_tensors, ITensor[ψ[v], op("I", s[v]), dag(prime(ψ[v]))]))[] + bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) + bpc = update(bpc) + env_tensors = environment(bpc, [PartitionVertex(v)]) + numerator = contract(vcat(env_tensors, ITensor[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))]))[] + denominator = contract(vcat(env_tensors, ITensor[ψ[v], op("I", s[v]), dag(prime(ψ[v]))]))[] - # @test abs.((numerator / denominator) - exact_sz) <= 1e-14 + @test abs.((numerator / denominator) - exact_sz) <= 1e-14 - # #Test updating the underlying tensornetwork in the cache - # v = first(vertices(ψψ)) - # new_tensor = randomITensor(inds(ψψ[v])) - # bpc = update_factor(bpc, v, new_tensor) - # ψψ_updated = tensornetwork(bpc) - # @test ψψ_updated[v] == new_tensor + #Test updating the underlying tensornetwork in the cache + v = first(vertices(ψψ)) + new_tensor = randomITensor(inds(ψψ[v])) + bpc = update_factor(bpc, v, new_tensor) + ψψ_updated = tensornetwork(bpc) + @test ψψ_updated[v] == new_tensor - # #Now test on a tree, should also be exact - # g = named_comb_tree((4, 4)) - # s = siteinds("S=1/2", g) - # χ = 2 - # Random.seed!(1564) - # ψ = random_tensornetwork(s; link_space=χ) + #Now test on a tree, should also be exact + g = named_comb_tree((4, 4)) + s = siteinds("S=1/2", g) + χ = 2 + Random.seed!(1564) + ψ = random_tensornetwork(s; link_space=χ) - # ψψ = ψ ⊗ prime(dag(ψ); sites=[]) + ψψ = ψ ⊗ prime(dag(ψ); sites=[]) - # v = (1, 3) + v = (1, 3) - # Oψ = copy(ψ) - # Oψ[v] = apply(op("Sz", s[v]), ψ[v]) - # exact_sz = inner(Oψ, ψ) / inner(ψ, ψ) + Oψ = copy(ψ) + Oψ[v] = apply(op("Sz", s[v]), ψ[v]) + exact_sz = inner(Oψ, ψ) / inner(ψ, ψ) - # bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) - # bpc = update(bpc) - # env_tensors = environment(bpc, [PartitionVertex(v)]) - # numerator = contract(vcat(env_tensors, ITensor[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))]))[] - # denominator = contract(vcat(env_tensors, ITensor[ψ[v], op("I", s[v]), dag(prime(ψ[v]))]))[] + bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) + bpc = update(bpc) + env_tensors = environment(bpc, [PartitionVertex(v)]) + numerator = contract(vcat(env_tensors, ITensor[ψ[v], op("Sz", s[v]), dag(prime(ψ[v]))]))[] + denominator = contract(vcat(env_tensors, ITensor[ψ[v], op("I", s[v]), dag(prime(ψ[v]))]))[] - # @test abs.((numerator / denominator) - exact_sz) <= 1e-14 + @test abs.((numerator / denominator) - exact_sz) <= 1e-14 - # # #Now test two-site expec taking on the partition function of the Ising model. Not exact, but close + #Now test two-site expec taking on the partition function of the Ising model. Not exact, but close g_dims = (3, 4) g = named_grid(g_dims) s = IndsNetwork(g; link_space=2) @@ -106,63 +108,62 @@ using Test: @test, @testset numerator = contract(vcat(env_tensors, ITensor[ψOψ[v] for v in vs]))[] denominator = contract(vcat(env_tensors, ITensor[ψψ[v] for v in vs]))[] - @show abs.((numerator / denominator) - actual_szsz) @test abs.((numerator / denominator) - actual_szsz) <= 0.05 - # # #Test forming a two-site RDM. Check it has the correct size, trace 1 and is PSD - # g_dims = (3, 3) - # g = named_grid(g_dims) - # s = siteinds("S=1/2", g) - # vs = [(2, 2), (2, 3)] - # χ = 3 - # ψ = random_tensornetwork(s; link_space=χ) - # ψψ = ψ ⊗ prime(dag(ψ); sites=[]) - - # bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) - # bpc = update(bpc; maxiter=20) - - # ψψsplit = split_index(ψψ, NamedEdge.([(v, 1) => (v, 2) for v in vs])) - # env_tensors = environment(bpc, [(v, 2) for v in vs]) - # rdm = contract(vcat(env_tensors, ITensor[ψψsplit[vp] for vp in [(v, 2) for v in vs]])) - - # rdm = array((rdm * combiner(inds(rdm; plev=0)...)) * combiner(inds(rdm; plev=1)...)) - # rdm /= tr(rdm) - - # eigs = eigvals(rdm) - # @test size(rdm) == (2^length(vs), 2^length(vs)) - # @test all(>=(0), real(eigs)) && all(==(0), imag(eigs)) - - # # # #Test more advanced block BP with MPS message tensors on a grid - # g_dims = (4, 3) - # g = named_grid(g_dims) - # s = siteinds("S=1/2", g) - # χ = 2 - # ψ = random_tensornetwork(s; link_space=χ) - # v = (2, 2) - - # ψψ = flatten_networks(ψ, dag(ψ); combine_linkinds=false, map_bra_linkinds=prime) - # Oψ = copy(ψ) - # Oψ[v] = apply(op("Sz", s[v]), ψ[v]) - # ψOψ = flatten_networks(ψ, dag(Oψ); combine_linkinds=false, map_bra_linkinds=prime) - - # combiners = linkinds_combiners(ψψ) - # ψψ = combine_linkinds(ψψ, combiners) - # ψOψ = combine_linkinds(ψOψ, combiners) - - # bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) - # bpc = update( - # bpc; - # message_update=ITensorNetworks.contract_density_matrix, - # message_update_kwargs=(; cutoff=1e-6, maxdim=4), - # ) - - # env_tensors = environment(bpc, [v]) - # numerator = contract(vcat(env_tensors, ITensor[ψOψ[v]]))[] - # denominator = contract(vcat(env_tensors, ITensor[ψψ[v]]))[] - - # exact_sz = - # contract_boundary_mps(ψOψ; cutoff=1e-16) / contract_boundary_mps(ψψ; cutoff=1e-16) - - # @test abs.((numerator / denominator) - exact_sz) <= 1e-5 + #Test forming a two-site RDM. Check it has the correct size, trace 1 and is PSD + g_dims = (3, 3) + g = named_grid(g_dims) + s = siteinds("S=1/2", g) + vs = [(2, 2), (2, 3)] + χ = 3 + ψ = random_tensornetwork(s; link_space=χ) + ψψ = ψ ⊗ prime(dag(ψ); sites=[]) + + bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) + bpc = update(bpc; maxiter=20) + + ψψsplit = split_index(ψψ, NamedEdge.([(v, 1) => (v, 2) for v in vs])) + env_tensors = environment(bpc, [(v, 2) for v in vs]) + rdm = contract(vcat(env_tensors, ITensor[ψψsplit[vp] for vp in [(v, 2) for v in vs]])) + + rdm = array((rdm * combiner(inds(rdm; plev=0)...)) * combiner(inds(rdm; plev=1)...)) + rdm /= tr(rdm) + + eigs = eigvals(rdm) + @test size(rdm) == (2^length(vs), 2^length(vs)) + @test all(>=(0), real(eigs)) && all(==(0), imag(eigs)) + + #Test more advanced block BP with MPS message tensors on a grid + g_dims = (4, 3) + g = named_grid(g_dims) + s = siteinds("S=1/2", g) + χ = 2 + ψ = random_tensornetwork(s; link_space=χ) + v = (2, 2) + + ψψ = flatten_networks(ψ, dag(ψ); combine_linkinds=false, map_bra_linkinds=prime) + Oψ = copy(ψ) + Oψ[v] = apply(op("Sz", s[v]), ψ[v]) + ψOψ = flatten_networks(ψ, dag(Oψ); combine_linkinds=false, map_bra_linkinds=prime) + + combiners = linkinds_combiners(ψψ) + ψψ = combine_linkinds(ψψ, combiners) + ψOψ = combine_linkinds(ψOψ, combiners) + + bpc = BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) + message_update_func(tns; kwargs...) = + Vector{ITensor}(first(contract(ITensorNetwork(tns); alg="density_matrix", kwargs...))) + bpc = update( + bpc; message_update=message_update_func, message_update_kwargs=(; cutoff=1e-6, maxdim=4) + ) + + env_tensors = environment(bpc, [v]) + numerator = contract(vcat(env_tensors, ITensor[ψOψ[v]]))[] + denominator = contract(vcat(env_tensors, ITensor[ψψ[v]]))[] + + exact_sz = + contract_boundary_mps(ψOψ; cutoff=1e-16) / contract_boundary_mps(ψψ; cutoff=1e-16) + + @test abs.((numerator / denominator) - exact_sz) <= 1e-5 end end diff --git a/test/test_itensornetwork.jl b/test/test_itensornetwork.jl index ab67d9b7..93189137 100644 --- a/test/test_itensornetwork.jl +++ b/test/test_itensornetwork.jl @@ -41,6 +41,8 @@ using ITensorNetworks: inner_network, internalinds, linkinds, + neighbor_itensors, + norm_sqr, norm_sqr_network, orthogonalize, random_tensornetwork, @@ -50,8 +52,6 @@ using NamedGraphs: NamedEdge, incident_edges, named_comb_tree, named_grid using Random: Random, randn! using Test: @test, @test_broken, @testset -using ITensorNetworks: norm_sqr, neighbor_itensors - @testset "ITensorNetwork tests" begin @testset "ITensorNetwork Basics" begin Random.seed!(1234) From fd5b9a3a358ccaebff8ea21a8e1a792434c8b9a6 Mon Sep 17 00:00:00 2001 From: Joey Date: Wed, 10 Apr 2024 11:04:55 -0400 Subject: [PATCH 34/35] Fix Bug in test BP --- src/abstractitensornetwork.jl | 1 - test/test_belief_propagation.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index bfcb9397..7a4d2925 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -741,7 +741,6 @@ end # TODO: We should make this use the QuadraticFormNetwork constructor here. # Parts of the code (tests relying on norm_sqr being two layer and the gauging code # which relies on specific message tensors) currently would break in that case so we need to resolve -# We could have the option in the Form constructors to pre-contract the operator into the bra or ket function norm_sqr_network(ψ::AbstractITensorNetwork) return disjoint_union("bra" => dag(prime(ψ; sites=[])), "ket" => ψ) end diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index 6c380852..20b73fdc 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -11,7 +11,6 @@ using ITensorNetworks: combine_linkinds, contract, contract_boundary_mps, - contract_density_matrix, contraction_sequence, environment, flatten_networks, From 002662360b65994db12d3745275ec8a8254505fd Mon Sep 17 00:00:00 2001 From: Joey Date: Fri, 12 Apr 2024 11:10:18 -0400 Subject: [PATCH 35/35] Construct correct identity_network in BiLinearFormNetwork --- src/formnetworks/bilinearformnetwork.jl | 2 +- src/itensornetwork.jl | 12 ++++++++---- test/test_forms.jl | 1 - 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/formnetworks/bilinearformnetwork.jl b/src/formnetworks/bilinearformnetwork.jl index b7408e4c..bcb59704 100644 --- a/src/formnetworks/bilinearformnetwork.jl +++ b/src/formnetworks/bilinearformnetwork.jl @@ -59,7 +59,7 @@ function BilinearFormNetwork( ) @assert issetequal(externalinds(bra), externalinds(ket)) operator_inds = union_all_inds(siteinds(ket), dual_site_index_map(siteinds(ket))) - O = delta_network(operator_inds) + O = ITensorNetwork(Op("I"), operator_inds) return BilinearFormNetwork(O, bra, ket; dual_site_index_map, kwargs...) end diff --git a/src/itensornetwork.jl b/src/itensornetwork.jl index 96f4b604..a2d4af3f 100644 --- a/src/itensornetwork.jl +++ b/src/itensornetwork.jl @@ -164,10 +164,14 @@ function generic_state(a::AbstractArray, inds::Vector) end function generic_state(x::Op, inds::NamedTuple) # TODO: Figure out what to do if there is more than one site. - @assert length(inds.siteinds) == 2 - i = inds.siteinds[findfirst(i -> plev(i) == 0, inds.siteinds)] - @assert i' ∈ inds.siteinds - site_tensors = [op(x.which_op, i)] + if !isempty(inds.siteinds) + @assert length(inds.siteinds) == 2 + i = inds.siteinds[findfirst(i -> plev(i) == 0, inds.siteinds)] + @assert i' ∈ inds.siteinds + site_tensors = [op(x.which_op, i)] + else + site_tensors = [] + end link_tensors = [[onehot(i => 1) for i in inds.linkinds[e]] for e in keys(inds.linkinds)] return contract(reduce(vcat, link_tensors; init=site_tensors)) end diff --git a/test/test_forms.jl b/test/test_forms.jl index 4ff3e2d0..a9a0e453 100644 --- a/test/test_forms.jl +++ b/test/test_forms.jl @@ -8,7 +8,6 @@ using ITensorNetworks: QuadraticFormNetwork, bra_network, bra_vertex, - delta_network, dual_index_map, environment, externalinds,