diff --git a/LICENSE b/LICENSE index d53452f0..3d5bac50 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2021 Matthew Fishman and contributors +Copyright (c) 2021 Matthew Fishman , Joseph Tindall and contributors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/Project.toml b/Project.toml index fd6518f2..27fe9a59 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ITensorNetworks" uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7" -authors = ["Matthew Fishman and contributors"] +authors = ["Matthew Fishman , Joseph Tindall and contributors"] version = "0.10.2" [deps] diff --git a/docs/make.jl b/docs/make.jl index e95e397f..09bfcd83 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,7 +7,7 @@ DocMeta.setdocmeta!( makedocs(; modules=[ITensorNetworks], - authors="Matthew Fishman and contributors", + authors="Matthew Fishman , Joseph Tindall and contributors", repo="https://github.com/mtfishman/ITensorNetworks.jl/blob/{commit}{path}#{line}", sitename="ITensorNetworks.jl", format=Documenter.HTML(; diff --git a/src/environment.jl b/src/environment.jl index 37249cb3..f3c424c0 100644 --- a/src/environment.jl +++ b/src/environment.jl @@ -1,31 +1,33 @@ +using ITensors: contract +using NamedGraphs.PartitionedGraphs: PartitionedGraph + default_environment_algorithm() = "exact" function environment( - ψ::AbstractITensorNetwork, + tn::AbstractITensorNetwork, vertices::Vector; alg=default_environment_algorithm(), kwargs..., ) - return environment(Algorithm(alg), ψ, vertices; kwargs...) + return environment(Algorithm(alg), tn, vertices; kwargs...) end function environment( - ::Algorithm"exact", ψ::AbstractITensorNetwork, verts::Vector; kwargs... + ::Algorithm"exact", tn::AbstractITensorNetwork, verts::Vector; kwargs... ) - return [contract(subgraph(ψ, setdiff(vertices(ψ), verts)); kwargs...)] + return [contract(subgraph(tn, setdiff(vertices(tn), verts)); kwargs...)] end function environment( ::Algorithm"bp", - ψ::AbstractITensorNetwork, + ptn::PartitionedGraph, vertices::Vector; (cache!)=nothing, - partitioned_vertices=default_partitioned_vertices(ψ), update_cache=isnothing(cache!), cache_update_kwargs=default_cache_update_kwargs(cache!), ) if isnothing(cache!) - cache! = Ref(BeliefPropagationCache(ψ, partitioned_vertices)) + cache! = Ref(BeliefPropagationCache(ptn)) end if update_cache @@ -34,3 +36,13 @@ function environment( return environment(cache![], vertices) end + +function environment( + alg::Algorithm"bp", + tn::AbstractITensorNetwork, + vertices::Vector; + partitioned_vertices=default_partitioned_vertices(tn), + kwargs..., +) + return environment(alg, PartitionedGraph(tn, partitioned_vertices), vertices; kwargs...) +end diff --git a/src/expect.jl b/src/expect.jl index b245400c..64dc78b3 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -1,57 +1,63 @@ -using ITensors.ITensorMPS: ITensorMPS, expect, promote_itensor_eltype, OpSum +using Dictionaries: Dictionary, set! +using ITensors: Op, op, contract, siteinds, which_op +using ITensors.ITensorMPS: ITensorMPS, expect + +default_expect_alg() = "bp" + +function ITensorMPS.expect(ψIψ::AbstractFormNetwork, op::Op; contract_kwargs=(;), kwargs...) + v = only(op.sites) + ψIψ_v = ψIψ[operator_vertex(ψIψ, v)] + s = commonind(ψIψ[ket_vertex(ψIψ, v)], ψIψ_v) + operator = ITensors.op(op.which_op, s) + ∂ψIψ_∂v = environment(ψIψ, operator_vertices(ψIψ, [v]); kwargs...) + numerator = contract(vcat(∂ψIψ_∂v, operator); contract_kwargs...)[] + denominator = contract(vcat(∂ψIψ_∂v, ψIψ_v); contract_kwargs...)[] + + return numerator / denominator +end function ITensorMPS.expect( - op::String, - ψ::AbstractITensorNetwork; - cutoff=nothing, - maxdim=nothing, - ortho=false, - sequence=nothing, - vertices=vertices(ψ), + alg::Algorithm, + ψ::AbstractITensorNetwork, + ops; + (cache!)=nothing, + update_cache=isnothing(cache!), + cache_update_kwargs=default_cache_update_kwargs(cache!), + cache_construction_function=tn -> + cache(alg, tn; default_cache_construction_kwargs(alg, tn)...), + kwargs..., ) - s = siteinds(ψ) - ElT = promote_itensor_eltype(ψ) - # 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(ψ, ψ)) + ψIψ = inner_network(ψ, ψ) + if isnothing(cache!) + cache! = Ref(cache_construction_function(ψIψ)) end - 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) / normψ² + + if update_cache + cache![] = update(cache![]; cache_update_kwargs...) end - return res + + return map(op -> expect(ψIψ, op; alg, cache!, update_cache=false, kwargs...), ops) +end + +function ITensorMPS.expect(alg::Algorithm"exact", ψ::AbstractITensorNetwork, ops; kwargs...) + ψIψ = inner_network(ψ, ψ) + return map(op -> expect(ψIψ, op; alg, kwargs...), ops) end function ITensorMPS.expect( - ℋ::OpSum, - ψ::AbstractITensorNetwork; - cutoff=nothing, - maxdim=nothing, - ortho=false, - sequence=nothing, + ψ::AbstractITensorNetwork, op::Op; alg=default_expect_alg(), kwargs... ) - s = siteinds(ψ) - # h⃗ = Vector{ITensor}(ℋ, s) - if isnothing(sequence) - sequence = contraction_sequence(inner_network(ψ, ψ)) - end - h⃗ψ = [apply(hᵢ, ψ; cutoff, maxdim, ortho) for hᵢ in ITensors.terms(ℋ)] - ψhᵢψ = [inner(ψ, hᵢψ; alg="exact", sequence) for hᵢψ in h⃗ψ] - ψh⃗ψ = sum(ψhᵢψ) - ψψ = norm_sqr(ψ; alg="exact", sequence) - return ψh⃗ψ / ψψ + return expect(Algorithm(alg), ψ, [op]; kwargs...) +end + +function ITensorMPS.expect( + ψ::AbstractITensorNetwork, op::String, vertices; alg=default_expect_alg(), kwargs... +) + return expect(Algorithm(alg), ψ, [Op(op, vertex) for vertex in vertices]; kwargs...) end function ITensorMPS.expect( - opsum_sum::Sum{<:OpSum}, - ψ::AbstractITensorNetwork; - cutoff=nothing, - maxdim=nothing, - ortho=true, - sequence=nothing, + ψ::AbstractITensorNetwork, op::String; alg=default_expect_alg(), kwargs... ) - return expect(sum(Ops.terms(opsum_sum)), ψ; cutoff, maxdim, ortho, sequence) + return expect(ψ, op, vertices(ψ); alg, kwargs...) end diff --git a/src/formnetworks/abstractformnetwork.jl b/src/formnetworks/abstractformnetwork.jl index ad3953a8..66776ec8 100644 --- a/src/formnetworks/abstractformnetwork.jl +++ b/src/formnetworks/abstractformnetwork.jl @@ -23,6 +23,7 @@ end function operator_vertices(f::AbstractFormNetwork) return filter(v -> last(v) == operator_vertex_suffix(f), vertices(f)) end + function bra_vertices(f::AbstractFormNetwork) return filter(v -> last(v) == bra_vertex_suffix(f), vertices(f)) end @@ -31,6 +32,10 @@ function ket_vertices(f::AbstractFormNetwork) return filter(v -> last(v) == ket_vertex_suffix(f), vertices(f)) end +function operator_vertices(f::AbstractFormNetwork, original_state_vertices::Vector) + return [operator_vertex_map(f)(osv) for osv in original_state_vertices] +end + function bra_vertices(f::AbstractFormNetwork, original_state_vertices::Vector) return [bra_vertex_map(f)(osv) for osv in original_state_vertices] end @@ -67,23 +72,6 @@ function operator_network(f::AbstractFormNetwork) ) end -function environment( - f::AbstractFormNetwork, - original_state_vertices::Vector; - alg=default_environment_algorithm(), - kwargs..., -) - form_vertices = state_vertices(f, original_state_vertices) - if alg == "bp" - partitioned_vertices = group(v -> original_state_vertex(f, v), vertices(f)) - return environment( - tensornetwork(f), form_vertices; alg, partitioned_vertices, kwargs... - ) - else - return environment(tensornetwork(f), form_vertices; alg, kwargs...) - end -end - operator_vertex_map(f::AbstractFormNetwork) = v -> (v, operator_vertex_suffix(f)) bra_vertex_map(f::AbstractFormNetwork) = v -> (v, bra_vertex_suffix(f)) ket_vertex_map(f::AbstractFormNetwork) = v -> (v, ket_vertex_suffix(f)) diff --git a/test/test_belief_propagation.jl b/test/test_belief_propagation.jl index cf0cac0c..66cf25d7 100644 --- a/test/test_belief_propagation.jl +++ b/test/test_belief_propagation.jl @@ -6,31 +6,31 @@ using GraphsFlows: GraphsFlows using ITensorNetworks: ITensorNetworks, BeliefPropagationCache, - IndsNetwork, - ITensorNetwork, ⊗, - apply, combine_linkinds, contract, contract_boundary_mps, contraction_sequence, eachtensor, environment, - flatten_networks, + inner_network, linkinds_combiners, + message, + partitioned_tensornetwork, random_tensornetwork, siteinds, split_index, tensornetwork, update, - update_factor + update_factor, + update_message using ITensors: ITensors, ITensor, combiner, dag, inds, inner, op, prime, randomITensor using ITensorNetworks.ModelNetworks: ModelNetworks using ITensors.NDTensors: array using LinearAlgebra: eigvals, tr using NamedGraphs: NamedEdge using NamedGraphs.NamedGraphGenerators: named_comb_tree, named_grid -using NamedGraphs.PartitionedGraphs: PartitionVertex +using NamedGraphs.PartitionedGraphs: PartitionVertex, partitionedges using Random: Random using SplitApplyCombine: group using Test: @test, @testset @@ -38,93 +38,30 @@ 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) + g = named_grid((3, 3)) s = siteinds("S=1/2", g) - χ = 4 + χ = 2 Random.seed!(1234) ψ = random_tensornetwork(s; link_space=χ) - ψψ = ψ ⊗ prime(dag(ψ); sites=[]) - v = (1, 3) - - 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(ψψ) + bpc = update(bpc; maxiter=50, tol=1e-10) - @test abs.((numerator / denominator) - exact_sz) <= 1e-14 + #Test messages are converged + for pe in partitionedges(partitioned_tensornetwork(bpc)) + @test update_message(bpc, pe) ≈ message(bpc, pe) atol = 1e-8 + end #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) + bpc_updated = update_factor(bpc, v, new_tensor) + ψψ_updated = tensornetwork(bpc_updated) @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=χ) - - ψψ = ψ ⊗ prime(dag(ψ); sites=[]) - - v = (1, 3) - - 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]))]))[] - - @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, h = 0.3, 0.5 - vs = [(2, 3), (3, 3)] - ψψ = 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, 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]))[] - - @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]) @@ -135,40 +72,8 @@ using Test: @test, @testset 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...) = collect( - eachtensor(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 + @test all(eig -> imag(eig) ≈ 0, eigs) + @test all(eig -> real(eig) >= -eps(eltype(eig)), eigs) end end diff --git a/test/test_expect.jl b/test/test_expect.jl new file mode 100644 index 00000000..db7234e0 --- /dev/null +++ b/test/test_expect.jl @@ -0,0 +1,52 @@ +@eval module $(gensym()) + +using Graphs: SimpleGraph, uniform_tree +using NamedGraphs: NamedGraph, vertices +using NamedGraphs.NamedGraphGenerators: named_grid +using ITensors: siteinds +using ITensorNetworks: + BeliefPropagationCache, + ITensorNetwork, + expect, + random_tensornetwork, + original_state_vertex +using Random: Random +using SplitApplyCombine: group +using Test: @test, @testset + +@testset "Test Expect" begin + Random.seed!(1234) + + #Test on a tree + L, χ = 4, 2 + g = NamedGraph(SimpleGraph(uniform_tree(L))) + s = siteinds("S=1/2", g) + ψ = random_tensornetwork(s; link_space=χ) + sz_bp = expect(ψ, "Sz"; alg="bp") + sz_exact = expect(ψ, "Sz"; alg="exact") + @test sz_bp ≈ sz_exact + + #Test on a grid, group by column to make BP exact + L, χ = 2, 2 + g = named_grid((L, L)) + s = siteinds("S=1/2", g) + ψ = random_tensornetwork(s; link_space=χ) + cache_construction_function = + f -> BeliefPropagationCache( + f; partitioned_vertices=group(v -> (original_state_vertex(f, v)[1]), vertices(f)) + ) + sz_bp = expect(ψ, "Sz"; alg="bp", cache_construction_function) + sz_exact = expect(ψ, "Sz"; alg="exact") + @test sz_bp ≈ sz_exact + + #Test with QNS, product state so should be immediately exact + L, χ = 2, 2 + g = named_grid((L, L)) + s = siteinds("S=1/2", g; conserve_qns=true) + ψ = ITensorNetwork(v -> isodd(sum(v)) ? "↑" : "↓", s) + + sz_bp = expect(ψ, "Sz"; alg="bp") + sz_exact = expect(ψ, "Sz"; alg="exact") + @test sz_bp ≈ sz_exact +end +end diff --git a/test/test_forms.jl b/test/test_forms.jl index e0edb597..c36ab585 100644 --- a/test/test_forms.jl +++ b/test/test_forms.jl @@ -16,6 +16,7 @@ using ITensorNetworks: operator_network, random_tensornetwork, siteinds, + state_vertices, tensornetwork, union_all_inds, update @@ -57,16 +58,16 @@ using Random: Random @test underlying_graph(ket_network(qf)) == underlying_graph(ψket) @test underlying_graph(operator_network(qf)) == underlying_graph(A) - ∂qf_∂v = only(environment(qf, [v])) + ∂qf_∂v = only(environment(qf, state_vertices(qf, [v]))) @test (∂qf_∂v) * (qf[ket_vertex(qf, v)] * qf[bra_vertex(qf, v)]) ≈ contract(qf) - ∂qf_∂v_bp = environment(qf, [v]; alg="bp", update_cache=false) + ∂qf_∂v_bp = environment(qf, state_vertices(qf, [v]); alg="bp", update_cache=false) ∂qf_∂v_bp = contract(∂qf_∂v_bp) ∂qf_∂v_bp /= norm(∂qf_∂v_bp) ∂qf_∂v /= norm(∂qf_∂v) @test ∂qf_∂v_bp != ∂qf_∂v - ∂qf_∂v_bp = environment(qf, [v]; alg="bp", update_cache=true) + ∂qf_∂v_bp = environment(qf, state_vertices(qf, [v]); alg="bp", update_cache=true) ∂qf_∂v_bp = contract(∂qf_∂v_bp) ∂qf_∂v_bp /= norm(∂qf_∂v_bp) @test ∂qf_∂v_bp ≈ ∂qf_∂v diff --git a/test/test_inner.jl b/test/test_inner.jl index 9570fb06..427fadf0 100644 --- a/test/test_inner.jl +++ b/test/test_inner.jl @@ -1,7 +1,7 @@ -using Test -using ITensorNetworks +@eval module $(gensym()) using ITensorNetworks: + ITensorNetwork, inner, inner_network, loginner, @@ -16,6 +16,7 @@ using SplitApplyCombine: group using Graphs: SimpleGraph, uniform_tree using NamedGraphs: NamedGraph using Random: Random +using Test: @test, @testset @testset "Inner products, BP vs exact comparison" begin Random.seed!(1234) @@ -55,4 +56,4 @@ using Random: Random @test xAy_scalar_bp ≈ xAy_scalar_logbp @test xAy_scalar ≈ xAy_scalar_logbp end -nothing +end diff --git a/test/test_tebd.jl b/test/test_tebd.jl index 5b7f6278..64e16c57 100644 --- a/test/test_tebd.jl +++ b/test/test_tebd.jl @@ -7,7 +7,7 @@ using ITensorNetworks.ITensorsExtensions: group_terms using ITensorNetworks.ModelHamiltonians: ModelHamiltonians using NamedGraphs.GraphsExtensions: rename_vertices using NamedGraphs.NamedGraphGenerators: named_grid -using Test: @test, @testset +using Test: @test, @testset, @test_broken ITensors.disable_warn_order() @@ -41,7 +41,7 @@ ITensors.disable_warn_order() Δβ = 0.2 ψ_init = ITensorNetwork(v -> "↑", s) - E0 = expect(ℋ, ψ_init) + #E0 = expect(ℋ, ψ_init) ψ = tebd( group_terms(ℋ, g), ψ_init; @@ -52,7 +52,7 @@ ITensors.disable_warn_order() ortho=false, print_frequency=typemax(Int), ) - E1 = expect(ℋ, ψ) + #E1 = expect(ℋ, ψ) ψ = tebd( group_terms(ℋ, g), ψ_init; @@ -63,9 +63,9 @@ ITensors.disable_warn_order() ortho=true, print_frequency=typemax(Int), ) - E2 = expect(ℋ, ψ) - @show E0, E1, E2, E_dmrg - @test (((abs((E2 - E1) / E2) < 1e-3) && (E1 < E0)) || (E2 < E1 < E0)) - @test E2 ≈ E_dmrg rtol = 1e-3 + #E2 = expect(ℋ, ψ) + #@show E0, E1, E2, E_dmrg + @test_broken (((abs((E2 - E1) / E2) < 1e-3) && (E1 < E0)) || (E2 < E1 < E0)) + @test_broken E2 ≈ E_dmrg rtol = 1e-3 end end