From f813653c4eb5be0711986909b4c00716d66fc744 Mon Sep 17 00:00:00 2001 From: Joseph Tindall <51231103+JoeyT1994@users.noreply.github.com> Date: Fri, 12 Apr 2024 12:47:08 -0400 Subject: [PATCH] Implement `inner` using BP (#147) --- README.md | 12 +- src/ITensorNetworks.jl | 3 +- src/ITensorsExtensions/ITensorsExtensions.jl | 82 ++++++++++ src/abstractitensornetwork.jl | 66 ++------ .../binary_tree_partition.jl | 2 +- src/caches/beliefpropagationcache.jl | 59 +++++++- src/contract.jl | 59 ++++++-- src/expect.jl | 12 +- src/formnetworks/bilinearformnetwork.jl | 32 +++- src/formnetworks/quadraticformnetwork.jl | 27 +++- src/gauging.jl | 8 +- src/inner.jl | 143 ++++++++++++++++++ src/itensornetwork.jl | 12 +- src/utility.jl | 1 + test/test_additensornetworks.jl | 28 +--- test/test_apply.jl | 28 ++-- test/test_belief_propagation.jl | 33 ++-- test/test_forms.jl | 10 +- test/test_gauging.jl | 14 +- test/test_inner.jl | 58 +++++++ test/test_itensornetwork.jl | 42 ++--- test/test_tebd.jl | 15 +- test/test_tno.jl | 16 +- 23 files changed, 554 insertions(+), 208 deletions(-) create mode 100644 src/ITensorsExtensions/ITensorsExtensions.jl create mode 100644 src/inner.jl create mode 100644 test/test_inner.jl diff --git a/README.md b/README.md index 2f86bee0..fe933ec8 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=723|"1×1,1×2"), (dim=2|id=712|"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=656|"S=1/2,Site,n=3"), (dim=2|id=190|"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/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 5e58d1b7..bfe4b8a2 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -13,7 +13,6 @@ include("opsum.jl") include("sitetype.jl") include("abstractitensornetwork.jl") include("contraction_sequences.jl") -include("expect.jl") include("tebd.jl") include("itensornetwork.jl") include("mincut.jl") @@ -64,6 +63,8 @@ include("solvers/contract.jl") include("solvers/linsolve.jl") include("solvers/sweep_plans/sweep_plans.jl") include("apply.jl") +include("inner.jl") +include("expect.jl") include("environment.jl") include("exports.jl") include("ModelHamiltonians/ModelHamiltonians.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 e700ed6d..3e044abe 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -735,67 +735,23 @@ 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=[]))) - # TODO: Use a different constructor here? - 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; kwargs...) + return BilinearFormNetwork(x, y; kwargs...) 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..., + x::AbstractITensorNetwork, A::AbstractITensorNetwork, y::AbstractITensorNetwork; 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 + return BilinearFormNetwork(A, x, y; kwargs...) 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)[] +# 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 +function norm_sqr_network(ψ::AbstractITensorNetwork) + return disjoint_union("bra" => dag(prime(ψ; sites=[])), "ket" => ψ) 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) - -norm_sqr_network(ψ::AbstractITensorNetwork; kwargs...) = inner_network(ψ, ψ; kwargs...) - # # Printing # @@ -942,7 +898,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) @@ -966,3 +922,5 @@ function ITensorMPS.add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork end Base.:+(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork) = add(tn1, tn2) + +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 7b37dff5..b6657c19 100644 --- a/src/approx_itensornetwork/binary_tree_partition.jl +++ b/src/approx_itensornetwork/binary_tree_partition.jl @@ -130,6 +130,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/caches/beliefpropagationcache.jl b/src/caches/beliefpropagationcache.jl index 5e7f8e43..35769012 100644 --- a/src/caches/beliefpropagationcache.jl +++ b/src/caches/beliefpropagationcache.jl @@ -6,10 +6,11 @@ using NamedGraphs: PartitionVertex using LinearAlgebra: diag using ITensors: dir 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() +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...) @@ -21,12 +22,20 @@ end return default_bp_maxiter(undirected_graph(underlying_graph(g))) end default_partitioned_vertices(ψ::AbstractITensorNetwork) = group(v -> v, vertices(ψ)) +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 (; partitioned_vertices=default_partitioned_vertices(ψ)) +end -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) - return 0.5 * - norm((denseblocks(lhs) / sum(diag(lhs))) - (denseblocks(rhs) / sum(diag(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} @@ -47,8 +56,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_tensornetwork(bp_cache::BeliefPropagationCache) @@ -118,7 +133,7 @@ function environment( ) 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 environment( @@ -216,11 +231,11 @@ function update( kwargs..., ) compute_error = !isnothing(tol) - diff = compute_error ? Ref(0.0) : nothing 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 @@ -251,3 +266,31 @@ end function update_factor(bp_cache, vertex, factor) return update_factors(bp_cache, [vertex], ITensor[factor]) end + +function region_scalar(bp_cache::BeliefPropagationCache, pv::PartitionVertex) + incoming_mts = environment(bp_cache, [pv]) + local_state = factor(bp_cache, pv) + return contract(vcat(incoming_mts, local_state))[] +end + +function region_scalar(bp_cache::BeliefPropagationCache, pe::PartitionEdge) + return contract(vcat(message(bp_cache, pe), message(bp_cache, reverse(pe))))[] +end + +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=partitionedges(partitioned_tensornetwork(bp_cache)), +) + return [region_scalar(bp_cache, pe) for pe in pes] +end + +function scalar_factors_quotient(bp_cache::BeliefPropagationCache) + return vertex_scalars(bp_cache), edge_scalars(bp_cache) +end diff --git a/src/contract.jl b/src/contract.jl index f358bb57..a5f3fdd7 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -1,10 +1,10 @@ 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! -function NDTensors.contract(tn::AbstractITensorNetwork; alg::String="exact", kwargs...) +function NDTensors.contract(tn::AbstractITensorNetwork; alg="exact", kwargs...) return contract(Algorithm(alg), tn; kwargs...) end @@ -24,15 +24,52 @@ function NDTensors.contract( return approx_tensornetwork(alg, tn, output_structure; kwargs...) end -function contract_density_matrix( - contract_list::Vector{ITensor}; normalize=true, contractor_kwargs... +function ITensors.scalar(alg::Algorithm, tn::AbstractITensorNetwork; kwargs...) + return contract(alg, tn; kwargs...)[] +end + +function ITensors.scalar(tn::AbstractITensorNetwork; alg="exact", kwargs...) + return scalar(Algorithm(alg), tn; kwargs...) +end + +function logscalar(tn::AbstractITensorNetwork; alg="exact", kwargs...) + return logscalar(Algorithm(alg), tn; kwargs...) +end + +function logscalar(alg::Algorithm"exact", tn::AbstractITensorNetwork; kwargs...) + s = scalar(alg, tn; kwargs...) + s = real(s) < 0 ? complex(s) : s + return log(s) +end + +function logscalar( + alg::Algorithm, + tn::AbstractITensorNetwork; + (cache!)=nothing, + cache_construction_kwargs=default_cache_construction_kwargs(alg, tn), + update_cache=isnothing(cache!), + cache_update_kwargs=default_cache_update_kwargs(cache!), ) - tn, _ = contract( - ITensorNetwork(contract_list); alg="density_matrix", contractor_kwargs... - ) - out = Vector{ITensor}(tn) - if normalize - out .= normalize!.(copy.(out)) + if isnothing(cache!) + 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_quotient(cache![]) + 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 out + + return sum(log.(numerator_terms)) - sum(log.((denominator_terms))) +end + +function ITensors.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 5f1432d1..b245400c 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -14,13 +14,13 @@ function ITensorMPS.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(ψ, ψ)) end - normψ² = norm_sqr(ψ; sequence) + normψ² = norm_sqr(ψ; alg="exact", sequence) 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) / normψ² end return res end @@ -36,12 +36,12 @@ function ITensorMPS.expect( s = siteinds(ψ) # h⃗ = Vector{ITensor}(ℋ, s) if isnothing(sequence) - sequence = contraction_sequence(inner_network(ψ, ψ; flatten=true)) + sequence = contraction_sequence(inner_network(ψ, ψ)) 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", sequence) for hᵢψ in h⃗ψ] ψh⃗ψ = sum(ψhᵢψ) - ψψ = norm_sqr(ψ; sequence) + ψψ = norm_sqr(ψ; alg="exact", sequence) return ψh⃗ψ / ψψ end diff --git a/src/formnetworks/bilinearformnetwork.jl b/src/formnetworks/bilinearformnetwork.jl index e55f74b6..bcb59704 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,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, ) + 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 @@ -44,23 +52,31 @@ function Base.copy(blf::BilinearFormNetwork) end 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)) - O = delta_network(operator_inds) - return BilinearFormNetwork(O, bra, ket; kwargs...) + @assert issetequal(externalinds(bra), externalinds(ket)) + operator_inds = union_all_inds(siteinds(ket), dual_site_index_map(siteinds(ket))) + O = ITensorNetwork(Op("I"), operator_inds) + 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 a5dfca5a..d0254501 100644 --- a/src/formnetworks/quadraticformnetwork.jl +++ b/src/formnetworks/quadraticformnetwork.jl @@ -41,8 +41,14 @@ 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 +58,25 @@ 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/gauging.jl b/src/gauging.jl index 449a8cbd..4d2c4f6a 100644 --- a/src/gauging.jl +++ b/src/gauging.jl @@ -25,8 +25,8 @@ function Base.copy(ψ::VidalITensorNetwork) end function default_norm_cache(ψ::ITensorNetwork) - ψψ = norm_network(ψ) - return BeliefPropagationCache(ψψ, group(v -> v[1], vertices(ψψ))) + ψψ = norm_sqr_network(ψ) + return BeliefPropagationCache(ψψ, group(v -> first(v), vertices(ψψ))) end function ITensorNetwork( @@ -51,7 +51,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 @@ -80,7 +80,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 new file mode 100644 index 00000000..166a2c6f --- /dev/null +++ b/src/inner.jl @@ -0,0 +1,143 @@ +using ITensors: inner, scalar, loginner +using LinearAlgebra: norm, norm_sqr + +default_contract_alg(tns::Tuple) = "bp" + +function ITensors.inner( + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + alg=default_contract_alg((ϕ, ψ)), + kwargs..., +) + return inner(Algorithm(alg), ϕ, ψ; kwargs...) +end + +function ITensors.inner( + ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + alg=default_contract_alg((ϕ, A, ψ)), + kwargs..., +) + return inner(Algorithm(alg), ϕ, A, ψ; kwargs...) +end + +function ITensors.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 + +function ITensors.inner( + alg::Algorithm"exact", + ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + sequence=nothing, + contraction_sequence_kwargs=(;), + kwargs..., +) + tn = inner_network(ϕ, A, ψ; kwargs...) + if isnothing(sequence) + sequence = contraction_sequence(tn; contraction_sequence_kwargs...) + end + return scalar(tn; sequence) +end + +function ITensors.loginner( + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + alg=default_contract_alg((ϕ, ψ)), + kwargs..., +) + return loginner(Algorithm(alg), ϕ, ψ; kwargs...) +end + +function ITensors.loginner( + ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + alg=default_contract_alg((ϕ, A, ψ)), + kwargs..., +) + return loginner(Algorithm(alg), ϕ, A, ψ; kwargs...) +end + +function ITensors.loginner( + alg::Algorithm"exact", ϕ::AbstractITensorNetwork, ψ::AbstractITensorNetwork; kwargs... +) + return log(inner(alg, ϕ, ψ); kwargs...) +end + +function ITensors.loginner( + alg::Algorithm"exact", + ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + kwargs..., +) + return log(inner(alg, ϕ, A, ψ); kwargs...) +end + +function ITensors.loginner( + alg::Algorithm"bp", + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + dual_link_index_map=sim, + kwargs..., +) + tn = inner_network(ϕ, ψ; dual_link_index_map) + return logscalar(alg, tn; kwargs...) +end + +function ITensors.loginner( + alg::Algorithm"bp", + ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + dual_link_index_map=sim, + kwargs..., +) + tn = inner_network(ϕ, A, ψ; dual_link_index_map) + return logscalar(alg, tn; kwargs...) +end + +function ITensors.inner( + alg::Algorithm"bp", + ϕ::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + dual_link_index_map=sim, + kwargs..., +) + tn = inner_network(ϕ, ψ; dual_link_index_map) + return scalar(alg, tn; kwargs...) +end + +function ITensors.inner( + alg::Algorithm"bp", + ϕ::AbstractITensorNetwork, + A::AbstractITensorNetwork, + ψ::AbstractITensorNetwork; + dual_link_index_map=sim, + kwargs..., +) + tn = inner_network(ϕ, A, ψ; dual_link_index_map) + return scalar(alg, tn; kwargs...) +end + +# TODO: rename `sqnorm` to match https://github.com/JuliaStats/Distances.jl, +# or `norm_sqr` to match `LinearAlgebra.norm_sqr` +LinearAlgebra.norm_sqr(ψ::AbstractITensorNetwork; kwargs...) = inner(ψ, ψ; kwargs...) + +function LinearAlgebra.norm(ψ::AbstractITensorNetwork; kwargs...) + return sqrt(abs(real(norm_sqr(ψ; 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/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 e2ab2d89..279f3b2f 100644 --- a/test/test_additensornetworks.jl +++ b/test/test_additensornetworks.jl @@ -2,13 +2,14 @@ using Graphs: rem_edge!, vertices using NamedGraphs: NamedEdge, hexagonal_lattice_graph, named_grid using ITensorNetworks: ITensorNetwork, inner_network, random_tensornetwork, siteinds -using ITensors: ITensors, apply, contract, op +using ITensors: ITensors, apply, op, scalar, inner +using LinearAlgebra: norm_sqr using Random: Random using Test: @test, @testset @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(v -> "↑", s) ψ2 = ITensorNetwork(v -> "↓", s) @@ -22,11 +23,9 @@ using Test: @test, @testset ψψ_GHZ = inner_network(ψ_GHZ, ψ_GHZ) ψOψ_GHZ = inner_network(ψ_GHZ, Oψ_GHZ) - @test contract(ψOψ_GHZ)[] / 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))) @@ -46,22 +45,11 @@ using Test: @test, @testset 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 = - (contract(ψOψ_1)[] + contract(ψOψ_2)[] + 2 * contract(ψ1Oψ2)[]) / - (contract(ψψ_1)[] + contract(ψψ_2)[] + 2 * contract(ψ1ψ2)[]) - expec_method2 = contract(ψOψ_12)[] / 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 062cbb2f..d4c408e0 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, random_tensornetwork, siteinds, update -using ITensors: ITensors +using ITensors: ITensors, inner, op using NamedGraphs: PartitionVertex, named_grid using Random: Random using SplitApplyCombine: group @@ -27,8 +26,7 @@ using Test: @test, @testset χ = 2 ψ = random_tensornetwork(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(ψψ))) bp_cache = update(bp_cache; maxiter=20) @@ -39,12 +37,14 @@ using Test: @test, @testset #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 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( @@ -68,14 +68,16 @@ using Test: @test, @testset envisposdef=true, ) fSBP = - contract_inner(ψOSBP, ψOexact) / - sqrt(contract_inner(ψOexact, ψOexact) * contract_inner(ψOSBP, ψOSBP)) + inner(ψOSBP, ψOexact; alg=inner_alg) / + sqrt(inner(ψOexact, ψOexact; alg=inner_alg) * inner(ψOSBP, ψOSBP; alg=inner_alg)) fVidal = - contract_inner(ψOVidal_symm, ψOexact) / - sqrt(contract_inner(ψOexact, ψOexact) * contract_inner(ψOVidal_symm, ψOVidal_symm)) + inner(ψOVidal_symm, ψOexact; alg=inner_alg) / sqrt( + inner(ψOexact, ψOexact; alg=inner_alg) * + inner(ψOVidal_symm, ψOVidal_symm; alg=inner_alg), + ) fGBP = - contract_inner(ψOGBP, ψOexact) / - sqrt(contract_inner(ψOexact, ψOexact) * contract_inner(ψOGBP, ψOGBP)) + 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_belief_propagation.jl b/test/test_belief_propagation.jl index fd029f3a..20b73fdc 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -5,11 +5,11 @@ using ITensorNetworks: ITensorNetworks, BeliefPropagationCache, IndsNetwork, + ITensorNetwork, ⊗, apply, combine_linkinds, contract, - contract_inner, contract_boundary_mps, contraction_sequence, environment, @@ -21,8 +21,8 @@ using ITensorNetworks: tensornetwork, update, update_factor +using ITensors: ITensors, ITensor, combiner, dag, inds, inner, op, prime, randomITensor using ITensorNetworks.ModelNetworks: ModelNetworks -using ITensors: ITensors, ITensor, combiner, dag, inds, op, prime, randomITensor using ITensors.NDTensors: array using LinearAlgebra: eigvals, tr using NamedGraphs: NamedEdge, PartitionVertex, named_comb_tree, named_grid @@ -30,9 +30,8 @@ using Random: Random using SplitApplyCombine: group using Test: @test, @testset -ITensors.disable_warn_order() - @testset "belief_propagation" begin + ITensors.disable_warn_order() #First test on an MPS, should be exact g_dims = (1, 6) @@ -48,7 +47,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) @@ -78,7 +77,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) @@ -88,21 +87,21 @@ ITensors.disable_warn_order() @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) - 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]))[] @@ -110,7 +109,7 @@ ITensors.disable_warn_order() @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 + #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) @@ -133,7 +132,7 @@ ITensors.disable_warn_order() @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 + #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) @@ -151,10 +150,10 @@ ITensors.disable_warn_order() ψ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=ITensorNetworks.contract_density_matrix, - message_update_kwargs=(; cutoff=1e-6, maxdim=4), + bpc; message_update=message_update_func, message_update_kwargs=(; cutoff=1e-6, maxdim=4) ) env_tensors = environment(bpc, [v]) diff --git a/test/test_forms.jl b/test/test_forms.jl index 1d940cfd..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, @@ -27,13 +26,12 @@ using Random: Random @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 = 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_gauging.jl b/test/test_gauging.jl index bd8af9cb..1c7bff7d 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, random_tensornetwork, siteinds, update -using ITensors: diagITensor, inds +using ITensors: diagITensor, inds, inner using ITensors.NDTensors: vector using LinearAlgebra: diag using NamedGraphs: named_grid @@ -28,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}() @@ -37,11 +38,10 @@ using Test: @test, @testset 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) + 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 diff --git a/test/test_inner.jl b/test/test_inner.jl new file mode 100644 index 00000000..9570fb06 --- /dev/null +++ b/test/test_inner.jl @@ -0,0 +1,58 @@ +using Test +using ITensorNetworks + +using ITensorNetworks: + 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 +using NamedGraphs: NamedGraph +using Random: Random + +@testset "Inner products, BP vs exact comparison" begin + Random.seed!(1234) + L = 4 + χ = 2 + g = NamedGraph(SimpleGraph(uniform_tree(L))) + s = siteinds("S=1/2", g) + 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) + xy_scalar = scalar(xy) + xy_scalar_bp = scalar(xy; alg="bp") + xy_scalar_logbp = exp(logscalar(xy; alg="bp")) + + @test xy_scalar ≈ xy_scalar_bp + @test xy_scalar_bp ≈ xy_scalar_logbp + @test xy_scalar ≈ xy_scalar_logbp + + #Now lets 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 = 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")) + + @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 62f53d07..18845dd5 100644 --- a/test/test_itensornetwork.jl +++ b/test/test_itensornetwork.jl @@ -26,6 +26,7 @@ using ITensors: hascommoninds, hasinds, inds, + inner, itensor, onehot, order, @@ -44,6 +45,9 @@ using ITensorNetworks: inner_network, internalinds, linkinds, + neighbor_itensors, + norm_sqr, + norm_sqr_network, orthogonalize, random_tensornetwork, siteinds, @@ -136,10 +140,10 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64}) g = named_grid(dims) s = siteinds("S=1/2", g) ψ = ITensorNetwork(v -> "↑", s) - 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 @@ -148,15 +152,15 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64}) s = siteinds("S=1/2", g) ψ = ITensorNetwork(v -> "↑", s) 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 (eltype=$elt)" for elt in elts, @@ -263,26 +267,26 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64}) @testset "orthogonalize" begin tn = random_tensornetwork(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 @@ -319,7 +323,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64}) 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)) diff --git a/test/test_tebd.jl b/test/test_tebd.jl index 99aa9bbc..fe7185f1 100644 --- a/test/test_tebd.jl +++ b/test/test_tebd.jl @@ -11,7 +11,7 @@ using Test: @test, @testset ITensors.disable_warn_order() @testset "Ising TEBD" begin - dims = (4, 4) + dims = (2, 3) n = prod(dims) g = named_grid(dims) @@ -39,11 +39,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(v -> "↑", s) - E0 = expect(ℋ, ψ_init; sequence=inner_sequence) + E0 = expect(ℋ, ψ_init) ψ = tebd( group_terms(ℋ, g), ψ_init; @@ -54,7 +51,7 @@ ITensors.disable_warn_order() ortho=false, print_frequency=typemax(Int), ) - E1 = expect(ℋ, ψ; sequence=inner_sequence) + E1 = expect(ℋ, ψ) ψ = tebd( group_terms(ℋ, g), ψ_init; @@ -65,9 +62,9 @@ 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 end diff --git a/test/test_tno.jl b/test/test_tno.jl index 1512768c..30811130 100644 --- a/test/test_tno.jl +++ b/test/test_tno.jl @@ -2,15 +2,14 @@ using Graphs: vertices using ITensorNetworks: apply, - contract_inner, flatten_networks, group_commuting_itensors, gate_group_to_tno, get_tnos, random_tensornetwork, siteinds +using ITensors: ITensor, inner, noprime using ITensorNetworks.ModelHamiltonians: ModelHamiltonians -using ITensors: ITensor, noprime using NamedGraphs: named_grid using Test: @test, @testset @@ -53,13 +52,12 @@ using Test: @test, @testset ψ_tno[v] = 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