From 2c29c6cfe0301ce36d50d399981debd342734b22 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Thu, 18 Apr 2024 09:37:06 -0400 Subject: [PATCH] Fix more tests --- src/abstractitensornetwork.jl | 6 ++- src/contract.jl | 6 ++- src/contraction_sequences.jl | 5 ++- src/formnetworks/abstractformnetwork.jl | 1 + src/gauging.jl | 2 +- src/treetensornetworks/opsum_to_ttn.jl | 50 ++++++++++++++++--------- test/test_opsum_to_ttn.jl | 34 ++++++++--------- 7 files changed, 63 insertions(+), 41 deletions(-) diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl index 2dc983ee..634f0ffe 100644 --- a/src/abstractitensornetwork.jl +++ b/src/abstractitensornetwork.jl @@ -38,7 +38,8 @@ using ITensors: swaptags using ITensors.ITensorMPS: ITensorMPS, add, linkdim, linkinds, siteinds using LinearAlgebra: LinearAlgebra, factorize -using NamedGraphs: NamedGraphs, NamedGraph, not_implemented, vertex_to_parent_vertex +using NamedGraphs: + NamedGraphs, NamedGraph, not_implemented, parent_vertex_to_vertex, vertex_to_parent_vertex using NamedGraphs.GraphsExtensions: ⊔, directed_graph, incident_edges, rename_vertices, vertextype using NDTensors: NDTensors, dim @@ -95,6 +96,9 @@ DataGraphs.underlying_graph(tn::AbstractITensorNetwork) = underlying_graph(data_ function NamedGraphs.vertex_to_parent_vertex(tn::AbstractITensorNetwork, vertex) return vertex_to_parent_vertex(underlying_graph(tn), vertex) end +function NamedGraphs.parent_vertex_to_vertex(tn::AbstractITensorNetwork, parent_vertex) + return parent_vertex_to_vertex(underlying_graph(tn), parent_vertex) +end # # Iteration diff --git a/src/contract.jl b/src/contract.jl index a5f3fdd7..309c297c 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -9,7 +9,11 @@ function NDTensors.contract(tn::AbstractITensorNetwork; alg="exact", kwargs...) end function NDTensors.contract( - alg::Algorithm"exact", tn::AbstractITensorNetwork; sequence=vertices(tn), kwargs... + alg::Algorithm"exact", + tn::AbstractITensorNetwork; + contraction_sequence_kwargs=(;), + sequence=contraction_sequence(tn; contraction_sequence_kwargs...), + kwargs..., ) sequence_linear_index = deepmap(v -> vertex_to_parent_vertex(tn, v), sequence) return contract(Vector{ITensor}(tn); sequence=sequence_linear_index, kwargs...) diff --git a/src/contraction_sequences.jl b/src/contraction_sequences.jl index e60575b4..eb10d465 100644 --- a/src/contraction_sequences.jl +++ b/src/contraction_sequences.jl @@ -2,6 +2,7 @@ using Graphs: vertices using ITensors: ITensor, contract using ITensors.ContractionSequenceOptimization: deepmap, optimal_contraction_sequence using ITensors.NDTensors: Algorithm, @Algorithm_str +using NamedGraphs: parent_vertex_to_vertex using NamedGraphs.Keys: Key function contraction_sequence(tn::Vector{ITensor}; alg="optimal", kwargs...) @@ -10,8 +11,8 @@ end function contraction_sequence(tn::AbstractITensorNetwork; kwargs...) seq_linear_index = contraction_sequence(Vector{ITensor}(tn); kwargs...) - # TODO: Use Functors.fmap? - return deepmap(n -> Key(vertices(tn)[n]), seq_linear_index) + # TODO: Use `Functors.fmap` or `StructWalk`? + return deepmap(n -> Key(parent_vertex_to_vertex(tn, n)), seq_linear_index) end function contraction_sequence(::Algorithm"optimal", tn::Vector{ITensor}) diff --git a/src/formnetworks/abstractformnetwork.jl b/src/formnetworks/abstractformnetwork.jl index 17f647eb..1f235746 100644 --- a/src/formnetworks/abstractformnetwork.jl +++ b/src/formnetworks/abstractformnetwork.jl @@ -1,4 +1,5 @@ using Graphs: induced_subgraph + default_bra_vertex_suffix() = "bra" default_ket_vertex_suffix() = "ket" default_operator_vertex_suffix() = "operator" diff --git a/src/gauging.jl b/src/gauging.jl index f76c91ac..c0a21864 100644 --- a/src/gauging.jl +++ b/src/gauging.jl @@ -4,7 +4,7 @@ using NamedGraphs.PartitionedGraphs: partitionedge function default_bond_tensors(ψ::ITensorNetwork) return DataGraph( - underlying_graph(ψ); edge_data_eltype=Nothing, vertex_data_eltype=ITensor + underlying_graph(ψ); vertex_data_eltype=Nothing, edge_data_eltype=ITensor ) end diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index e257c953..bb94eb84 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -62,13 +62,18 @@ function ttn_svd( # traverse tree outwards from root vertex vs = _default_vertex_ordering(sites, root_vertex) - # ToDo: Add check in ttn_svd that the ordering matches that of find_index_in_tree, which is used in sorteachterm #fermion-sign! - es = _default_edge_ordering(sites, root_vertex) # store edges in fixed ordering relative to root + # TODO: Add check in ttn_svd that the ordering matches that of find_index_in_tree, which is used in sorteachterm #fermion-sign! + # store edges in fixed ordering relative to root + es = _default_edge_ordering(sites, root_vertex) # some things to keep track of - degrees = Dict(v => degree(sites, v) for v in vs) # rank of every TTN tensor in network - Vs = Dict(e => Dict{QN,Matrix{coefficient_type}}() for e in es) # link isometries for SVD compression of TTN - inmaps = Dict{Pair{edgetype_sites,QN},Dict{Vector{Op},Int}}() # map from term in Hamiltonian to incoming channel index for every edge - outmaps = Dict{Pair{edgetype_sites,QN},Dict{Vector{Op},Int}}() # map from term in Hamiltonian to outgoing channel index for every edge + # rank of every TTN tensor in network + degrees = Dict(v => degree(sites, v) for v in vs) + # link isometries for SVD compression of TTN + Vs = Dict(e => Dict{QN,Matrix{coefficient_type}}() for e in es) + # map from term in Hamiltonian to incoming channel index for every edge + inmaps = Dict{Pair{edgetype_sites,QN},Dict{Vector{Op},Int}}() + # map from term in Hamiltonian to outgoing channel index for every edge + outmaps = Dict{Pair{edgetype_sites,QN},Dict{Vector{Op},Int}}() op_cache = Dict{Pair{String,vertextype_sites},ITensor}() @@ -96,14 +101,16 @@ function ttn_svd( is_internal[v] = isempty(sites[v]) if isempty(sites[v]) # FIXME: This logic only works for trivial flux, breaks for nonzero flux - # ToDo: add assert or fix and add test! + # TODO: add assert or fix and add test! sites[v] = [Index(Hflux => 1)] end end + # bond coefficients for incoming edge channels inbond_coefs = Dict( e => Dict{QN,Vector{ITensorMPS.MatElem{coefficient_type}}}() for e in es - ) # bond coefficients for incoming edge channels - site_coef_done = Prod{Op}[] # list of terms for which the coefficient has been added to a site factor + ) + # list of terms for which the coefficient has been added to a site factor + site_coef_done = Prod{Op}[] # temporary symbolic representation of TTN Hamiltonian tempTTN = Dict(v => QNArrElem{Scaled{coefficient_type,Prod{Op}},degrees[v]}[] for v in vs) @@ -203,7 +210,8 @@ function ttn_svd( coutmap = get!( outmaps, edges[dout] => outgoing_qns[edges[dout]], Dict{Vector{Op},Int}() ) - T_inds[dout] = ITensorMPS.posInLink!(coutmap, outgoing[edges[dout]]) # add outgoing channel + # add outgoing channel + T_inds[dout] = ITensorMPS.posInLink!(coutmap, outgoing[edges[dout]]) T_qns[dout] = outgoing_qns[edges[dout]] end # if term starts at this site, add its coefficient as a site factor @@ -211,7 +219,8 @@ function ttn_svd( if (isnothing(dim_in) || T_inds[dim_in] == -1) && ITensors.argument(term) ∉ site_coef_done site_coef = ITensors.coefficient(term) - site_coef = convert(coefficient_type, site_coef) # required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real + # required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real + site_coef = convert(coefficient_type, site_coef) push!(site_coef_done, ITensors.argument(term)) end # add onsite identity for interactions passing through vertex @@ -268,7 +277,7 @@ function ttn_svd( for v in vs # redo the whole thing like before - # ToDo: use neighborhood instead of going through all edges, see above + # TODO: use neighborhood instead of going through all edges, see above edges = align_and_reorder_edges(incident_edges(sites, v), es) dim_in = findfirst(e -> dst(e) == v, edges) dims_out = findall(e -> src(e) == v, edges) @@ -337,7 +346,8 @@ function ttn_svd( for ((b, q_op), m) in blocks Op = computeSiteProd(sites, Prod(q_op)) - if hasqns(Op) # FIXME: this may not be safe, we may want to check for the equivalent (zero tensor?) case in the dense case as well + if hasqns(Op) + # FIXME: this may not be safe, we may want to check for the equivalent (zero tensor?) case in the dense case as well iszero(nnzblocks(Op)) && continue end sq = flux(Op) @@ -363,7 +373,7 @@ function ttn_svd( if is_internal[v] H[v] += iT else - #ToDo: Remove this assert since it seems to be costly + #TODO: Remove this assert since it seems to be costly #if hasqns(iT) # @assert flux(iT * Op) == Hflux #end @@ -375,18 +385,22 @@ function ttn_svd( # add starting and ending identity operators idT = zeros(coefficient_type, linkdims...) if isnothing(dim_in) - idT[ones(Int, degrees[v])...] = 1.0 # only one real starting identity + # only one real starting identity + idT[ones(Int, degrees[v])...] = 1.0 end # ending identities are a little more involved if !isnothing(dim_in) - idT[linkdims...] = 1.0 # place identity if all channels end + # place identity if all channels end + idT[linkdims...] = 1.0 # place identity from start of incoming channel to start of each single outgoing channel, and end all other channels idT_end_inds = [linkdims...] - idT_end_inds[dim_in] = 1 #this should really be an int + #this should really be an int + idT_end_inds[dim_in] = 1 for dout in dims_out idT_end_inds[dout] = 1 idT[idT_end_inds...] = 1.0 - idT_end_inds[dout] = linkdims[dout] # reset + # reset + idT_end_inds[dout] = linkdims[dout] end end T = itensor(idT, _linkinds) diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl index f035790f..ba150fa6 100644 --- a/test/test_opsum_to_ttn.jl +++ b/test/test_opsum_to_ttn.jl @@ -35,7 +35,6 @@ end @testset "OpSum to TTN" begin # small comb tree auto_fermion_enabled = ITensors.using_auto_fermion() - ITensors.disable_auto_fermion() # ToDo: remove when autofermion incompatibility with no QNs is fixed tooth_lengths = fill(2, 3) c = named_comb_tree(tooth_lengths) @@ -43,7 +42,7 @@ end # linearized version linear_order = [4, 1, 2, 5, 3, 6] - vmap = Dictionary(vertices(is)[linear_order], 1:length(linear_order)) + vmap = Dictionary(collect(vertices(is))[linear_order], eachindex(linear_order)) sites = only.(collect(vertex_data(is)))[linear_order] # test with next-to-nearest-neighbor Ising Hamiltonian @@ -63,7 +62,7 @@ end @testset "Svd approach" for root_vertex in leaf_vertices(is) # get TTN Hamiltonian directly - Hsvd = ttn(H, is; root_vertex=root_vertex, cutoff=1e-10) + Hsvd = ttn(H, is; root_vertex, cutoff=1e-10) # get corresponding MPO Hamiltonian Hline = ITensorMPS.MPO(relabel_sites(H, vmap), sites) # compare resulting dense Hamiltonians @@ -73,8 +72,7 @@ end end @test Tttno ≈ Tmpo rtol = 1e-6 - # this breaks for longer range interactions - Hsvd_lr = ttn(Hlr, is; root_vertex=root_vertex, algorithm="svd", cutoff=1e-10) + Hsvd_lr = ttn(Hlr, is; root_vertex, cutoff=1e-10) Hline_lr = ITensorMPS.MPO(relabel_sites(Hlr, vmap), sites) @disable_warn_order begin Tttno_lr = prod(Hline_lr) @@ -89,7 +87,9 @@ end @testset "Multiple onsite terms (regression test for issue #62)" begin auto_fermion_enabled = ITensors.using_auto_fermion() - ITensors.disable_auto_fermion() # ToDo: remove when autofermion incompatibility with no QNs is fixed + if !auto_fermion_enabled + ITensors.enable_auto_fermion() + end grid_dims = (2, 1) g = named_grid(grid_dims) s = siteinds("S=1/2", g) @@ -120,7 +120,7 @@ end # linearized version linear_order = [4, 1, 2, 5, 3, 6] - vmap = Dictionary(vertices(is)[linear_order], 1:length(linear_order)) + vmap = Dictionary(collect(vertices(is))[linear_order], eachindex(linear_order)) sites = only.(collect(vertex_data(is)))[linear_order] # test with next-to-nearest-neighbor Ising Hamiltonian @@ -140,7 +140,7 @@ end @testset "Svd approach" for root_vertex in leaf_vertices(is) # get TTN Hamiltonian directly - Hsvd = ttn(H, is; root_vertex=root_vertex, cutoff=1e-10) + Hsvd = ttn(H, is; root_vertex, cutoff=1e-10) # get corresponding MPO Hamiltonian Hline = ITensorMPS.MPO(relabel_sites(H, vmap), sites) # compare resulting sparse Hamiltonians @@ -151,8 +151,7 @@ end end @test Tttno ≈ Tmpo rtol = 1e-6 - # this breaks for longer range interactions ###not anymore - Hsvd_lr = ttn(Hlr, is; root_vertex=root_vertex, algorithm="svd", cutoff=1e-10) + Hsvd_lr = ttn(Hlr, is; root_vertex, cutoff=1e-10) Hline_lr = ITensorMPS.MPO(relabel_sites(Hlr, vmap), sites) @disable_warn_order begin Tttno_lr = prod(Hline_lr) @@ -184,7 +183,7 @@ end @testset "Svd approach" for root_vertex in leaf_vertices(is) # get TTN Hamiltonian directly - Hsvd = ttn(H, is; root_vertex=root_vertex, cutoff=1e-10) + Hsvd = ttn(H, is; root_vertex, cutoff=1e-10) # get corresponding MPO Hamiltonian sites = [only(is[v]) for v in reverse(post_order_dfs_vertices(c, root_vertex))] vmap = Dictionary(reverse(post_order_dfs_vertices(c, root_vertex)), 1:length(sites)) @@ -201,7 +200,8 @@ end @test norm(Tttno) > 0 @test norm(Tmpo) ≈ norm(Tttno) rtol = 1e-6 - @test_broken Tmpo ≈ Tttno # ToDo fix comparison for fermionic tensors + # TODO: fix comparison for fermionic tensors + @test_broken Tmpo ≈ Tttno # In the meantime: matricize tensors and convert to dense Matrix to compare element by element dTmm = to_matrix(Tmpo) dTtm = to_matrix(Tttno) @@ -236,14 +236,14 @@ end # linearized version linear_order = [4, 1, 2, 5, 3, 6] - vmap = Dictionary(getindices(vertices(is), linear_order), eachindex(linear_order)) + vmap = Dictionary(collect(vertices(is))[linear_order], eachindex(linear_order)) sites = only.(filter(d -> !isempty(d), collect(vertex_data(is_missing_site))))[linear_order] J1 = -1 J2 = 2 h = 0.5 # connectivity of the Hamiltonian is that of the original comb graph - H = ModelHamiltonians.heisenberg(c; J1=J1, J2=J2, h=h) + H = ModelHamiltonians.heisenberg(c; J1, J2, h) # add combination of longer range interactions Hlr = copy(H) @@ -254,7 +254,7 @@ end @testset "Svd approach" for root_vertex in leaf_vertices(is) # get TTN Hamiltonian directly - Hsvd = ttn(H, is_missing_site; root_vertex=root_vertex, cutoff=1e-10) + Hsvd = ttn(H, is_missing_site; root_vertex, cutoff=1e-10) # get corresponding MPO Hamiltonian Hline = ITensorMPS.MPO(relabel_sites(H, vmap), sites) @@ -265,9 +265,7 @@ end end @test Tttno ≈ Tmpo rtol = 1e-6 - Hsvd_lr = ttn( - Hlr, is_missing_site; root_vertex=root_vertex, algorithm="svd", cutoff=1e-10 - ) + Hsvd_lr = ttn(Hlr, is_missing_site; root_vertex, cutoff=1e-10) Hline_lr = ITensorMPS.MPO(relabel_sites(Hlr, vmap), sites) @disable_warn_order begin Tttno_lr = prod(Hline_lr)