Fix more tests

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

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...)

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})

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"

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
 

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,15 +210,17 @@ 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
       site_coef = one(coefficient_type)
       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)

test/test_opsum_to_ttn.jl

@@ -35,15 +35,14 @@ 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)
 
     is = siteinds("S=1/2", c)
 
     # 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)