From 5ab9951eea7a2a37e04df8784b0f89ca1a091a0d Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Mon, 8 Jan 2024 12:51:00 -0500
Subject: [PATCH 01/32] Add blocksparse Hamiltonian constructor for trees from
 OpSum, qn_svdTTN.

---
 src/treetensornetworks/opsum_to_ttn.jl | 396 ++++++++++++++++++++++++-
 test/test_opsum_to_ttn.jl              |  65 ++++
 2 files changed, 446 insertions(+), 15 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 33fc0af2..58ef1ab0 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -71,7 +71,6 @@ function svdTTN(
     edge_in = (isnothing(dim_in) ? [] : edges[dim_in])
     dims_out = findall(e -> src(e) == v, edges)
     edges_out = edges[dims_out]
-
     # sanity check, leaves only have single incoming or outgoing edge
     @assert !isempty(dims_out) || !isnothing(dim_in)
     (isempty(dims_out) || isnothing(dim_in)) && @assert is_leaf(sites, v)
@@ -186,9 +185,10 @@ function svdTTN(
           z = ct
           temp_inds = copy(T_inds)
           for (i, d) in enumerate(normal_dims)
+            #@show size(Vv[d]), T_inds[d],c[i]
             V_factor = Vv[d][T_inds[d], c[i]]
             z *= (d == dim_in ? conj(V_factor) : V_factor) # conjugate incoming isemetry factor
-            temp_inds[d] = 1 + c[i]
+            temp_inds[d] = 1 + c[i]     ##the offset by one comes from the 1,...,1 padding (starting/ending...)
           end
           T[temp_inds...] += z
         end
@@ -221,6 +221,349 @@ function svdTTN(
   return H
 end
 
+
+function qn_svdTTN(
+  os::OpSum{C}, sites::IndsNetwork{VT,<:Index}, root_vertex::VT; kwargs...
+)::TTN where {C,VT}
+  mindim::Int = get(kwargs, :mindim, 1)
+  maxdim::Int = get(kwargs, :maxdim, 10000)
+  cutoff::Float64 = get(kwargs, :cutoff, 1e-15)
+
+  ValType = ITensors.determineValType(ITensors.terms(os))
+
+  # traverse tree outwards from root vertex
+  vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
+  es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
+  # some things to keep track of
+  ranks = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
+  Vs = Dict(e => Dict{QN,Matrix{ValType}}() for e in es)                                  # link isometries for SVD compression of TTN
+  #Vs = Dict(e => Matrix{ValType}(undef, 1, 1) for e in es)    ##no qn implementation                              # link isometries for SVD compression of TTN
+  #inmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                    # map from term in Hamiltonian to incoming channel index for every edge
+  #outmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   # map from term in Hamiltonian to outgoing channel index for every edge
+  inmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()
+  outmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()
+  
+  ##ToDo implement op cache? Not used?
+  op_cache = Dict{Pair{String,VT},ITensor}()   ###FIXME needs to be Pair{String,Vertex}
+  #@show VT
+  #@show op_cache
+  ##copied calcQN from ITensors MPO implementation  
+  function calcQN(term::Vector{Op})
+    q = QN()
+    for st in term
+      op_tensor = get(op_cache, ITensors.which_op(st) => ITensors.site(st), nothing)
+      if op_tensor === nothing
+        op_tensor = op(sites[ITensors.site(st)], ITensors.which_op(st); ITensors.params(st)...)
+        op_cache[ITensors.which_op(st) => ITensors.site(st)] = op_tensor
+      end
+      q -= flux(op_tensor)
+    end
+    return q
+  end
+  Hflux = -calcQN(terms(first(terms(os))))
+  ###FIXME: add assert that Hflux is the zero-element for addition
+  ##CHECK: PERFORMANCE/FOOTPRINT? inbond_coeffs equivalent in MPO implementation is inside the loop
+  ##we may want to generate only for the adjacent edges inside loop over vertices instead?
+  inbond_coefs = Dict(e => Dict{QN,Vector{ITensors.MatElem{ValType}}}() 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
+
+  # temporary symbolic representation of TTN Hamiltonian
+  tempTTN = Dict(v => QNArrElem{Scaled{C,Prod{Op}},ranks[v]}[] for v in vs)
+
+  # build compressed finite state machine representation
+  for v in vs
+    # for every vertex, find all edges that contain this vertex
+    edges = filter(e -> dst(e) == v || src(e) == v, es)
+    # use the corresponding ordering as index order for tensor elements at this site
+    dim_in = findfirst(e -> dst(e) == v, edges)
+    edge_in = (isnothing(dim_in) ? [] : edges[dim_in])
+    dims_out = findall(e -> src(e) == v, edges)
+    edges_out = edges[dims_out]
+
+
+    # sanity check, leaves only have single incoming or outgoing edge
+    @assert !isempty(dims_out) || !isnothing(dim_in)
+    (isempty(dims_out) || isnothing(dim_in)) && @assert is_leaf(sites, v)
+
+    for term in os
+      # loop over OpSum and pick out terms that act on current vertex
+      crosses_vertex(term, sites, v) || continue
+
+      # filter out factors that come in from the direction of the incoming edge
+      incoming = filter(
+        t -> edge_in ∈ edge_path(sites, ITensors.site(t), v), ITensors.terms(term)
+      )
+      # also store all non-incoming factors in standard order, used for channel merging
+      not_incoming = filter(
+        t -> edge_in ∉ edge_path(sites, ITensors.site(t), v), ITensors.terms(term)
+      )
+      # filter out factor that acts on current vertex
+      onsite = filter(t -> (ITensors.site(t) == v), ITensors.terms(term))
+      # for every outgoing edge, filter out factors that go out along that edge
+      outgoing = Dict(
+        e => filter(t -> e ∈ edge_path(sites, v, ITensors.site(t)), ITensors.terms(term))
+        for e in edges_out
+      )
+
+      ##compute QNs
+      incoming_qn=calcQN(incoming)
+      not_incoming_qn=calcQN(not_incoming)  ###does this work?
+      outgoing_qns=Dict(e => calcQN(outgoing[e]) for e in edges_out)
+
+      
+      ##sanity check that qn of not_incoming is equal to sum of outgoing_qns?
+      site_qn=calcQN(onsite)
+      #@show incoming_qn, not_incoming_qn, values(outgoing_qns), Hflux
+      #@show Hflux , not_incoming_qn + incoming_qn +Hflux
+      #@assert sum([values(outgoing_qns)...,Hflux])  == -site_qn - incoming_qn + Hflux
+      
+      #@show sum(values(outgoing_qns)),site_qn, incoming_qn
+      # translate into tensor entry #(within a QNBlock)
+      T_inds = MVector{ranks[v]}(fill(-1, ranks[v]))
+      T_qns = MVector{ranks[v]}(fill(QN(), ranks[v]))
+      
+      bond_row = -1
+      bond_col = -1
+      if !isempty(incoming)
+        coutmap=get!(outmaps,edge_in=>-not_incoming_qn,Dict{Vector{Op},Int}())
+        cinmap=get!(inmaps,edge_in=>incoming_qn,Dict{Vector{Op},Int}())
+        
+        bond_row = ITensors.posInLink!(cinmap, incoming)
+        bond_col = ITensors.posInLink!(coutmap, not_incoming) # get incoming channel
+        bond_coef = convert(ValType, ITensors.coefficient(term))
+        q_inbond_coefs = get!(inbond_coefs[edge_in],incoming_qn, ITensors.MatElem{ValType}[])
+        push!(q_inbond_coefs, ITensors.MatElem(bond_row, bond_col, bond_coef))
+        T_inds[dim_in] = bond_col
+        T_qns[dim_in] = incoming_qn
+      end
+      for dout in dims_out
+        #@show outgoing_qns[edges[dout]]
+        #@show keys(outmaps)
+        coutmap=get!(outmaps,edges[dout]=>-outgoing_qns[edges[dout]],Dict{Vector{Op},Int}())
+        T_inds[dout] = ITensors.posInLink!(coutmap, outgoing[edges[dout]]) # add outgoing channel
+        T_qns[dout] = -outgoing_qns[edges[dout]]  ###minus sign to account for outgoing arrow direction
+      end
+      # if term starts at this site, add its coefficient as a site factor
+      site_coef = one(C)
+      if (isnothing(dim_in) || T_inds[dim_in] == -1) &&
+        ITensors.argument(term) ∉ site_coef_done
+        site_coef = ITensors.coefficient(term)
+        push!(site_coef_done, ITensors.argument(term))
+      end
+      # add onsite identity for interactions passing through vertex
+      if isempty(onsite)
+        if !ITensors.using_auto_fermion() && isfermionic(incoming, sites)
+          error("No verified fermion support for automatic TTN constructor!")
+        else
+          push!(onsite, Op("Id", v))
+        end
+      end
+      # save indices and value of symbolic tensor entry
+      el = QNArrElem(T_qns, T_inds, site_coef * Prod(onsite))
+      push!(tempTTN[v], el)
+    end
+    ITensors.remove_dups!(tempTTN[v])
+    # manual truncation: isometry on incoming edge
+    if !isnothing(dim_in) && !isempty(inbond_coefs[edges[dim_in]])
+      for (q, mat) in inbond_coefs[edges[dim_in]]
+        M = ITensors.toMatrix(mat)
+        U, S, V = svd(M)
+        P = S .^ 2
+        truncate!(P; maxdim=maxdim, cutoff=cutoff, mindim=mindim)
+        tdim = length(P)
+        nc = size(M, 2) ###CHECK: verify whether this is still correct, diverges from opsum_to_mpo_qn
+        Vs[edges[dim_in]][q] = Matrix{ValType}(V[1:nc, 1:tdim]) ###CHECK: verify whether this is still correct, diverges from opsum_to_mpo_qn
+      end
+    end
+  end
+
+  # compress this tempTTN representation into dense form
+
+  #link_space = dictionary([
+  #  e => Index((isempty(outmaps[e]) ? 0 : size(Vs[e], 2)) + 2, edge_tag(e)) for e in es
+  #])
+  ###FIXME: that's an ugly constructor... ###QNIndex(undef) not defined
+  link_space = Dict{NamedGraphs.NamedEdge{VT},Index}()
+  linkdir = ITensors.using_auto_fermion() ? ITensors.In : ITensors.Out ###FIXME: ??
+  for e in es ###this one might also be tricky, in case there's no Vq defined on the edge
+    ###also tricky w.r.t. the link direction
+    ### ---> should we just infer whether the edge is incoming or outgoing and choose linkdir accordingly?
+    ### ===> this is basically decided by daggering in the MPO code (we should be able to stick to this pattern)
+    qi = Vector{Pair{QN,Int}}()
+    push!(qi, QN() => 1)
+    for (q, Vq) in Vs[e]
+      cols = size(Vq, 2)
+      if ITensors.using_auto_fermion() # <fermions>
+        push!(qi, (-q) => cols) ###the minus sign comes from the opposite linkdir ?
+      else
+        push!(qi, q => cols)
+      end
+    end
+    push!(qi, Hflux => 1)
+    link_space[e] = Index(qi...; tags=edge_tag(e), dir=linkdir) ##FIXME: maybe a single linkdir is not the right thing to do here
+  end
+
+  H = TTN(sites)
+  ##Copy from ITensors qn_svdMPO
+  # Find location where block of Index i
+  # matches QN q, but *not* 1 or dim(i)
+  # which are special ending/starting states
+  function qnblock(i::Index, q::QN)
+    for b in 2:(nblocks(i) - 1)
+      flux(i, Block(b)) == q && return b
+    end
+    return error("Could not find block of QNIndex with matching QN")
+  end
+  qnblockdim(i::Index, q::QN) = blockdim(i, qnblock(i, q))
+  
+  for v in vs
+    #@show v
+    # redo the whole thing like before
+    edges = filter(e -> dst(e) == v || src(e) == v, es)
+    dim_in = findfirst(e -> dst(e) == v, edges)
+    dims_out = findall(e -> src(e) == v, edges)
+    
+    #edge_in = (isnothing(dim_in) ? [] : edges[dim_in])
+    #edges_out = edges[dims_out]
+    # slice isometries at this vertex
+    Vv = [Vs[e] for e in edges]
+
+    linkinds = [link_space[e] for e in edges]
+    
+
+    #####
+    # until here things are pretty straightforward at first sight
+    #####
+    #begin_block = Dict{Tuple{QN,Vector{Op}},Array{ValType,ranks[v]}}()
+    #cont_block = Dict{Tuple{QN,Vector{Op}},Array{ValType,ranks[v]}}()
+    #end_block = Dict{Tuple{QN,Vector{Op}},Array{ValType,ranks[v]}}()
+    #onsite_block = Dict{Tuple{QN,Vector{Op}},Array{ValType,ranks[v]}}()
+    #general_block = Dict{Tuple{QN,Vector{Op}},Array{ValType,ranks[v]}}()
+    blocks = Dict{Tuple{Block{ranks[v]},Vector{Op}},Array{ValType,ranks[v]}}()
+    ##construct blocks
+    linkdims = dim.(linkinds)
+    for el in tempTTN[v]
+      t = el.val
+      (abs(coefficient(t)) > eps()) || continue
+      block_helper_inds=fill(-1,ranks[v])
+      T_inds=el.idxs
+      #@show T_inds
+      T_qns=el.qn_idxs
+      ct = convert(ValType, coefficient(t))
+      sublinkdims= [(T_inds[i] == -1 ? 1 : qnblockdim(linkinds[i], T_qns[i])) for i in 1:ranks[v]]
+      zero_arr() = zeros(ValType, sublinkdims...)
+      terminal_dims = findall(d -> T_inds[d] == -1, 1:ranks[v])   # directions in which term starts or ends
+      normal_dims = findall(d -> T_inds[d] ≠ -1, 1:ranks[v])      # normal dimensions, do truncation thingies
+      T_inds[terminal_dims] .= 1                                  # start in channel 1  ###??
+      block_helper_inds[terminal_dims].=1
+      for dout in filter(d -> d ∈ terminal_dims, dims_out)
+        T_inds[dout] = sublinkdims[dout]
+        @assert T_inds[dout]==1                         # end in channel linkdims[d] for each dimension d
+        block_helper_inds[dout]=nblocks(linkinds[dout])
+      end
+      #non-trivial helper inds
+      for d in normal_dims
+        #if d in dims_out
+        block_helper_inds[d]=qnblock(linkinds[d], T_qns[d])
+      end
+
+      @assert all(block_helper_inds .!= -1)
+      theblock=Block(Tuple(block_helper_inds))
+      
+      ### T_qns takes into account incoming/outgoing arrow direction, while Vs are stored for incoming arrow
+      ### flip sign of QN for outgoing arrows, to match incoming arrow definition of Vs
+      iT_qns=deepcopy(T_qns)
+      for d in dims_out
+        iT_qns[d]=-iT_qns[d]
+      end
+
+      if isempty(normal_dims)
+
+        M = get!(blocks, (theblock, terms(t)),zero_arr())
+        @assert product(size(M))==1 ###FIXME: figure out how to do this assert properly
+        M[] += ct
+      else
+        M = get!(blocks, (theblock, terms(t)),zero_arr())
+        dim_ranges = Tuple(size(Vv[d][iT_qns[d]],2) for d in normal_dims)
+        for c in CartesianIndices(dim_ranges)
+          z = ct
+          temp_inds = copy(T_inds)
+           for (i, d) in enumerate(normal_dims)            
+            V_factor = Vv[d][iT_qns[d]][T_inds[d], c[i]]
+            z *= (d == dim_in ? conj(V_factor) : V_factor) # conjugate incoming isometry factor
+            temp_inds[d] = c[i]
+          end
+          M[temp_inds...] += z 
+        end
+      end
+
+    end
+    
+    
+    H[v] = ITensor()
+
+    _linkinds=copy(linkinds)
+    if !isnothing(dim_in)
+      _linkinds[dim_in]=dag(_linkinds[dim_in])
+    end
+
+    for ((b,q_op),m) in blocks
+      #@show q_op
+      #op_prod = q_op[2]
+      ###FIXME: make this work if there's no physical degree of freedom at site
+      Op = computeSiteProd(sites, Prod(q_op))  ###FIXME: is this implemented?
+      (nnzblocks(Op) == 0) && continue  ###FIXME: this one may be breaking for no physical indices on site
+      sq = flux(Op)
+      if !isnothing(sq)
+        if ITensors.using_auto_fermion()
+          @assert false
+          ###assuming we want index ordering dim_in, phys..., dims_out...
+          ###but are assembling the tensor by tensor product (dim_in,dims_out...) x (phys...)
+          ###we need to undo the autofermion logic.
+          ###this requires the perm and the quantum numbers
+          ###perm is 1,3,2 (switching phys and dims_out dimensions)
+          ###quantum numbers need to be retrieved from block and linkinds
+          ###basically first(space(linkinds[i])[j]) for (i,j) in block
+          ###this logic requires dim_in to be 1 always, i don't think this is guaranteed
+          ###however, maybe it is required when using autofermion?
+          ###so we would have to figure out a perm to bring it into that format
+          
+          perm=(1,3,2)
+
+        end
+      end
+      T = ITensors.BlockSparseTensor(ValType, [b], _linkinds)  ###why the brackets around b? does the constructor dispatch on a vector of blocks?
+      T[b] .= m
+      H[v] += (itensor(T)*Op)
+      
+    end
+    
+    linkdims=dim.(linkinds)
+    # add starting and ending identity operators
+    idT = zeros(ValType, linkdims...)
+    if isnothing(dim_in)
+      idT[ones(Int, ranks[v])...] = 1.0 # only one real starting identity
+    end
+    # ending identities are a little more involved
+    if !isnothing(dim_in)
+      idT[linkdims...] = 1.0 # place identity if all channels end
+      # 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
+      for dout in dims_out
+        idT_end_inds[dout] = 1
+        idT[idT_end_inds...] = 1.0
+        idT_end_inds[dout] = linkdims[dout] # reset
+      end
+    end
+    T = itensor(idT, _linkinds)  ###definitely need to dagger link_in , or do that beforehand?
+    H[v] += T * ITensorNetworks.computeSiteProd(sites, Prod([(Op("Id", v))]))
+  end
+
+  return H
+end
+
+
 # 
 # Tree adaptations of functionalities in ITensors.jl/src/physics/autompo/opsum_to_mpo_generic.jl
 # 
@@ -251,6 +594,7 @@ function computeSiteProd(sites::IndsNetwork{V,<:Index}, ops::Prod{Op})::ITensor
   return T
 end
 
+###CHECK/ToDo: Understand the fermion logic on trees...
 # changed `isless_site` to use tree vertex ordering relative to root
 function sorteachterm(os::OpSum, sites::IndsNetwork{V,<:Index}, root_vertex::V) where {V}
   os = copy(os)
@@ -341,20 +685,13 @@ function TTN(
   os = ITensors.sortmergeterms(os) # not exported
 
   if hasqns(first(first(vertex_data(sites))))
-    if !is_path_graph(sites)
-      error(
-        "OpSum → TTN constructor for QN conserving tensor networks only works for path/linear graphs right now.",
-      )
+    ###added feature
+    T=qn_svdTTN(os,sites,root_vertex; kwargs...)
+    if splitblocks
+      error("splitblocks not yet implemented for AbstractTreeTensorNetwork.")
+      T = ITensors.splitblocks(linkinds, T) # TODO: make this work
     end
-    # Use `ITensors.MPO` for now until general TTN constructor
-    # works for QNs.
-    # TODO: Check it is a path graph and get a linear arrangement!
-    sites_linear_vertices = [only(sites[v]) for v in vertices(sites)]
-    vertices_to_linear_vertices = Dictionary(vertices(sites), eachindex(vertices(sites)))
-    os_linear_vertices = replace_vertices(os, vertices_to_linear_vertices)
-    mpo = MPO(os_linear_vertices, sites_linear_vertices)
-    tn = TTN(Dictionary(vertices(sites), [mpo[v] for v in 1:nv(sites)]))
-    return tn
+    return T
   end
   T = svdTTN(os, sites, root_vertex; kwargs...)
   if splitblocks
@@ -424,6 +761,35 @@ function Base.isless(a1::ArrElem{T,N}, a2::ArrElem{T,N})::Bool where {T,N}
   return a1.val < a2.val
 end
 
+#################################
+# QNArrElem (sparse array with QNs) #
+#################################
+
+struct QNArrElem{T,N}
+  qn_idxs::MVector{N,QN}
+  idxs::MVector{N,Int}
+  val::T
+end
+
+function Base.:(==)(a1::QNArrElem{T,N}, a2::QNArrElem{T,N})::Bool where {T,N}
+  return (a1.idxs == a2.idxs && a1.val == a2.val && a1.qn_idxs == a2.qn_idxs)
+end
+
+function Base.isless(a1::QNArrElem{T,N}, a2::QNArrElem{T,N})::Bool where {T,N}
+  ###two separate loops s.t. the MPS case reduces to the ITensors Implementation of QNMatElem
+  for n in 1:N
+    if a1.qn_idxs[n] != a2.qn_idxs[n]
+      return a1.qn_idxs[n] < a2.qn_idxs[n]
+    end
+  end
+  for n in 1:N
+    if a1.idxs[n] != a2.idxs[n]
+      return a1.idxs[n] < a2.idxs[n]
+    end
+  end
+  return a1.val < a2.val
+end
+
 # 
 # Sparse finite state machine construction
 # 
diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index 267eeaa5..096c4e0c 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -4,11 +4,15 @@ using ITensorNetworks
 using Random
 using Test
 
+
+
 @testset "OpSum to TTN converter" begin
   @testset "OpSum to TTN" begin
     # small comb tree
     tooth_lengths = fill(2, 3)
     c = named_comb_tree(tooth_lengths)
+    @show c
+    
     is = siteinds("S=1/2", c)
 
     # linearized version
@@ -96,4 +100,65 @@ using Test
 
     @test H1 + H2 ≈ H3 rtol = 1e-6
   end
+  
+  @testset "OpSum to TTN QN" begin
+    # small comb tree
+    tooth_lengths = fill(2, 3)
+    c = named_comb_tree(tooth_lengths)
+    is = siteinds("S=1/2", c;conserve_qns=true)
+    is_noqns=copy(is)
+    for v in vertices(is)
+      is_noqns[v]=removeqns(is_noqns[v])
+    end
+
+    # linearized version
+    linear_order = [4, 1, 2, 5, 3, 6]
+    vmap = Dictionary(vertices(is)[linear_order], 1:length(linear_order))
+    sites = only.(collect(vertex_data(is)))[linear_order]
+
+    # test with next-to-nearest-neighbor Ising Hamiltonian
+    J1 = -1
+    J2 = 2
+    h = 0.5
+    H = heisenberg(c; J1=J1, J2=J2)
+
+    # add combination of longer range interactions
+    Hlr = copy(H)
+    Hlr += 5, "Z", (1, 2), "Z", (2, 2)
+    Hlr += -4, "Z", (1, 1), "Z", (2, 2)
+    Hlr += 2.0, "Z", (2, 2), "Z", (3, 2)
+    Hlr += -1.0, "Z", (1, 2), "Z", (3, 1)
+
+    # root_vertex = (1, 2)
+    # println(leaf_vertices(is))
+
+    @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)
+      # get corresponding MPO Hamiltonian
+      Hline = MPO(relabel_sites(H, vmap), sites)
+      # compare resulting sparse Hamiltonians
+      
+      
+      @disable_warn_order begin
+        Tmpo = prod(Hline)
+        Tttno = contract(Hsvd)
+        
+      end
+      @test Tttno ≈ Tmpo rtol = 1e-6
+      
+      # this breaks for longer range interactions ###not anymore
+      Hsvd_lr = TTN(Hlr, is; root_vertex=root_vertex, method=:svd, cutoff=1e-10)
+      Hline_lr = MPO(relabel_sites(Hlr, vmap), sites)
+      @disable_warn_order begin
+        Tttno_lr = prod(Hline_lr)
+        Tmpo_lr = contract(Hsvd_lr)
+      end
+      @test Tttno_lr ≈ Tmpo_lr rtol = 1e-6
+      
+    end
+    
+    
+  end
+  
 end

From e16f4401c47097aa2c2d1f33ba0cdcd174928d62 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Mon, 8 Jan 2024 15:31:29 -0500
Subject: [PATCH 02/32] Fix bug in assert for onsite terms. Clean up comments.

---
 src/treetensornetworks/opsum_to_ttn.jl | 99 +++++++++-----------------
 1 file changed, 34 insertions(+), 65 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 58ef1ab0..2e542a2f 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -236,18 +236,13 @@ function qn_svdTTN(
   es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
   # some things to keep track of
   ranks = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
-  Vs = Dict(e => Dict{QN,Matrix{ValType}}() for e in es)                                  # link isometries for SVD compression of TTN
-  #Vs = Dict(e => Matrix{ValType}(undef, 1, 1) for e in es)    ##no qn implementation                              # link isometries for SVD compression of TTN
-  #inmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                    # map from term in Hamiltonian to incoming channel index for every edge
-  #outmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   # map from term in Hamiltonian to outgoing channel index for every edge
-  inmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()
-  outmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()
+  Vs = Dict(e => Dict{QN,Matrix{ValType}}() for e in es)                                    # link isometries for SVD compression of TTN
+  #inmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   
+  #outmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   
+  inmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()                  # map from term in Hamiltonian to incoming channel index for every edge
+  outmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()                 # map from term in Hamiltonian to outgoing channel index for every edge
   
-  ##ToDo implement op cache? Not used?
-  op_cache = Dict{Pair{String,VT},ITensor}()   ###FIXME needs to be Pair{String,Vertex}
-  #@show VT
-  #@show op_cache
-  ##copied calcQN from ITensors MPO implementation  
+  op_cache = Dict{Pair{String,VT},ITensor}()
   function calcQN(term::Vector{Op})
     q = QN()
     for st in term
@@ -260,13 +255,10 @@ function qn_svdTTN(
     end
     return q
   end
+  
   Hflux = -calcQN(terms(first(terms(os))))
-  ###FIXME: add assert that Hflux is the zero-element for addition
-  ##CHECK: PERFORMANCE/FOOTPRINT? inbond_coeffs equivalent in MPO implementation is inside the loop
-  ##we may want to generate only for the adjacent edges inside loop over vertices instead?
   inbond_coefs = Dict(e => Dict{QN,Vector{ITensors.MatElem{ValType}}}() 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
-
   # temporary symbolic representation of TTN Hamiltonian
   tempTTN = Dict(v => QNArrElem{Scaled{C,Prod{Op}},ranks[v]}[] for v in vs)
 
@@ -280,7 +272,6 @@ function qn_svdTTN(
     dims_out = findall(e -> src(e) == v, edges)
     edges_out = edges[dims_out]
 
-
     # sanity check, leaves only have single incoming or outgoing edge
     @assert !isempty(dims_out) || !isnothing(dim_in)
     (isempty(dims_out) || isnothing(dim_in)) && @assert is_leaf(sites, v)
@@ -304,31 +295,25 @@ function qn_svdTTN(
         e => filter(t -> e ∈ edge_path(sites, v, ITensors.site(t)), ITensors.terms(term))
         for e in edges_out
       )
-
-      ##compute QNs
+      
+      # compute QNs
       incoming_qn=calcQN(incoming)
       not_incoming_qn=calcQN(not_incoming)  ###does this work?
       outgoing_qns=Dict(e => calcQN(outgoing[e]) for e in edges_out)
-
-      
-      ##sanity check that qn of not_incoming is equal to sum of outgoing_qns?
       site_qn=calcQN(onsite)
-      #@show incoming_qn, not_incoming_qn, values(outgoing_qns), Hflux
-      #@show Hflux , not_incoming_qn + incoming_qn +Hflux
-      #@assert sum([values(outgoing_qns)...,Hflux])  == -site_qn - incoming_qn + Hflux
       
-      #@show sum(values(outgoing_qns)),site_qn, incoming_qn
-      # translate into tensor entry #(within a QNBlock)
+      # initialize QNArrayElement indices and quantum numbers 
       T_inds = MVector{ranks[v]}(fill(-1, ranks[v]))
       T_qns = MVector{ranks[v]}(fill(QN(), ranks[v]))
-      
+      # initialize ArrayElement indices for inbond_coefs
       bond_row = -1
       bond_col = -1
       if !isempty(incoming)
+        # get the correct map from edge=>QN to term and channel
         coutmap=get!(outmaps,edge_in=>-not_incoming_qn,Dict{Vector{Op},Int}())
         cinmap=get!(inmaps,edge_in=>incoming_qn,Dict{Vector{Op},Int}())
-        
-        bond_row = ITensors.posInLink!(cinmap, incoming)
+        # this checks if term exists on edge=>QN ( otherwise insert it) and returns it's index
+        bond_row = ITensors.posInLink!(cinmap, incoming)  
         bond_col = ITensors.posInLink!(coutmap, not_incoming) # get incoming channel
         bond_coef = convert(ValType, ITensors.coefficient(term))
         q_inbond_coefs = get!(inbond_coefs[edge_in],incoming_qn, ITensors.MatElem{ValType}[])
@@ -417,37 +402,24 @@ function qn_svdTTN(
   qnblockdim(i::Index, q::QN) = blockdim(i, qnblock(i, q))
   
   for v in vs
-    #@show v
     # redo the whole thing like before
     edges = filter(e -> dst(e) == v || src(e) == v, es)
     dim_in = findfirst(e -> dst(e) == v, edges)
     dims_out = findall(e -> src(e) == v, edges)
     
-    #edge_in = (isnothing(dim_in) ? [] : edges[dim_in])
-    #edges_out = edges[dims_out]
     # slice isometries at this vertex
     Vv = [Vs[e] for e in edges]
 
     linkinds = [link_space[e] for e in edges]
+    linkdims = dim.(linkinds)
     
-
-    #####
-    # until here things are pretty straightforward at first sight
-    #####
-    #begin_block = Dict{Tuple{QN,Vector{Op}},Array{ValType,ranks[v]}}()
-    #cont_block = Dict{Tuple{QN,Vector{Op}},Array{ValType,ranks[v]}}()
-    #end_block = Dict{Tuple{QN,Vector{Op}},Array{ValType,ranks[v]}}()
-    #onsite_block = Dict{Tuple{QN,Vector{Op}},Array{ValType,ranks[v]}}()
-    #general_block = Dict{Tuple{QN,Vector{Op}},Array{ValType,ranks[v]}}()
+    # construct blocks
     blocks = Dict{Tuple{Block{ranks[v]},Vector{Op}},Array{ValType,ranks[v]}}()
-    ##construct blocks
-    linkdims = dim.(linkinds)
     for el in tempTTN[v]
       t = el.val
       (abs(coefficient(t)) > eps()) || continue
-      block_helper_inds=fill(-1,ranks[v])
+      block_helper_inds=fill(-1,ranks[v]) # we manipulate T_inds later, and loose track of ending/starting information, so keep track of it here
       T_inds=el.idxs
-      #@show T_inds
       T_qns=el.qn_idxs
       ct = convert(ValType, coefficient(t))
       sublinkdims= [(T_inds[i] == -1 ? 1 : qnblockdim(linkinds[i], T_qns[i])) for i in 1:ranks[v]]
@@ -457,30 +429,30 @@ function qn_svdTTN(
       T_inds[terminal_dims] .= 1                                  # start in channel 1  ###??
       block_helper_inds[terminal_dims].=1
       for dout in filter(d -> d ∈ terminal_dims, dims_out)
-        T_inds[dout] = sublinkdims[dout]
-        @assert T_inds[dout]==1                         # end in channel linkdims[d] for each dimension d
+        T_inds[dout] = sublinkdims[dout]                          # end in channel linkdims[d] for each dimension d
+        @assert T_inds[dout]==1                         
         block_helper_inds[dout]=nblocks(linkinds[dout])
       end
-      #non-trivial helper inds
+      
+      # set non-trivial helper inds
       for d in normal_dims
-        #if d in dims_out
         block_helper_inds[d]=qnblock(linkinds[d], T_qns[d])
       end
 
-      @assert all(block_helper_inds .!= -1)
+      @assert all(block_helper_inds .!= -1) # check that all block 
+      # make Block 
       theblock=Block(Tuple(block_helper_inds))
       
-      ### T_qns takes into account incoming/outgoing arrow direction, while Vs are stored for incoming arrow
-      ### flip sign of QN for outgoing arrows, to match incoming arrow definition of Vs
+      # T_qns takes into account incoming/outgoing arrow direction, while Vs are stored for incoming arrow
+      # flip sign of QN for outgoing arrows, to match incoming arrow definition of Vs
       iT_qns=deepcopy(T_qns)
       for d in dims_out
         iT_qns[d]=-iT_qns[d]
       end
 
       if isempty(normal_dims)
-
         M = get!(blocks, (theblock, terms(t)),zero_arr())
-        @assert product(size(M))==1 ###FIXME: figure out how to do this assert properly
+        @assert reduce((*),size(M))==1
         M[] += ct
       else
         M = get!(blocks, (theblock, terms(t)),zero_arr())
@@ -508,8 +480,6 @@ function qn_svdTTN(
     end
 
     for ((b,q_op),m) in blocks
-      #@show q_op
-      #op_prod = q_op[2]
       ###FIXME: make this work if there's no physical degree of freedom at site
       Op = computeSiteProd(sites, Prod(q_op))  ###FIXME: is this implemented?
       (nnzblocks(Op) == 0) && continue  ###FIXME: this one may be breaking for no physical indices on site
@@ -527,12 +497,11 @@ function qn_svdTTN(
           ###this logic requires dim_in to be 1 always, i don't think this is guaranteed
           ###however, maybe it is required when using autofermion?
           ###so we would have to figure out a perm to bring it into that format
-          
-          perm=(1,3,2)
-
-        end
+          ##compute perm factor to bring into this format, combine with perm factor that results from MPS-like logic
+          #perm=(1,3,2)
+          end
       end
-      T = ITensors.BlockSparseTensor(ValType, [b], _linkinds)  ###why the brackets around b? does the constructor dispatch on a vector of blocks?
+      T = ITensors.BlockSparseTensor(ValType, [b], _linkinds) 
       T[b] .= m
       H[v] += (itensor(T)*Op)
       
@@ -556,7 +525,7 @@ function qn_svdTTN(
         idT_end_inds[dout] = linkdims[dout] # reset
       end
     end
-    T = itensor(idT, _linkinds)  ###definitely need to dagger link_in , or do that beforehand?
+    T = itensor(idT, _linkinds)
     H[v] += T * ITensorNetworks.computeSiteProd(sites, Prod([(Op("Id", v))]))
   end
 
@@ -568,7 +537,7 @@ end
 # Tree adaptations of functionalities in ITensors.jl/src/physics/autompo/opsum_to_mpo_generic.jl
 # 
 
-# TODO: fix quantum number and fermion support, definitely broken
+# TODO: fix fermion support, definitely broken
 
 # needed an extra `only` compared to ITensors version since IndsNetwork has Vector{<:Index}
 # as vertex data
@@ -761,9 +730,9 @@ function Base.isless(a1::ArrElem{T,N}, a2::ArrElem{T,N})::Bool where {T,N}
   return a1.val < a2.val
 end
 
-#################################
+#####################################
 # QNArrElem (sparse array with QNs) #
-#################################
+#####################################
 
 struct QNArrElem{T,N}
   qn_idxs::MVector{N,QN}

From c62ed5c161c952e8839a91ecbe9b516105e52002 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Mon, 8 Jan 2024 15:32:29 -0500
Subject: [PATCH 03/32] Use on-site terms in test.

---
 test/test_opsum_to_ttn.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index 096c4e0c..f667f715 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -120,7 +120,7 @@ using Test
     J1 = -1
     J2 = 2
     h = 0.5
-    H = heisenberg(c; J1=J1, J2=J2)
+    H = heisenberg(c; J1=J1, J2=J2,h=h)
 
     # add combination of longer range interactions
     Hlr = copy(H)

From a9e84a1737af3ac3a66249667e990267e4e025db Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Mon, 8 Jan 2024 15:34:38 -0500
Subject: [PATCH 04/32] Format.

---
 src/treetensornetworks/opsum_to_ttn.jl | 105 +++++++++++++------------
 test/test_opsum_to_ttn.jl              |  24 ++----
 2 files changed, 62 insertions(+), 67 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 2e542a2f..6a2d780c 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -221,7 +221,6 @@ function svdTTN(
   return H
 end
 
-
 function qn_svdTTN(
   os::OpSum{C}, sites::IndsNetwork{VT,<:Index}, root_vertex::VT; kwargs...
 )::TTN where {C,VT}
@@ -241,21 +240,23 @@ function qn_svdTTN(
   #outmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   
   inmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()                  # map from term in Hamiltonian to incoming channel index for every edge
   outmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()                 # map from term in Hamiltonian to outgoing channel index for every edge
-  
+
   op_cache = Dict{Pair{String,VT},ITensor}()
   function calcQN(term::Vector{Op})
     q = QN()
     for st in term
       op_tensor = get(op_cache, ITensors.which_op(st) => ITensors.site(st), nothing)
       if op_tensor === nothing
-        op_tensor = op(sites[ITensors.site(st)], ITensors.which_op(st); ITensors.params(st)...)
+        op_tensor = op(
+          sites[ITensors.site(st)], ITensors.which_op(st); ITensors.params(st)...
+        )
         op_cache[ITensors.which_op(st) => ITensors.site(st)] = op_tensor
       end
       q -= flux(op_tensor)
     end
     return q
   end
-  
+
   Hflux = -calcQN(terms(first(terms(os))))
   inbond_coefs = Dict(e => Dict{QN,Vector{ITensors.MatElem{ValType}}}() 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
@@ -295,13 +296,13 @@ function qn_svdTTN(
         e => filter(t -> e ∈ edge_path(sites, v, ITensors.site(t)), ITensors.terms(term))
         for e in edges_out
       )
-      
+
       # compute QNs
-      incoming_qn=calcQN(incoming)
-      not_incoming_qn=calcQN(not_incoming)  ###does this work?
-      outgoing_qns=Dict(e => calcQN(outgoing[e]) for e in edges_out)
-      site_qn=calcQN(onsite)
-      
+      incoming_qn = calcQN(incoming)
+      not_incoming_qn = calcQN(not_incoming)  ###does this work?
+      outgoing_qns = Dict(e => calcQN(outgoing[e]) for e in edges_out)
+      site_qn = calcQN(onsite)
+
       # initialize QNArrayElement indices and quantum numbers 
       T_inds = MVector{ranks[v]}(fill(-1, ranks[v]))
       T_qns = MVector{ranks[v]}(fill(QN(), ranks[v]))
@@ -310,13 +311,15 @@ function qn_svdTTN(
       bond_col = -1
       if !isempty(incoming)
         # get the correct map from edge=>QN to term and channel
-        coutmap=get!(outmaps,edge_in=>-not_incoming_qn,Dict{Vector{Op},Int}())
-        cinmap=get!(inmaps,edge_in=>incoming_qn,Dict{Vector{Op},Int}())
+        coutmap = get!(outmaps, edge_in => -not_incoming_qn, Dict{Vector{Op},Int}())
+        cinmap = get!(inmaps, edge_in => incoming_qn, Dict{Vector{Op},Int}())
         # this checks if term exists on edge=>QN ( otherwise insert it) and returns it's index
-        bond_row = ITensors.posInLink!(cinmap, incoming)  
+        bond_row = ITensors.posInLink!(cinmap, incoming)
         bond_col = ITensors.posInLink!(coutmap, not_incoming) # get incoming channel
         bond_coef = convert(ValType, ITensors.coefficient(term))
-        q_inbond_coefs = get!(inbond_coefs[edge_in],incoming_qn, ITensors.MatElem{ValType}[])
+        q_inbond_coefs = get!(
+          inbond_coefs[edge_in], incoming_qn, ITensors.MatElem{ValType}[]
+        )
         push!(q_inbond_coefs, ITensors.MatElem(bond_row, bond_col, bond_coef))
         T_inds[dim_in] = bond_col
         T_qns[dim_in] = incoming_qn
@@ -324,7 +327,9 @@ function qn_svdTTN(
       for dout in dims_out
         #@show outgoing_qns[edges[dout]]
         #@show keys(outmaps)
-        coutmap=get!(outmaps,edges[dout]=>-outgoing_qns[edges[dout]],Dict{Vector{Op},Int}())
+        coutmap = get!(
+          outmaps, edges[dout] => -outgoing_qns[edges[dout]], Dict{Vector{Op},Int}()
+        )
         T_inds[dout] = ITensors.posInLink!(coutmap, outgoing[edges[dout]]) # add outgoing channel
         T_qns[dout] = -outgoing_qns[edges[dout]]  ###minus sign to account for outgoing arrow direction
       end
@@ -400,86 +405,86 @@ function qn_svdTTN(
     return error("Could not find block of QNIndex with matching QN")
   end
   qnblockdim(i::Index, q::QN) = blockdim(i, qnblock(i, q))
-  
+
   for v in vs
     # redo the whole thing like before
     edges = filter(e -> dst(e) == v || src(e) == v, es)
     dim_in = findfirst(e -> dst(e) == v, edges)
     dims_out = findall(e -> src(e) == v, edges)
-    
+
     # slice isometries at this vertex
     Vv = [Vs[e] for e in edges]
 
     linkinds = [link_space[e] for e in edges]
     linkdims = dim.(linkinds)
-    
+
     # construct blocks
     blocks = Dict{Tuple{Block{ranks[v]},Vector{Op}},Array{ValType,ranks[v]}}()
     for el in tempTTN[v]
       t = el.val
       (abs(coefficient(t)) > eps()) || continue
-      block_helper_inds=fill(-1,ranks[v]) # we manipulate T_inds later, and loose track of ending/starting information, so keep track of it here
-      T_inds=el.idxs
-      T_qns=el.qn_idxs
+      block_helper_inds = fill(-1, ranks[v]) # we manipulate T_inds later, and loose track of ending/starting information, so keep track of it here
+      T_inds = el.idxs
+      T_qns = el.qn_idxs
       ct = convert(ValType, coefficient(t))
-      sublinkdims= [(T_inds[i] == -1 ? 1 : qnblockdim(linkinds[i], T_qns[i])) for i in 1:ranks[v]]
+      sublinkdims = [
+        (T_inds[i] == -1 ? 1 : qnblockdim(linkinds[i], T_qns[i])) for i in 1:ranks[v]
+      ]
       zero_arr() = zeros(ValType, sublinkdims...)
       terminal_dims = findall(d -> T_inds[d] == -1, 1:ranks[v])   # directions in which term starts or ends
       normal_dims = findall(d -> T_inds[d] ≠ -1, 1:ranks[v])      # normal dimensions, do truncation thingies
       T_inds[terminal_dims] .= 1                                  # start in channel 1  ###??
-      block_helper_inds[terminal_dims].=1
+      block_helper_inds[terminal_dims] .= 1
       for dout in filter(d -> d ∈ terminal_dims, dims_out)
         T_inds[dout] = sublinkdims[dout]                          # end in channel linkdims[d] for each dimension d
-        @assert T_inds[dout]==1                         
-        block_helper_inds[dout]=nblocks(linkinds[dout])
+        @assert T_inds[dout] == 1
+        block_helper_inds[dout] = nblocks(linkinds[dout])
       end
-      
+
       # set non-trivial helper inds
       for d in normal_dims
-        block_helper_inds[d]=qnblock(linkinds[d], T_qns[d])
+        block_helper_inds[d] = qnblock(linkinds[d], T_qns[d])
       end
 
       @assert all(block_helper_inds .!= -1) # check that all block 
       # make Block 
-      theblock=Block(Tuple(block_helper_inds))
-      
+      theblock = Block(Tuple(block_helper_inds))
+
       # T_qns takes into account incoming/outgoing arrow direction, while Vs are stored for incoming arrow
       # flip sign of QN for outgoing arrows, to match incoming arrow definition of Vs
-      iT_qns=deepcopy(T_qns)
+      iT_qns = deepcopy(T_qns)
       for d in dims_out
-        iT_qns[d]=-iT_qns[d]
+        iT_qns[d] = -iT_qns[d]
       end
 
       if isempty(normal_dims)
-        M = get!(blocks, (theblock, terms(t)),zero_arr())
-        @assert reduce((*),size(M))==1
+        M = get!(blocks, (theblock, terms(t)), zero_arr())
+        @assert reduce((*), size(M)) == 1
         M[] += ct
       else
-        M = get!(blocks, (theblock, terms(t)),zero_arr())
-        dim_ranges = Tuple(size(Vv[d][iT_qns[d]],2) for d in normal_dims)
+        M = get!(blocks, (theblock, terms(t)), zero_arr())
+        dim_ranges = Tuple(size(Vv[d][iT_qns[d]], 2) for d in normal_dims)
         for c in CartesianIndices(dim_ranges)
           z = ct
           temp_inds = copy(T_inds)
-           for (i, d) in enumerate(normal_dims)            
+          for (i, d) in enumerate(normal_dims)
             V_factor = Vv[d][iT_qns[d]][T_inds[d], c[i]]
             z *= (d == dim_in ? conj(V_factor) : V_factor) # conjugate incoming isometry factor
             temp_inds[d] = c[i]
           end
-          M[temp_inds...] += z 
+          M[temp_inds...] += z
         end
       end
-
     end
-    
-    
+
     H[v] = ITensor()
 
-    _linkinds=copy(linkinds)
+    _linkinds = copy(linkinds)
     if !isnothing(dim_in)
-      _linkinds[dim_in]=dag(_linkinds[dim_in])
+      _linkinds[dim_in] = dag(_linkinds[dim_in])
     end
 
-    for ((b,q_op),m) in blocks
+    for ((b, q_op), m) in blocks
       ###FIXME: make this work if there's no physical degree of freedom at site
       Op = computeSiteProd(sites, Prod(q_op))  ###FIXME: is this implemented?
       (nnzblocks(Op) == 0) && continue  ###FIXME: this one may be breaking for no physical indices on site
@@ -499,15 +504,14 @@ function qn_svdTTN(
           ###so we would have to figure out a perm to bring it into that format
           ##compute perm factor to bring into this format, combine with perm factor that results from MPS-like logic
           #perm=(1,3,2)
-          end
+        end
       end
-      T = ITensors.BlockSparseTensor(ValType, [b], _linkinds) 
+      T = ITensors.BlockSparseTensor(ValType, [b], _linkinds)
       T[b] .= m
-      H[v] += (itensor(T)*Op)
-      
+      H[v] += (itensor(T) * Op)
     end
-    
-    linkdims=dim.(linkinds)
+
+    linkdims = dim.(linkinds)
     # add starting and ending identity operators
     idT = zeros(ValType, linkdims...)
     if isnothing(dim_in)
@@ -532,7 +536,6 @@ function qn_svdTTN(
   return H
 end
 
-
 # 
 # Tree adaptations of functionalities in ITensors.jl/src/physics/autompo/opsum_to_mpo_generic.jl
 # 
@@ -655,7 +658,7 @@ function TTN(
 
   if hasqns(first(first(vertex_data(sites))))
     ###added feature
-    T=qn_svdTTN(os,sites,root_vertex; kwargs...)
+    T = qn_svdTTN(os, sites, root_vertex; kwargs...)
     if splitblocks
       error("splitblocks not yet implemented for AbstractTreeTensorNetwork.")
       T = ITensors.splitblocks(linkinds, T) # TODO: make this work
diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index f667f715..871457a6 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -4,15 +4,13 @@ using ITensorNetworks
 using Random
 using Test
 
-
-
 @testset "OpSum to TTN converter" begin
   @testset "OpSum to TTN" begin
     # small comb tree
     tooth_lengths = fill(2, 3)
     c = named_comb_tree(tooth_lengths)
     @show c
-    
+
     is = siteinds("S=1/2", c)
 
     # linearized version
@@ -100,15 +98,15 @@ using Test
 
     @test H1 + H2 ≈ H3 rtol = 1e-6
   end
-  
+
   @testset "OpSum to TTN QN" begin
     # small comb tree
     tooth_lengths = fill(2, 3)
     c = named_comb_tree(tooth_lengths)
-    is = siteinds("S=1/2", c;conserve_qns=true)
-    is_noqns=copy(is)
+    is = siteinds("S=1/2", c; conserve_qns=true)
+    is_noqns = copy(is)
     for v in vertices(is)
-      is_noqns[v]=removeqns(is_noqns[v])
+      is_noqns[v] = removeqns(is_noqns[v])
     end
 
     # linearized version
@@ -120,7 +118,7 @@ using Test
     J1 = -1
     J2 = 2
     h = 0.5
-    H = heisenberg(c; J1=J1, J2=J2,h=h)
+    H = heisenberg(c; J1=J1, J2=J2, h=h)
 
     # add combination of longer range interactions
     Hlr = copy(H)
@@ -138,15 +136,13 @@ using Test
       # get corresponding MPO Hamiltonian
       Hline = MPO(relabel_sites(H, vmap), sites)
       # compare resulting sparse Hamiltonians
-      
-      
+
       @disable_warn_order begin
         Tmpo = prod(Hline)
         Tttno = contract(Hsvd)
-        
       end
       @test Tttno ≈ Tmpo rtol = 1e-6
-      
+
       # this breaks for longer range interactions ###not anymore
       Hsvd_lr = TTN(Hlr, is; root_vertex=root_vertex, method=:svd, cutoff=1e-10)
       Hline_lr = MPO(relabel_sites(Hlr, vmap), sites)
@@ -155,10 +151,6 @@ using Test
         Tmpo_lr = contract(Hsvd_lr)
       end
       @test Tttno_lr ≈ Tmpo_lr rtol = 1e-6
-      
     end
-    
-    
   end
-  
 end

From ada435b287dc6025ffa175e0ea60af6a38bc0363 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Mon, 8 Jan 2024 15:36:01 -0500
Subject: [PATCH 05/32] Remove debugging output from test.

---
 test/test_opsum_to_ttn.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index 871457a6..81abda57 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -9,7 +9,6 @@ using Test
     # small comb tree
     tooth_lengths = fill(2, 3)
     c = named_comb_tree(tooth_lengths)
-    @show c
 
     is = siteinds("S=1/2", c)
 

From eb556edb2774f3500d667907e6dfa6b1e62a6fc4 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Mon, 8 Jan 2024 15:57:48 -0500
Subject: [PATCH 06/32] Add function barrier, similar to (qn_)svdMPO in
 ITensors.

---
 src/treetensornetworks/opsum_to_ttn.jl | 26 +++++++++++++++++++++++++-
 1 file changed, 25 insertions(+), 1 deletion(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 6a2d780c..505a427f 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -41,6 +41,18 @@ Hamiltonian, compressing shared interaction channels.
 """
 function svdTTN(
   os::OpSum{C}, sites::IndsNetwork{VT,<:Index}, root_vertex::VT; kwargs...
+)::TTN where {C,VT}
+  # Function barrier to improve type stability
+  ValType = ITensors.determineValType(terms(os))
+  return svdTTN(ValType, os, sites, root_vertex; kwargs...)
+end
+
+function svdTTN(
+  ValType::Type{<:Number},
+  os::OpSum{C},
+  sites::IndsNetwork{VT,<:Index},
+  root_vertex::VT;
+  kwargs...,
 )::TTN where {C,VT}
   mindim::Int = get(kwargs, :mindim, 1)
   maxdim::Int = get(kwargs, :maxdim, 10000)
@@ -223,12 +235,24 @@ end
 
 function qn_svdTTN(
   os::OpSum{C}, sites::IndsNetwork{VT,<:Index}, root_vertex::VT; kwargs...
+)::TTN where {C,VT}
+  # Function barrier to improve type stability
+  ValType = ITensors.determineValType(terms(os))
+  return qn_svdTTN(ValType, os, sites, root_vertex; kwargs...)
+end
+
+function qn_svdTTN(
+  ValType::Type{<:Number},
+  os::OpSum{C},
+  sites::IndsNetwork{VT,<:Index},
+  root_vertex::VT;
+  kwargs...,
 )::TTN where {C,VT}
   mindim::Int = get(kwargs, :mindim, 1)
   maxdim::Int = get(kwargs, :maxdim, 10000)
   cutoff::Float64 = get(kwargs, :cutoff, 1e-15)
 
-  ValType = ITensors.determineValType(ITensors.terms(os))
+  #ValType = ITensors.determineValType(ITensors.terms(os))  #now included as argument in function signature
 
   # traverse tree outwards from root vertex
   vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root

From ee8a6c10f5a006703133e24d338efb7e11c80646 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Mon, 8 Jan 2024 16:36:52 -0500
Subject: [PATCH 07/32] Process both QN and non-QN OpSum-to-TTN conversion via
 qn_svdTTN.

---
 src/treetensornetworks/opsum_to_ttn.jl | 47 ++++++++++++++++++++------
 1 file changed, 37 insertions(+), 10 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 505a427f..02e9f114 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -248,11 +248,23 @@ function qn_svdTTN(
   root_vertex::VT;
   kwargs...,
 )::TTN where {C,VT}
+  
   mindim::Int = get(kwargs, :mindim, 1)
   maxdim::Int = get(kwargs, :maxdim, 10000)
   cutoff::Float64 = get(kwargs, :cutoff, 1e-15)
-
   #ValType = ITensors.determineValType(ITensors.terms(os))  #now included as argument in function signature
+  
+  # check for qns on the site indices
+  #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
+  #FIXME: how are sites handled where some sites have QNs and some don't?
+  
+  thishasqns=false
+  for site in vertices(sites)
+    if hasqns(sites[site])
+      thishasqns=true
+      break
+    end
+  end
 
   # traverse tree outwards from root vertex
   vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
@@ -276,7 +288,9 @@ function qn_svdTTN(
         )
         op_cache[ITensors.which_op(st) => ITensors.site(st)] = op_tensor
       end
-      q -= flux(op_tensor)
+      if !isnothing(flux(op_tensor))
+        q -= flux(op_tensor)
+      end
     end
     return q
   end
@@ -511,7 +525,9 @@ function qn_svdTTN(
     for ((b, q_op), m) in blocks
       ###FIXME: make this work if there's no physical degree of freedom at site
       Op = computeSiteProd(sites, Prod(q_op))  ###FIXME: is this implemented?
-      (nnzblocks(Op) == 0) && continue  ###FIXME: this one may be breaking for no physical indices on site
+      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
+        (nnzblocks(Op) == 0) && continue  ###FIXME: this one may be breaking for no physical indices on site
+      end
       sq = flux(Op)
       if !isnothing(sq)
         if ITensors.using_auto_fermion()
@@ -532,7 +548,12 @@ function qn_svdTTN(
       end
       T = ITensors.BlockSparseTensor(ValType, [b], _linkinds)
       T[b] .= m
-      H[v] += (itensor(T) * Op)
+      if !thishasqns
+        iT=removeqns(itensor(T))
+      else
+        iT=itensor(T)
+      end
+      H[v] += (iT * Op)
     end
 
     linkdims = dim.(linkinds)
@@ -554,6 +575,9 @@ function qn_svdTTN(
       end
     end
     T = itensor(idT, _linkinds)
+    if !thishasqns
+      T=removeqns(T)
+    end
     H[v] += T * ITensorNetworks.computeSiteProd(sites, Prod([(Op("Id", v))]))
   end
 
@@ -679,7 +703,9 @@ function TTN(
   os = deepcopy(os)
   os = sorteachterm(os, sites, root_vertex)
   os = ITensors.sortmergeterms(os) # not exported
-
+  
+  T = qn_svdTTN(os, sites, root_vertex; kwargs...)
+  #=
   if hasqns(first(first(vertex_data(sites))))
     ###added feature
     T = qn_svdTTN(os, sites, root_vertex; kwargs...)
@@ -689,11 +715,12 @@ function TTN(
     end
     return T
   end
-  T = svdTTN(os, sites, root_vertex; kwargs...)
-  if splitblocks
-    error("splitblocks not yet implemented for AbstractTreeTensorNetwork.")
-    T = ITensors.splitblocks(linkinds, T) # TODO: make this work
-  end
+  =#
+  #T = svdTTN(os, sites, root_vertex; kwargs...)
+  #if splitblocks
+  #  error("splitblocks not yet implemented for AbstractTreeTensorNetwork.")
+  #  T = ITensors.splitblocks(linkinds, T) # TODO: make this work
+  #end
   return T
 end
 

From 7f3de410774c71b0c70bf8b47fb1f979bfb071c7 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Mon, 8 Jan 2024 16:43:09 -0500
Subject: [PATCH 08/32] Minor cosmetic edits.

---
 src/treetensornetworks/opsum_to_ttn.jl | 25 +++++++++----------------
 1 file changed, 9 insertions(+), 16 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 02e9f114..71f5682a 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -248,23 +248,16 @@ function qn_svdTTN(
   root_vertex::VT;
   kwargs...,
 )::TTN where {C,VT}
-  
   mindim::Int = get(kwargs, :mindim, 1)
   maxdim::Int = get(kwargs, :maxdim, 10000)
   cutoff::Float64 = get(kwargs, :cutoff, 1e-15)
   #ValType = ITensors.determineValType(ITensors.terms(os))  #now included as argument in function signature
-  
+
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
   #FIXME: how are sites handled where some sites have QNs and some don't?
-  
-  thishasqns=false
-  for site in vertices(sites)
-    if hasqns(sites[site])
-      thishasqns=true
-      break
-    end
-  end
+
+  thishasqns = any(v -> hasqns(sites[v]), vertices(sites))
 
   # traverse tree outwards from root vertex
   vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
@@ -484,7 +477,7 @@ function qn_svdTTN(
         block_helper_inds[d] = qnblock(linkinds[d], T_qns[d])
       end
 
-      @assert all(block_helper_inds .!= -1) # check that all block 
+      @assert all(block_helper_inds .!= -1) # check that all block indices are set
       # make Block 
       theblock = Block(Tuple(block_helper_inds))
 
@@ -497,7 +490,7 @@ function qn_svdTTN(
 
       if isempty(normal_dims)
         M = get!(blocks, (theblock, terms(t)), zero_arr())
-        @assert reduce((*), size(M)) == 1
+        @assert length(M) == 1
         M[] += ct
       else
         M = get!(blocks, (theblock, terms(t)), zero_arr())
@@ -549,9 +542,9 @@ function qn_svdTTN(
       T = ITensors.BlockSparseTensor(ValType, [b], _linkinds)
       T[b] .= m
       if !thishasqns
-        iT=removeqns(itensor(T))
+        iT = removeqns(itensor(T))
       else
-        iT=itensor(T)
+        iT = itensor(T)
       end
       H[v] += (iT * Op)
     end
@@ -576,7 +569,7 @@ function qn_svdTTN(
     end
     T = itensor(idT, _linkinds)
     if !thishasqns
-      T=removeqns(T)
+      T = removeqns(T)
     end
     H[v] += T * ITensorNetworks.computeSiteProd(sites, Prod([(Op("Id", v))]))
   end
@@ -703,7 +696,7 @@ function TTN(
   os = deepcopy(os)
   os = sorteachterm(os, sites, root_vertex)
   os = ITensors.sortmergeterms(os) # not exported
-  
+
   T = qn_svdTTN(os, sites, root_vertex; kwargs...)
   #=
   if hasqns(first(first(vertex_data(sites))))

From 600364f39970d27c10a4090984cdff24fdba3111 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Mon, 8 Jan 2024 17:17:24 -0500
Subject: [PATCH 09/32] Accept suggested change for conversion to qnless
 ITensor.

---
 src/treetensornetworks/opsum_to_ttn.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 71f5682a..20b8511b 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -541,10 +541,9 @@ function qn_svdTTN(
       end
       T = ITensors.BlockSparseTensor(ValType, [b], _linkinds)
       T[b] .= m
+      iT = itensor(T)
       if !thishasqns
-        iT = removeqns(itensor(T))
-      else
-        iT = itensor(T)
+        iT = removeqns(iT)
       end
       H[v] += (iT * Op)
     end

From a5934c353f3d9230c25213ad178d39fa9933ecfc Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 15:27:15 -0500
Subject: [PATCH 10/32] Rename ValType to coefficient_type.

---
 src/treetensornetworks/opsum_to_ttn.jl | 80 +++++++++++++-------------
 1 file changed, 41 insertions(+), 39 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 20b8511b..992b4506 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -43,12 +43,12 @@ function svdTTN(
   os::OpSum{C}, sites::IndsNetwork{VT,<:Index}, root_vertex::VT; kwargs...
 )::TTN where {C,VT}
   # Function barrier to improve type stability
-  ValType = ITensors.determineValType(terms(os))
-  return svdTTN(ValType, os, sites, root_vertex; kwargs...)
+  coefficient_type = ITensors.determineValType(terms(os))
+  return svdTTN(coefficient_type, os, sites, root_vertex; kwargs...)
 end
 
 function svdTTN(
-  ValType::Type{<:Number},
+  coefficient_type::Type{<:Number},
   os::OpSum{C},
   sites::IndsNetwork{VT,<:Index},
   root_vertex::VT;
@@ -58,17 +58,17 @@ function svdTTN(
   maxdim::Int = get(kwargs, :maxdim, 10000)
   cutoff::Float64 = get(kwargs, :cutoff, 1e-15)
 
-  ValType = ITensors.determineValType(ITensors.terms(os))
+  coefficient_type = ITensors.determineValType(ITensors.terms(os))
 
   # traverse tree outwards from root vertex
   vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
   es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
   # some things to keep track of
   ranks = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
-  Vs = Dict(e => Matrix{ValType}(undef, 1, 1) for e in es)                                  # link isometries for SVD compression of TTN
+  Vs = Dict(e => Matrix{coefficient_type}(undef, 1, 1) for e in es)                                  # link isometries for SVD compression of TTN
   inmaps = Dict(e => Dict{Vector{Op},Int}() for e in es)                                    # map from term in Hamiltonian to incoming channel index for every edge
   outmaps = Dict(e => Dict{Vector{Op},Int}() for e in es)                                   # map from term in Hamiltonian to outgoing channel index for every edge
-  inbond_coefs = Dict(e => ITensors.MatElem{ValType}[] for e in es)                         # bond coefficients for incoming edge channels
+  inbond_coefs = Dict(e => ITensors.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
 
   # temporary symbolic representation of TTN Hamiltonian
@@ -114,7 +114,7 @@ function svdTTN(
       if !isempty(incoming)
         bond_row = ITensors.posInLink!(inmaps[edge_in], incoming)
         bond_col = ITensors.posInLink!(outmaps[edge_in], not_incoming) # get incoming channel
-        bond_coef = convert(ValType, ITensors.coefficient(term))
+        bond_coef = convert(coefficient_type, ITensors.coefficient(term))
         push!(inbond_coefs[edge_in], ITensors.MatElem(bond_row, bond_col, bond_coef))
         T_inds[dim_in] = bond_col
       end
@@ -149,7 +149,7 @@ function svdTTN(
       truncate!(P; maxdim=maxdim, cutoff=cutoff, mindim=mindim)
       tdim = length(P)
       nc = size(M, 2)
-      Vs[edges[dim_in]] = Matrix{ValType}(V[1:nc, 1:tdim])
+      Vs[edges[dim_in]] = Matrix{coefficient_type}(V[1:nc, 1:tdim])
     end
   end
 
@@ -180,8 +180,8 @@ function svdTTN(
       T_inds = el.idxs
       t = el.val
       (abs(coefficient(t)) > eps()) || continue
-      T = zeros(ValType, linkdims...)
-      ct = convert(ValType, coefficient(t))
+      T = zeros(coefficient_type, linkdims...)
+      ct = convert(coefficient_type, coefficient(t))
       terminal_dims = findall(d -> T_inds[d] == -1, 1:ranks[v])   # directions in which term starts or ends
       normal_dims = findall(d -> T_inds[d] ≠ -1, 1:ranks[v])      # normal dimensions, do truncation thingies
       T_inds[terminal_dims] .= 1                                  # start in channel 1
@@ -210,7 +210,7 @@ function svdTTN(
     end
 
     # add starting and ending identity operators
-    idT = zeros(ValType, linkdims...)
+    idT = zeros(coefficient_type, linkdims...)
     if isnothing(dim_in)
       idT[ones(Int, ranks[v])...] = 1.0 # only one real starting identity
     end
@@ -237,21 +237,21 @@ function qn_svdTTN(
   os::OpSum{C}, sites::IndsNetwork{VT,<:Index}, root_vertex::VT; kwargs...
 )::TTN where {C,VT}
   # Function barrier to improve type stability
-  ValType = ITensors.determineValType(terms(os))
-  return qn_svdTTN(ValType, os, sites, root_vertex; kwargs...)
+  coefficient_type = ITensors.determineValType(terms(os))
+  return qn_svdTTN(coefficient_type, os, sites, root_vertex; kwargs...)
 end
 
 function qn_svdTTN(
-  ValType::Type{<:Number},
+  coefficient_type::Type{<:Number},
   os::OpSum{C},
-  sites::IndsNetwork{VT,<:Index},
+  sites::IndsNetwork,
   root_vertex::VT;
   kwargs...,
 )::TTN where {C,VT}
   mindim::Int = get(kwargs, :mindim, 1)
   maxdim::Int = get(kwargs, :maxdim, 10000)
   cutoff::Float64 = get(kwargs, :cutoff, 1e-15)
-  #ValType = ITensors.determineValType(ITensors.terms(os))  #now included as argument in function signature
+  #coefficient_type = ITensors.determineValType(ITensors.terms(os))  #now included as argument in function signature
 
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
@@ -260,11 +260,12 @@ function qn_svdTTN(
   thishasqns = any(v -> hasqns(sites[v]), vertices(sites))
 
   # traverse tree outwards from root vertex
+  ### FIXME: generalize to internal vertices
   vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
   es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
   # some things to keep track of
   ranks = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
-  Vs = Dict(e => Dict{QN,Matrix{ValType}}() for e in es)                                    # link isometries for SVD compression of TTN
+  Vs = Dict(e => Dict{QN,Matrix{coefficient_type}}() for e in es)                                    # link isometries for SVD compression of TTN
   #inmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   
   #outmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   
   inmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()                  # map from term in Hamiltonian to incoming channel index for every edge
@@ -289,7 +290,7 @@ function qn_svdTTN(
   end
 
   Hflux = -calcQN(terms(first(terms(os))))
-  inbond_coefs = Dict(e => Dict{QN,Vector{ITensors.MatElem{ValType}}}() for e in es)                         # bond coefficients for incoming edge channels
+  inbond_coefs = Dict(e => Dict{QN,Vector{ITensors.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
   # temporary symbolic representation of TTN Hamiltonian
   tempTTN = Dict(v => QNArrElem{Scaled{C,Prod{Op}},ranks[v]}[] for v in vs)
@@ -347,9 +348,9 @@ function qn_svdTTN(
         # this checks if term exists on edge=>QN ( otherwise insert it) and returns it's index
         bond_row = ITensors.posInLink!(cinmap, incoming)
         bond_col = ITensors.posInLink!(coutmap, not_incoming) # get incoming channel
-        bond_coef = convert(ValType, ITensors.coefficient(term))
+        bond_coef = convert(coefficient_type, ITensors.coefficient(term))
         q_inbond_coefs = get!(
-          inbond_coefs[edge_in], incoming_qn, ITensors.MatElem{ValType}[]
+          inbond_coefs[edge_in], incoming_qn, ITensors.MatElem{coefficient_type}[]
         )
         push!(q_inbond_coefs, ITensors.MatElem(bond_row, bond_col, bond_coef))
         T_inds[dim_in] = bond_col
@@ -372,6 +373,7 @@ function qn_svdTTN(
         push!(site_coef_done, ITensors.argument(term))
       end
       # add onsite identity for interactions passing through vertex
+      ### FIXME: generalize to internal vertices, i.e. do not add Id if this is an internal site
       if isempty(onsite)
         if !ITensors.using_auto_fermion() && isfermionic(incoming, sites)
           error("No verified fermion support for automatic TTN constructor!")
@@ -393,7 +395,7 @@ function qn_svdTTN(
         truncate!(P; maxdim=maxdim, cutoff=cutoff, mindim=mindim)
         tdim = length(P)
         nc = size(M, 2) ###CHECK: verify whether this is still correct, diverges from opsum_to_mpo_qn
-        Vs[edges[dim_in]][q] = Matrix{ValType}(V[1:nc, 1:tdim]) ###CHECK: verify whether this is still correct, diverges from opsum_to_mpo_qn
+        Vs[edges[dim_in]][q] = Matrix{coefficient_type}(V[1:nc, 1:tdim]) ###CHECK: verify whether this is still correct, diverges from opsum_to_mpo_qn
       end
     end
   end
@@ -450,18 +452,18 @@ function qn_svdTTN(
     linkdims = dim.(linkinds)
 
     # construct blocks
-    blocks = Dict{Tuple{Block{ranks[v]},Vector{Op}},Array{ValType,ranks[v]}}()
+    blocks = Dict{Tuple{Block{ranks[v]},Vector{Op}},Array{coefficient_type,ranks[v]}}()
     for el in tempTTN[v]
       t = el.val
       (abs(coefficient(t)) > eps()) || continue
       block_helper_inds = fill(-1, ranks[v]) # we manipulate T_inds later, and loose track of ending/starting information, so keep track of it here
       T_inds = el.idxs
       T_qns = el.qn_idxs
-      ct = convert(ValType, coefficient(t))
+      ct = convert(coefficient_type, coefficient(t))
       sublinkdims = [
         (T_inds[i] == -1 ? 1 : qnblockdim(linkinds[i], T_qns[i])) for i in 1:ranks[v]
       ]
-      zero_arr() = zeros(ValType, sublinkdims...)
+      zero_arr() = zeros(coefficient_type, sublinkdims...)
       terminal_dims = findall(d -> T_inds[d] == -1, 1:ranks[v])   # directions in which term starts or ends
       normal_dims = findall(d -> T_inds[d] ≠ -1, 1:ranks[v])      # normal dimensions, do truncation thingies
       T_inds[terminal_dims] .= 1                                  # start in channel 1  ###??
@@ -539,7 +541,7 @@ function qn_svdTTN(
           #perm=(1,3,2)
         end
       end
-      T = ITensors.BlockSparseTensor(ValType, [b], _linkinds)
+      T = ITensors.BlockSparseTensor(coefficient_type, [b], _linkinds)
       T[b] .= m
       iT = itensor(T)
       if !thishasqns
@@ -550,7 +552,7 @@ function qn_svdTTN(
 
     linkdims = dim.(linkinds)
     # add starting and ending identity operators
-    idT = zeros(ValType, linkdims...)
+    idT = zeros(coefficient_type, linkdims...)
     if isnothing(dim_in)
       idT[ones(Int, ranks[v])...] = 1.0 # only one real starting identity
     end
@@ -829,10 +831,10 @@ function finite_state_machine(
   os = sorteachterm(os, sites, root_vertex)
   os = ITensors.sortmergeterms(os)
 
-  ValType = ITensors.determineValType(ITensors.terms(os))
+  coefficient_type = ITensors.determineValType(ITensors.terms(os))
 
   # sparse symbolic representation of the TTN Hamiltonian as a DataGraph of SparseArrays
-  sparseTTN = DataGraph{V,SparseArray{Sum{Scaled{ValType,Prod{Op}}}}}(
+  sparseTTN = DataGraph{V,SparseArray{Sum{Scaled{coefficient_type,Prod{Op}}}}}(
     underlying_graph(sites)
   )
 
@@ -847,7 +849,7 @@ function finite_state_machine(
 
   for v in vs
     # collect all nontrivial entries of the TTN tensor at vertex v
-    entries = Tuple{MVector{ranks[v],Int},Scaled{ValType,Prod{Op}}}[]
+    entries = Tuple{MVector{ranks[v],Int},Scaled{coefficient_type,Prod{Op}}}[]
 
     # for every vertex, find all edges that contain this vertex
     edges = filter(e -> dst(e) == v || src(e) == v, es)
@@ -912,7 +914,7 @@ function finite_state_machine(
     linkdims = Tuple([
       (isempty(linkmaps[e]) ? 0 : maximum(values(linkmaps[e]))) + 2 for e in edges
     ])
-    T = SparseArray{Sum{Scaled{ValType,Prod{Op}}},ranks[v]}(undef, linkdims)
+    T = SparseArray{Sum{Scaled{coefficient_type,Prod{Op}}},ranks[v]}(undef, linkdims)
     for (T_inds, t) in entries
       if !isnothing(dim_in)
         if T_inds[dim_in] == -1
@@ -931,17 +933,17 @@ function finite_state_machine(
     end
     # add starting and ending identity operators
     if isnothing(dim_in)
-      T[ones(Int, ranks[v])...] += one(ValType) * Prod([Op("Id", v)]) # only one real starting identity
+      T[ones(Int, ranks[v])...] += one(coefficient_type) * Prod([Op("Id", v)]) # only one real starting identity
     end
     # ending identities are a little more involved
     if !isnothing(dim_in)
-      T[linkdims...] += one(ValType) * Prod([Op("Id", v)]) # place identity if all channels end
+      T[linkdims...] += one(coefficient_type) * Prod([Op("Id", v)]) # place identity if all channels end
       # place identity from start of incoming channel to start of each single outgoing channel
       idT_end_inds = [linkdims...]
       idT_end_inds[dim_in] = 1
       for dout in dims_out
         idT_end_inds[dout] = 1
-        T[idT_end_inds...] += one(ValType) * Prod([Op("Id", v)])
+        T[idT_end_inds...] += one(coefficient_type) * Prod([Op("Id", v)])
         idT_end_inds[dout] = linkdims[dout] # reset
       end
     end
@@ -959,7 +961,7 @@ represenatation, without compression.
 function fsmTTN(
   os::OpSum{C}, sites::IndsNetwork{V,<:Index}, root_vertex::V; trunc=false, kwargs...
 )::TTN where {C,V}
-  ValType = ITensors.determineValType(ITensors.terms(os))
+  coefficient_type = ITensors.determineValType(ITensors.terms(os))
   # start from finite state machine
   fsm, edge_orders = finite_state_machine(os, sites, root_vertex)
   # some trickery to get link dimension for every edge
@@ -979,8 +981,8 @@ function fsmTTN(
     H[v] = ITensor()
     for (T_ind, t) in nonzero_pairs(fsm[v])
       any(map(x -> abs(coefficient(x)) > eps(), t)) || continue
-      T = zeros(ValType, linkdims...)
-      T[T_ind] += one(ValType)
+      T = zeros(coefficient_type, linkdims...)
+      T[T_ind] += one(coefficient_type)
       T = itensor(T, linkinds)
       H[v] += T * computeSiteSum(sites, t)
     end
@@ -1000,16 +1002,16 @@ end
 function computeSiteSum(
   sites::IndsNetwork{V,<:Index}, ops::Sum{Scaled{C,Prod{Op}}}
 )::ITensor where {V,C}
-  ValType = ITensors.determineValType(ITensors.terms(ops))
+  coefficient_type = ITensors.determineValType(ITensors.terms(ops))
   v = ITensors.site(ITensors.argument(ops[1])[1])
   T =
-    convert(ValType, coefficient(ops[1])) *
+    convert(coefficient_type, coefficient(ops[1])) *
     computeSiteProd(sites, ITensors.argument(ops[1]))
   for j in 2:length(ops)
     (ITensors.site(ITensors.argument(ops[j])[1]) != v) &&
       error("Mismatch of vertex labels in computeSiteSum")
     T +=
-      convert(ValType, coefficient(ops[j])) *
+      convert(coefficient_type, coefficient(ops[j])) *
       computeSiteProd(sites, ITensors.argument(ops[j]))
   end
   return T

From 684f14d6051dadc670f2d10c50dd604064a51826 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 15:32:35 -0500
Subject: [PATCH 11/32] Remove typing of IndsNetwork in function signatures..

---
 src/treetensornetworks/opsum_to_ttn.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 992b4506..34877269 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -40,7 +40,7 @@ Construct a dense TreeTensorNetwork from a symbolic OpSum representation of a
 Hamiltonian, compressing shared interaction channels.
 """
 function svdTTN(
-  os::OpSum{C}, sites::IndsNetwork{VT,<:Index}, root_vertex::VT; kwargs...
+  os::OpSum{C}, sites::IndsNetwork, root_vertex::VT; kwargs...
 )::TTN where {C,VT}
   # Function barrier to improve type stability
   coefficient_type = ITensors.determineValType(terms(os))
@@ -50,7 +50,7 @@ end
 function svdTTN(
   coefficient_type::Type{<:Number},
   os::OpSum{C},
-  sites::IndsNetwork{VT,<:Index},
+  sites::IndsNetwork,
   root_vertex::VT;
   kwargs...,
 )::TTN where {C,VT}
@@ -234,7 +234,7 @@ function svdTTN(
 end
 
 function qn_svdTTN(
-  os::OpSum{C}, sites::IndsNetwork{VT,<:Index}, root_vertex::VT; kwargs...
+  os::OpSum{C}, sites::IndsNetwork, root_vertex::VT; kwargs...
 )::TTN where {C,VT}
   # Function barrier to improve type stability
   coefficient_type = ITensors.determineValType(terms(os))

From 839de610e26c90701e9c9ea7d5d6c7b2b39a4bd8 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 15:48:53 -0500
Subject: [PATCH 12/32] Remove parametric type C (OpSum{C}) from function
 signatures and bodies.

---
 src/treetensornetworks/opsum_to_ttn.jl | 31 +++++++++++++-------------
 1 file changed, 16 insertions(+), 15 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 34877269..c3012033 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -40,8 +40,8 @@ Construct a dense TreeTensorNetwork from a symbolic OpSum representation of a
 Hamiltonian, compressing shared interaction channels.
 """
 function svdTTN(
-  os::OpSum{C}, sites::IndsNetwork, root_vertex::VT; kwargs...
-)::TTN where {C,VT}
+  os::OpSum, sites::IndsNetwork, root_vertex::VT; kwargs...
+)::TTN where {VT}
   # Function barrier to improve type stability
   coefficient_type = ITensors.determineValType(terms(os))
   return svdTTN(coefficient_type, os, sites, root_vertex; kwargs...)
@@ -49,11 +49,11 @@ end
 
 function svdTTN(
   coefficient_type::Type{<:Number},
-  os::OpSum{C},
+  os::OpSum,
   sites::IndsNetwork,
   root_vertex::VT;
   kwargs...,
-)::TTN where {C,VT}
+)::TTN where {VT}
   mindim::Int = get(kwargs, :mindim, 1)
   maxdim::Int = get(kwargs, :maxdim, 10000)
   cutoff::Float64 = get(kwargs, :cutoff, 1e-15)
@@ -72,7 +72,7 @@ function svdTTN(
   site_coef_done = Prod{Op}[]                                                               # list of terms for which the coefficient has been added to a site factor
 
   # temporary symbolic representation of TTN Hamiltonian
-  tempTTN = Dict(v => ArrElem{Scaled{C,Prod{Op}},ranks[v]}[] for v in vs)
+  tempTTN = Dict(v => ArrElem{Scaled{coefficient_type,Prod{Op}},ranks[v]}[] for v in vs)
 
   # build compressed finite state machine representation
   for v in vs
@@ -126,6 +126,8 @@ function svdTTN(
       if (isnothing(dim_in) || T_inds[dim_in] == -1) &&
         ITensors.argument(term) ∉ site_coef_done
         site_coef = ITensors.coefficient(term)
+        # FIXME: maybe put a try block here for a more helpful error message?
+        site_coef = convert(coefficient_type,site_coef) #required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real
         push!(site_coef_done, ITensors.argument(term))
       end
       # add onsite identity for interactions passing through vertex
@@ -197,7 +199,6 @@ function svdTTN(
           z = ct
           temp_inds = copy(T_inds)
           for (i, d) in enumerate(normal_dims)
-            #@show size(Vv[d]), T_inds[d],c[i]
             V_factor = Vv[d][T_inds[d], c[i]]
             z *= (d == dim_in ? conj(V_factor) : V_factor) # conjugate incoming isemetry factor
             temp_inds[d] = 1 + c[i]     ##the offset by one comes from the 1,...,1 padding (starting/ending...)
@@ -234,8 +235,8 @@ function svdTTN(
 end
 
 function qn_svdTTN(
-  os::OpSum{C}, sites::IndsNetwork, root_vertex::VT; kwargs...
-)::TTN where {C,VT}
+  os::OpSum, sites::IndsNetwork, root_vertex::VT; kwargs...
+)::TTN where {VT}
   # Function barrier to improve type stability
   coefficient_type = ITensors.determineValType(terms(os))
   return qn_svdTTN(coefficient_type, os, sites, root_vertex; kwargs...)
@@ -243,16 +244,15 @@ end
 
 function qn_svdTTN(
   coefficient_type::Type{<:Number},
-  os::OpSum{C},
+  os::OpSum,
   sites::IndsNetwork,
   root_vertex::VT;
   kwargs...,
-)::TTN where {C,VT}
+)::TTN where {VT}
   mindim::Int = get(kwargs, :mindim, 1)
   maxdim::Int = get(kwargs, :maxdim, 10000)
   cutoff::Float64 = get(kwargs, :cutoff, 1e-15)
   #coefficient_type = ITensors.determineValType(ITensors.terms(os))  #now included as argument in function signature
-
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
   #FIXME: how are sites handled where some sites have QNs and some don't?
@@ -293,7 +293,7 @@ function qn_svdTTN(
   inbond_coefs = Dict(e => Dict{QN,Vector{ITensors.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
   # temporary symbolic representation of TTN Hamiltonian
-  tempTTN = Dict(v => QNArrElem{Scaled{C,Prod{Op}},ranks[v]}[] for v in vs)
+  tempTTN = Dict(v => QNArrElem{Scaled{coefficient_type,Prod{Op}},ranks[v]}[] for v in vs)
 
   # build compressed finite state machine representation
   for v in vs
@@ -357,8 +357,6 @@ function qn_svdTTN(
         T_qns[dim_in] = incoming_qn
       end
       for dout in dims_out
-        #@show outgoing_qns[edges[dout]]
-        #@show keys(outmaps)
         coutmap = get!(
           outmaps, edges[dout] => -outgoing_qns[edges[dout]], Dict{Vector{Op},Int}()
         )
@@ -366,10 +364,12 @@ function qn_svdTTN(
         T_qns[dout] = -outgoing_qns[edges[dout]]  ###minus sign to account for outgoing arrow direction
       end
       # if term starts at this site, add its coefficient as a site factor
-      site_coef = one(C)
+      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)
+        # FIXME: maybe put a try block here for a more helpful error message?
+        site_coef = convert(coefficient_type,site_coef) #required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real
         push!(site_coef_done, ITensors.argument(term))
       end
       # add onsite identity for interactions passing through vertex
@@ -383,6 +383,7 @@ function qn_svdTTN(
       end
       # save indices and value of symbolic tensor entry
       el = QNArrElem(T_qns, T_inds, site_coef * Prod(onsite))
+      
       push!(tempTTN[v], el)
     end
     ITensors.remove_dups!(tempTTN[v])

From 8055abebfec54add9cdaa6fe39661c4aefc86d1b Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 16:13:53 -0500
Subject: [PATCH 13/32] Enumerate kwargs in function signature instead of
 getting them in function body.

---
 src/treetensornetworks/opsum_to_ttn.jl | 23 +++++++++++++----------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index c3012033..fc850c39 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -52,11 +52,11 @@ function svdTTN(
   os::OpSum,
   sites::IndsNetwork,
   root_vertex::VT;
-  kwargs...,
+  mindim::Int = 1,
+  maxdim::Int = 10000,      ###FIXME: better to use
+  cutoff::Float64 = 1e-15,  ###FIXME: better to use eps(coefficient_type) * 1e1 or something similar
 )::TTN where {VT}
-  mindim::Int = get(kwargs, :mindim, 1)
-  maxdim::Int = get(kwargs, :maxdim, 10000)
-  cutoff::Float64 = get(kwargs, :cutoff, 1e-15)
+  
 
   coefficient_type = ITensors.determineValType(ITensors.terms(os))
 
@@ -247,11 +247,10 @@ function qn_svdTTN(
   os::OpSum,
   sites::IndsNetwork,
   root_vertex::VT;
-  kwargs...,
+  mindim::Int = 1,
+  maxdim::Int = 10000,
+  cutoff::Float64 = 1e-15,  ###FIXME: better to use eps(coefficient_type) * 1e1 or something similar
 )::TTN where {VT}
-  mindim::Int = get(kwargs, :mindim, 1)
-  maxdim::Int = get(kwargs, :maxdim, 10000)
-  cutoff::Float64 = get(kwargs, :cutoff, 1e-15)
   #coefficient_type = ITensors.determineValType(ITensors.terms(os))  #now included as argument in function signature
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
@@ -689,6 +688,7 @@ function TTN(
   sites::IndsNetwork{V,<:Index};
   root_vertex::V=default_root_vertex(sites),
   splitblocks=false,
+  method=:svd,
   kwargs...,
 )::TTN where {V}
   length(ITensors.terms(os)) == 0 && error("OpSum has no terms")
@@ -698,8 +698,11 @@ function TTN(
   os = deepcopy(os)
   os = sorteachterm(os, sites, root_vertex)
   os = ITensors.sortmergeterms(os) # not exported
-
-  T = qn_svdTTN(os, sites, root_vertex; kwargs...)
+  if method==:svd
+    T = qn_svdTTN(os, sites, root_vertex; kwargs...)
+  else
+    error("Currently only SVD is supported as TTN constructor backend.")
+  end
   #=
   if hasqns(first(first(vertex_data(sites))))
     ###added feature

From 79a2e26d6398e2bfda12bc88cd5144a582400e7d Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 16:17:27 -0500
Subject: [PATCH 14/32] Renaming and cosmetic edits.

---
 src/treetensornetworks/opsum_to_ttn.jl | 58 +++++++++++++-------------
 1 file changed, 29 insertions(+), 29 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index fc850c39..9856a52c 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -64,7 +64,7 @@ function svdTTN(
   vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
   es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
   # some things to keep track of
-  ranks = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
+  degrees = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
   Vs = Dict(e => Matrix{coefficient_type}(undef, 1, 1) for e in es)                                  # link isometries for SVD compression of TTN
   inmaps = Dict(e => Dict{Vector{Op},Int}() for e in es)                                    # map from term in Hamiltonian to incoming channel index for every edge
   outmaps = Dict(e => Dict{Vector{Op},Int}() for e in es)                                   # map from term in Hamiltonian to outgoing channel index for every edge
@@ -72,7 +72,7 @@ function svdTTN(
   site_coef_done = Prod{Op}[]                                                               # list of terms for which the coefficient has been added to a site factor
 
   # temporary symbolic representation of TTN Hamiltonian
-  tempTTN = Dict(v => ArrElem{Scaled{coefficient_type,Prod{Op}},ranks[v]}[] for v in vs)
+  tempTTN = Dict(v => ArrElem{Scaled{coefficient_type,Prod{Op}},degrees[v]}[] for v in vs)
 
   # build compressed finite state machine representation
   for v in vs
@@ -108,7 +108,7 @@ function svdTTN(
       )
 
       # translate into tensor entry
-      T_inds = MVector{ranks[v]}(fill(-1, ranks[v]))
+      T_inds = MVector{degrees[v]}(fill(-1, degrees[v]))
       bond_row = -1
       bond_col = -1
       if !isempty(incoming)
@@ -148,7 +148,7 @@ function svdTTN(
       M = ITensors.toMatrix(inbond_coefs[edges[dim_in]])
       U, S, V = svd(M)
       P = S .^ 2
-      truncate!(P; maxdim=maxdim, cutoff=cutoff, mindim=mindim)
+      truncate!(P; maxdim, cutoff, mindim)
       tdim = length(P)
       nc = size(M, 2)
       Vs[edges[dim_in]] = Matrix{coefficient_type}(V[1:nc, 1:tdim])
@@ -184,8 +184,8 @@ function svdTTN(
       (abs(coefficient(t)) > eps()) || continue
       T = zeros(coefficient_type, linkdims...)
       ct = convert(coefficient_type, coefficient(t))
-      terminal_dims = findall(d -> T_inds[d] == -1, 1:ranks[v])   # directions in which term starts or ends
-      normal_dims = findall(d -> T_inds[d] ≠ -1, 1:ranks[v])      # normal dimensions, do truncation thingies
+      terminal_dims = findall(d -> T_inds[d] == -1, 1:degrees[v])   # directions in which term starts or ends
+      normal_dims = findall(d -> T_inds[d] ≠ -1, 1:degrees[v])      # normal dimensions, do truncation thingies
       T_inds[terminal_dims] .= 1                                  # start in channel 1
       for dout in filter(d -> d ∈ terminal_dims, dims_out)
         T_inds[dout] = linkdims[dout]                             # end in channel linkdims[d] for each dimension d
@@ -213,7 +213,7 @@ function svdTTN(
     # add starting and ending identity operators
     idT = zeros(coefficient_type, linkdims...)
     if isnothing(dim_in)
-      idT[ones(Int, ranks[v])...] = 1.0 # only one real starting identity
+      idT[ones(Int, degrees[v])...] = 1.0 # only one real starting identity
     end
     # ending identities are a little more involved
     if !isnothing(dim_in)
@@ -263,7 +263,7 @@ function qn_svdTTN(
   vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
   es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
   # some things to keep track of
-  ranks = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
+  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(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   
   #outmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   
@@ -271,7 +271,7 @@ function qn_svdTTN(
   outmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()                 # map from term in Hamiltonian to outgoing channel index for every edge
 
   op_cache = Dict{Pair{String,VT},ITensor}()
-  function calcQN(term::Vector{Op})
+  function calc_qn(term::Vector{Op})
     q = QN()
     for st in term
       op_tensor = get(op_cache, ITensors.which_op(st) => ITensors.site(st), nothing)
@@ -288,11 +288,11 @@ function qn_svdTTN(
     return q
   end
 
-  Hflux = -calcQN(terms(first(terms(os))))
+  Hflux = -calc_qn(terms(first(terms(os))))
   inbond_coefs = Dict(e => Dict{QN,Vector{ITensors.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
   # temporary symbolic representation of TTN Hamiltonian
-  tempTTN = Dict(v => QNArrElem{Scaled{coefficient_type,Prod{Op}},ranks[v]}[] for v in vs)
+  tempTTN = Dict(v => QNArrElem{Scaled{coefficient_type,Prod{Op}},degrees[v]}[] for v in vs)
 
   # build compressed finite state machine representation
   for v in vs
@@ -329,14 +329,14 @@ function qn_svdTTN(
       )
 
       # compute QNs
-      incoming_qn = calcQN(incoming)
-      not_incoming_qn = calcQN(not_incoming)  ###does this work?
-      outgoing_qns = Dict(e => calcQN(outgoing[e]) for e in edges_out)
-      site_qn = calcQN(onsite)
+      incoming_qn = calc_qn(incoming)
+      not_incoming_qn = calc_qn(not_incoming)  ###does this work?
+      outgoing_qns = Dict(e => calc_qn(outgoing[e]) for e in edges_out)
+      site_qn = calc_qn(onsite)
 
       # initialize QNArrayElement indices and quantum numbers 
-      T_inds = MVector{ranks[v]}(fill(-1, ranks[v]))
-      T_qns = MVector{ranks[v]}(fill(QN(), ranks[v]))
+      T_inds = MVector{degrees[v]}(fill(-1, degrees[v]))
+      T_qns = MVector{degrees[v]}(fill(QN(), degrees[v]))
       # initialize ArrayElement indices for inbond_coefs
       bond_row = -1
       bond_col = -1
@@ -392,7 +392,7 @@ function qn_svdTTN(
         M = ITensors.toMatrix(mat)
         U, S, V = svd(M)
         P = S .^ 2
-        truncate!(P; maxdim=maxdim, cutoff=cutoff, mindim=mindim)
+        truncate!(P; maxdim, cutoff, mindim)
         tdim = length(P)
         nc = size(M, 2) ###CHECK: verify whether this is still correct, diverges from opsum_to_mpo_qn
         Vs[edges[dim_in]][q] = Matrix{coefficient_type}(V[1:nc, 1:tdim]) ###CHECK: verify whether this is still correct, diverges from opsum_to_mpo_qn
@@ -452,20 +452,20 @@ function qn_svdTTN(
     linkdims = dim.(linkinds)
 
     # construct blocks
-    blocks = Dict{Tuple{Block{ranks[v]},Vector{Op}},Array{coefficient_type,ranks[v]}}()
+    blocks = Dict{Tuple{Block{degrees[v]},Vector{Op}},Array{coefficient_type,degrees[v]}}()
     for el in tempTTN[v]
       t = el.val
       (abs(coefficient(t)) > eps()) || continue
-      block_helper_inds = fill(-1, ranks[v]) # we manipulate T_inds later, and loose track of ending/starting information, so keep track of it here
+      block_helper_inds = fill(-1, degrees[v]) # we manipulate T_inds later, and loose track of ending/starting information, so keep track of it here
       T_inds = el.idxs
       T_qns = el.qn_idxs
       ct = convert(coefficient_type, coefficient(t))
       sublinkdims = [
-        (T_inds[i] == -1 ? 1 : qnblockdim(linkinds[i], T_qns[i])) for i in 1:ranks[v]
+        (T_inds[i] == -1 ? 1 : qnblockdim(linkinds[i], T_qns[i])) for i in 1:degrees[v]
       ]
       zero_arr() = zeros(coefficient_type, sublinkdims...)
-      terminal_dims = findall(d -> T_inds[d] == -1, 1:ranks[v])   # directions in which term starts or ends
-      normal_dims = findall(d -> T_inds[d] ≠ -1, 1:ranks[v])      # normal dimensions, do truncation thingies
+      terminal_dims = findall(d -> T_inds[d] == -1, 1:degrees[v])   # directions in which term starts or ends
+      normal_dims = findall(d -> T_inds[d] ≠ -1, 1:degrees[v])      # normal dimensions, do truncation thingies
       T_inds[terminal_dims] .= 1                                  # start in channel 1  ###??
       block_helper_inds[terminal_dims] .= 1
       for dout in filter(d -> d ∈ terminal_dims, dims_out)
@@ -554,7 +554,7 @@ function qn_svdTTN(
     # add starting and ending identity operators
     idT = zeros(coefficient_type, linkdims...)
     if isnothing(dim_in)
-      idT[ones(Int, ranks[v])...] = 1.0 # only one real starting identity
+      idT[ones(Int, degrees[v])...] = 1.0 # only one real starting identity
     end
     # ending identities are a little more involved
     if !isnothing(dim_in)
@@ -846,14 +846,14 @@ function finite_state_machine(
   vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
   es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
   # some things to keep track of
-  ranks = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
+  degrees = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
   linkmaps = Dict(e => Dict{Prod{Op},Int}() for e in es)                                    # map from term in Hamiltonian to edge channel index for every edge
   site_coef_done = Prod{Op}[]                                                               # list of Hamiltonian terms for which the coefficient has been added to a site factor
   edge_orders = DataGraph{V,Vector{edgetype(sites)}}(underlying_graph(sites))               # relate indices of sparse TTN tensor to incident graph edges for each site
 
   for v in vs
     # collect all nontrivial entries of the TTN tensor at vertex v
-    entries = Tuple{MVector{ranks[v],Int},Scaled{coefficient_type,Prod{Op}}}[]
+    entries = Tuple{MVector{degrees[v],Int},Scaled{coefficient_type,Prod{Op}}}[]
 
     # for every vertex, find all edges that contain this vertex
     edges = filter(e -> dst(e) == v || src(e) == v, es)
@@ -885,7 +885,7 @@ function finite_state_machine(
       )
 
       # translate into sparse tensor entry
-      T_inds = MVector{ranks[v]}(fill(-1, ranks[v]))
+      T_inds = MVector{degrees[v]}(fill(-1, degrees[v]))
       if !isnothing(dim_in) && !isempty(incoming)
         T_inds[dim_in] = ITensors.posInLink!(linkmaps[edge_in], ITensors.argument(term)) # get incoming channel
       end
@@ -918,7 +918,7 @@ function finite_state_machine(
     linkdims = Tuple([
       (isempty(linkmaps[e]) ? 0 : maximum(values(linkmaps[e]))) + 2 for e in edges
     ])
-    T = SparseArray{Sum{Scaled{coefficient_type,Prod{Op}}},ranks[v]}(undef, linkdims)
+    T = SparseArray{Sum{Scaled{coefficient_type,Prod{Op}}},degrees[v]}(undef, linkdims)
     for (T_inds, t) in entries
       if !isnothing(dim_in)
         if T_inds[dim_in] == -1
@@ -937,7 +937,7 @@ function finite_state_machine(
     end
     # add starting and ending identity operators
     if isnothing(dim_in)
-      T[ones(Int, ranks[v])...] += one(coefficient_type) * Prod([Op("Id", v)]) # only one real starting identity
+      T[ones(Int, degrees[v])...] += one(coefficient_type) * Prod([Op("Id", v)]) # only one real starting identity
     end
     # ending identities are a little more involved
     if !isnothing(dim_in)

From 027c33b73048d8707d02b24c2781bf3560950b02 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 16:46:39 -0500
Subject: [PATCH 15/32] Remove parametric type VT from signature. Changed
 OpCache key signature.

---
 src/treetensornetworks/opsum_to_ttn.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 9856a52c..a3075b6b 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -57,7 +57,7 @@ function svdTTN(
   cutoff::Float64 = 1e-15,  ###FIXME: better to use eps(coefficient_type) * 1e1 or something similar
 )::TTN where {VT}
   
-
+  edgetype_sites = edgetype(sites)
   coefficient_type = ITensors.determineValType(ITensors.terms(os))
 
   # traverse tree outwards from root vertex
@@ -251,11 +251,10 @@ function qn_svdTTN(
   maxdim::Int = 10000,
   cutoff::Float64 = 1e-15,  ###FIXME: better to use eps(coefficient_type) * 1e1 or something similar
 )::TTN where {VT}
-  #coefficient_type = ITensors.determineValType(ITensors.terms(os))  #now included as argument in function signature
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
   #FIXME: how are sites handled where some sites have QNs and some don't?
-
+  edgetype_sites = edgetype(sites)
   thishasqns = any(v -> hasqns(sites[v]), vertices(sites))
 
   # traverse tree outwards from root vertex
@@ -267,10 +266,11 @@ function qn_svdTTN(
   Vs = Dict(e => Dict{QN,Matrix{coefficient_type}}() for e in es)                                    # link isometries for SVD compression of TTN
   #inmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   
   #outmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   
-  inmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()                  # map from term in Hamiltonian to incoming channel index for every edge
-  outmaps = Dict{Pair{NamedGraphs.NamedEdge{VT},QN},Dict{Vector{Op},Int}}()                 # map from term in Hamiltonian to outgoing channel index for every edge
+  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
 
-  op_cache = Dict{Pair{String,VT},ITensor}()
+  op_cache = Dict{Pair{String,Any},ITensor}()
+  
   function calc_qn(term::Vector{Op})
     q = QN()
     for st in term
@@ -406,7 +406,7 @@ function qn_svdTTN(
   #  e => Index((isempty(outmaps[e]) ? 0 : size(Vs[e], 2)) + 2, edge_tag(e)) for e in es
   #])
   ###FIXME: that's an ugly constructor... ###QNIndex(undef) not defined
-  link_space = Dict{NamedGraphs.NamedEdge{VT},Index}()
+  link_space = Dict{edgetype_sites,Index}()
   linkdir = ITensors.using_auto_fermion() ? ITensors.In : ITensors.Out ###FIXME: ??
   for e in es ###this one might also be tricky, in case there's no Vq defined on the edge
     ###also tricky w.r.t. the link direction

From 28037f9693b1682ac6d0e4b4df72beb7d88afb20 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 16:48:27 -0500
Subject: [PATCH 16/32] Generalize kwarg defaults.

---
 src/treetensornetworks/opsum_to_ttn.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index a3075b6b..708f1702 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -53,8 +53,8 @@ function svdTTN(
   sites::IndsNetwork,
   root_vertex::VT;
   mindim::Int = 1,
-  maxdim::Int = 10000,      ###FIXME: better to use
-  cutoff::Float64 = 1e-15,  ###FIXME: better to use eps(coefficient_type) * 1e1 or something similar
+  maxdim::Int = typemax(Int),
+  cutoff::Float64 = eps(coefficient_type)*1e1
 )::TTN where {VT}
   
   edgetype_sites = edgetype(sites)
@@ -248,8 +248,8 @@ function qn_svdTTN(
   sites::IndsNetwork,
   root_vertex::VT;
   mindim::Int = 1,
-  maxdim::Int = 10000,
-  cutoff::Float64 = 1e-15,  ###FIXME: better to use eps(coefficient_type) * 1e1 or something similar
+  maxdim::Int = typemax(Int),
+  cutoff::Float64 = eps(coefficient_type)*1e1
 )::TTN where {VT}
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly

From a11a4e0c57dc19cf7bb2b931bf14c709247908ee Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 16:50:16 -0500
Subject: [PATCH 17/32] Format.

---
 src/treetensornetworks/opsum_to_ttn.jl | 31 +++++++++++++-------------
 1 file changed, 15 insertions(+), 16 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 708f1702..d1ae8984 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -39,9 +39,7 @@ end
 Construct a dense TreeTensorNetwork from a symbolic OpSum representation of a
 Hamiltonian, compressing shared interaction channels.
 """
-function svdTTN(
-  os::OpSum, sites::IndsNetwork, root_vertex::VT; kwargs...
-)::TTN where {VT}
+function svdTTN(os::OpSum, sites::IndsNetwork, root_vertex::VT; kwargs...)::TTN where {VT}
   # Function barrier to improve type stability
   coefficient_type = ITensors.determineValType(terms(os))
   return svdTTN(coefficient_type, os, sites, root_vertex; kwargs...)
@@ -52,11 +50,10 @@ function svdTTN(
   os::OpSum,
   sites::IndsNetwork,
   root_vertex::VT;
-  mindim::Int = 1,
-  maxdim::Int = typemax(Int),
-  cutoff::Float64 = eps(coefficient_type)*1e1
+  mindim::Int=1,
+  maxdim::Int=typemax(Int),
+  cutoff::Float64=eps(coefficient_type) * 1e1,
 )::TTN where {VT}
-  
   edgetype_sites = edgetype(sites)
   coefficient_type = ITensors.determineValType(ITensors.terms(os))
 
@@ -127,7 +124,7 @@ function svdTTN(
         ITensors.argument(term) ∉ site_coef_done
         site_coef = ITensors.coefficient(term)
         # FIXME: maybe put a try block here for a more helpful error message?
-        site_coef = convert(coefficient_type,site_coef) #required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real
+        site_coef = convert(coefficient_type, site_coef) #required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real
         push!(site_coef_done, ITensors.argument(term))
       end
       # add onsite identity for interactions passing through vertex
@@ -247,9 +244,9 @@ function qn_svdTTN(
   os::OpSum,
   sites::IndsNetwork,
   root_vertex::VT;
-  mindim::Int = 1,
-  maxdim::Int = typemax(Int),
-  cutoff::Float64 = eps(coefficient_type)*1e1
+  mindim::Int=1,
+  maxdim::Int=typemax(Int),
+  cutoff::Float64=eps(coefficient_type) * 1e1,
 )::TTN where {VT}
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
@@ -270,7 +267,7 @@ function qn_svdTTN(
   outmaps = Dict{Pair{edgetype_sites,QN},Dict{Vector{Op},Int}}()                 # map from term in Hamiltonian to outgoing channel index for every edge
 
   op_cache = Dict{Pair{String,Any},ITensor}()
-  
+
   function calc_qn(term::Vector{Op})
     q = QN()
     for st in term
@@ -289,7 +286,9 @@ function qn_svdTTN(
   end
 
   Hflux = -calc_qn(terms(first(terms(os))))
-  inbond_coefs = Dict(e => Dict{QN,Vector{ITensors.MatElem{coefficient_type}}}() for e in es)                         # bond coefficients for incoming edge channels
+  inbond_coefs = Dict(
+    e => Dict{QN,Vector{ITensors.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
   # temporary symbolic representation of TTN Hamiltonian
   tempTTN = Dict(v => QNArrElem{Scaled{coefficient_type,Prod{Op}},degrees[v]}[] for v in vs)
@@ -368,7 +367,7 @@ function qn_svdTTN(
         ITensors.argument(term) ∉ site_coef_done
         site_coef = ITensors.coefficient(term)
         # FIXME: maybe put a try block here for a more helpful error message?
-        site_coef = convert(coefficient_type,site_coef) #required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real
+        site_coef = convert(coefficient_type, site_coef) #required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real
         push!(site_coef_done, ITensors.argument(term))
       end
       # add onsite identity for interactions passing through vertex
@@ -382,7 +381,7 @@ function qn_svdTTN(
       end
       # save indices and value of symbolic tensor entry
       el = QNArrElem(T_qns, T_inds, site_coef * Prod(onsite))
-      
+
       push!(tempTTN[v], el)
     end
     ITensors.remove_dups!(tempTTN[v])
@@ -698,7 +697,7 @@ function TTN(
   os = deepcopy(os)
   os = sorteachterm(os, sites, root_vertex)
   os = ITensors.sortmergeterms(os) # not exported
-  if method==:svd
+  if method == :svd
     T = qn_svdTTN(os, sites, root_vertex; kwargs...)
   else
     error("Currently only SVD is supported as TTN constructor backend.")

From bc87d907b3013f3ccfe999ca4ad4baaa856d8777 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 16:53:22 -0500
Subject: [PATCH 18/32] use vertextype(sites) in op_cache key signature.

---
 src/treetensornetworks/opsum_to_ttn.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index d1ae8984..610c0c4a 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -252,6 +252,7 @@ function qn_svdTTN(
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
   #FIXME: how are sites handled where some sites have QNs and some don't?
   edgetype_sites = edgetype(sites)
+  vertextype_sites = vertextype(sites)
   thishasqns = any(v -> hasqns(sites[v]), vertices(sites))
 
   # traverse tree outwards from root vertex
@@ -266,7 +267,7 @@ function qn_svdTTN(
   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
 
-  op_cache = Dict{Pair{String,Any},ITensor}()
+  op_cache = Dict{Pair{String,vertextype_sites},ITensor}()
 
   function calc_qn(term::Vector{Op})
     q = QN()

From 07051d8107ea3e5c5f8bac8a9968e6876552d316 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 16:59:44 -0500
Subject: [PATCH 19/32] Fix kwarg defaults.

---
 src/treetensornetworks/opsum_to_ttn.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 610c0c4a..1e6ea3bd 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -52,7 +52,7 @@ function svdTTN(
   root_vertex::VT;
   mindim::Int=1,
   maxdim::Int=typemax(Int),
-  cutoff::Float64=eps(coefficient_type) * 1e1,
+  cutoff::eps(coefficient_type) * 10,
 )::TTN where {VT}
   edgetype_sites = edgetype(sites)
   coefficient_type = ITensors.determineValType(ITensors.terms(os))
@@ -246,7 +246,7 @@ function qn_svdTTN(
   root_vertex::VT;
   mindim::Int=1,
   maxdim::Int=typemax(Int),
-  cutoff::Float64=eps(coefficient_type) * 1e1,
+  cutoff::eps(coefficient_type) * 10,
 )::TTN where {VT}
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly

From b0cd3f64cff0f62570fcf32723370af9ffae5ff5 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 17:04:35 -0500
Subject: [PATCH 20/32] Actually remove parametric type VT from function
 signatures.

---
 src/treetensornetworks/opsum_to_ttn.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 1e6ea3bd..2a5d9f48 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -52,7 +52,7 @@ function svdTTN(
   root_vertex::VT;
   mindim::Int=1,
   maxdim::Int=typemax(Int),
-  cutoff::eps(coefficient_type) * 10,
+  cutoff=eps(coefficient_type) * 10,
 )::TTN where {VT}
   edgetype_sites = edgetype(sites)
   coefficient_type = ITensors.determineValType(ITensors.terms(os))
@@ -232,8 +232,8 @@ function svdTTN(
 end
 
 function qn_svdTTN(
-  os::OpSum, sites::IndsNetwork, root_vertex::VT; kwargs...
-)::TTN where {VT}
+  os::OpSum, sites::IndsNetwork, root_vertex; kwargs...
+)::TTN
   # Function barrier to improve type stability
   coefficient_type = ITensors.determineValType(terms(os))
   return qn_svdTTN(coefficient_type, os, sites, root_vertex; kwargs...)
@@ -243,11 +243,11 @@ function qn_svdTTN(
   coefficient_type::Type{<:Number},
   os::OpSum,
   sites::IndsNetwork,
-  root_vertex::VT;
+  root_vertex;
   mindim::Int=1,
   maxdim::Int=typemax(Int),
-  cutoff::eps(coefficient_type) * 10,
-)::TTN where {VT}
+  cutoff=eps(coefficient_type) * 10,
+)::TTN
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
   #FIXME: how are sites handled where some sites have QNs and some don't?

From ab4af0e53391f8deeb667c2a5c2fa1552443af34 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 17:09:56 -0500
Subject: [PATCH 21/32] Remove finite state machine code due to lack of support
 for QNs.

---
 src/treetensornetworks/opsum_to_ttn.jl | 180 -------------------------
 test/test_opsum_to_ttn.jl              |  87 ++++++++----
 2 files changed, 59 insertions(+), 208 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 2a5d9f48..cce18543 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -821,187 +821,7 @@ function Base.zero(::Type{S}) where {S<:Sum}
 end
 Base.zero(t::Sum) = zero(typeof(t))
 
-"""
-    finite_state_machine(os::OpSum{C}, sites::IndsNetwork{V,<:Index}, root_vertex::V) where {C,V}
-
-Finite state machine generator for ITensors.OpSum Hamiltonian defined on a tree graph. The
-site Index graph must be a tree graph, and the chosen root  must be a leaf vertex of this
-tree. Returns a DataGraph of SparseArrayKit.SparseArrays.
-"""
-function finite_state_machine(
-  os::OpSum{C}, sites::IndsNetwork{V,<:Index}, root_vertex::V
-) where {C,V}
-  os = deepcopy(os)
-  os = sorteachterm(os, sites, root_vertex)
-  os = ITensors.sortmergeterms(os)
-
-  coefficient_type = ITensors.determineValType(ITensors.terms(os))
-
-  # sparse symbolic representation of the TTN Hamiltonian as a DataGraph of SparseArrays
-  sparseTTN = DataGraph{V,SparseArray{Sum{Scaled{coefficient_type,Prod{Op}}}}}(
-    underlying_graph(sites)
-  )
-
-  # traverse tree outwards from root vertex
-  vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
-  es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
-  # some things to keep track of
-  degrees = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
-  linkmaps = Dict(e => Dict{Prod{Op},Int}() for e in es)                                    # map from term in Hamiltonian to edge channel index for every edge
-  site_coef_done = Prod{Op}[]                                                               # list of Hamiltonian terms for which the coefficient has been added to a site factor
-  edge_orders = DataGraph{V,Vector{edgetype(sites)}}(underlying_graph(sites))               # relate indices of sparse TTN tensor to incident graph edges for each site
-
-  for v in vs
-    # collect all nontrivial entries of the TTN tensor at vertex v
-    entries = Tuple{MVector{degrees[v],Int},Scaled{coefficient_type,Prod{Op}}}[]
-
-    # for every vertex, find all edges that contain this vertex
-    edges = filter(e -> dst(e) == v || src(e) == v, es)
-    # use the corresponding ordering as index order for tensor elements at this site
-    edge_orders[v] = edges
-    dim_in = findfirst(e -> dst(e) == v, edges)
-    edge_in = (isnothing(dim_in) ? [] : edges[dim_in])
-    dims_out = findall(e -> src(e) == v, edges)
-    edges_out = edges[dims_out]
-
-    # sanity check, leaves only have single incoming or outgoing edge
-    @assert !isempty(dims_out) || !isnothing(dim_in)
-    (isempty(dims_out) || isnothing(dim_in)) && @assert is_leaf(sites, v)
-
-    for term in os
-      # loop over OpSum and pick out terms that act on current vertex
-      crosses_vertex(term, sites, v) || continue
-
-      # filter out factors that come in from the direction of the incoming edge
-      incoming = filter(
-        t -> edge_in ∈ edge_path(sites, ITensors.site(t), v), ITensors.terms(term)
-      )
-      # filter out factor that acts on current vertex
-      onsite = filter(t -> (ITensors.site(t) == v), ITensors.terms(term))
-      # for every outgoing edge, filter out factors that go out along that edge
-      outgoing = Dict(
-        e => filter(t -> e ∈ edge_path(sites, v, ITensors.site(t)), ITensors.terms(term))
-        for e in edges_out
-      )
-
-      # translate into sparse tensor entry
-      T_inds = MVector{degrees[v]}(fill(-1, degrees[v]))
-      if !isnothing(dim_in) && !isempty(incoming)
-        T_inds[dim_in] = ITensors.posInLink!(linkmaps[edge_in], ITensors.argument(term)) # get incoming channel
-      end
-      for dout in dims_out
-        if !isempty(outgoing[edges[dout]])
-          T_inds[dout] = ITensors.posInLink!(linkmaps[edges[dout]], ITensors.argument(term)) # add outgoing channel
-        end
-      end
-      # if term starts at this site, add its coefficient as a site factor
-      site_coef = one(C)
-      if (isnothing(dim_in) || T_inds[dim_in] == -1) &&
-        ITensors.argument(term) ∉ site_coef_done
-        site_coef = ITensors.coefficient(term)
-        push!(site_coef_done, ITensors.argument(term))
-      end
-      # add onsite identity for interactions passing through vertex
-      if isempty(onsite)
-        if !ITensors.using_auto_fermion() && isfermionic(incoming, sites)
-          error("No verified fermion support for automatic TTN constructor!") # no verified support, just throw error
-        else
-          push!(onsite, Op("Id", v))
-        end
-      end
-      # save indices and value of sparse tensor entry
-      el = (T_inds, site_coef * Prod(onsite))
-      push!(entries, el)
-    end
-
-    # handle start and end of operator terms and convert to sparse array
-    linkdims = Tuple([
-      (isempty(linkmaps[e]) ? 0 : maximum(values(linkmaps[e]))) + 2 for e in edges
-    ])
-    T = SparseArray{Sum{Scaled{coefficient_type,Prod{Op}}},degrees[v]}(undef, linkdims)
-    for (T_inds, t) in entries
-      if !isnothing(dim_in)
-        if T_inds[dim_in] == -1
-          T_inds[dim_in] = 1 # always start in first channel
-        else
-          T_inds[dim_in] += 1 # shift regular channel
-        end
-      end
-      if !isempty(dims_out)
-        end_dims = filter(d -> T_inds[d] == -1, dims_out)
-        normal_dims = filter(d -> T_inds[d] != -1, dims_out)
-        T_inds[end_dims] .= linkdims[end_dims] # always end in last channel
-        T_inds[normal_dims] .+= 1 # shift regular channels
-      end
-      T[T_inds...] += t
-    end
-    # add starting and ending identity operators
-    if isnothing(dim_in)
-      T[ones(Int, degrees[v])...] += one(coefficient_type) * Prod([Op("Id", v)]) # only one real starting identity
-    end
-    # ending identities are a little more involved
-    if !isnothing(dim_in)
-      T[linkdims...] += one(coefficient_type) * Prod([Op("Id", v)]) # place identity if all channels end
-      # place identity from start of incoming channel to start of each single outgoing channel
-      idT_end_inds = [linkdims...]
-      idT_end_inds[dim_in] = 1
-      for dout in dims_out
-        idT_end_inds[dout] = 1
-        T[idT_end_inds...] += one(coefficient_type) * Prod([Op("Id", v)])
-        idT_end_inds[dout] = linkdims[dout] # reset
-      end
-    end
-    sparseTTN[v] = T
-  end
-  return sparseTTN, edge_orders
-end
-
-"""
-    fsmTTN(os::OpSum{C}, sites::IndsNetwork{V,<:Index}, root_vertex::V, kwargs...) where {C,V}
 
-Construct a dense TreeTensorNetwork from sparse finite state machine
-represenatation, without compression.
-"""
-function fsmTTN(
-  os::OpSum{C}, sites::IndsNetwork{V,<:Index}, root_vertex::V; trunc=false, kwargs...
-)::TTN where {C,V}
-  coefficient_type = ITensors.determineValType(ITensors.terms(os))
-  # start from finite state machine
-  fsm, edge_orders = finite_state_machine(os, sites, root_vertex)
-  # some trickery to get link dimension for every edge
-  link_space = Dict{edgetype(sites),Index}()
-  function get_linkind!(link_space, e)
-    if !haskey(link_space, e)
-      d = findfirst(x -> (x == e || x == reverse(e)), edge_orders[src(e)])
-      link_space[e] = Index(size(fsm[src(e)], d), edge_tag(e))
-    end
-    return link_space[e]
-  end
-  # compress finite state machine into dense form
-  H = TTN(sites)
-  for v in vertices(sites)
-    linkinds = [get_linkind!(link_space, e) for e in edge_orders[v]]
-    linkdims = dim.(linkinds)
-    H[v] = ITensor()
-    for (T_ind, t) in nonzero_pairs(fsm[v])
-      any(map(x -> abs(coefficient(x)) > eps(), t)) || continue
-      T = zeros(coefficient_type, linkdims...)
-      T[T_ind] += one(coefficient_type)
-      T = itensor(T, linkinds)
-      H[v] += T * computeSiteSum(sites, t)
-    end
-  end
-  # add option for numerical truncation, but throw warning as this can fail sometimes
-  if trunc
-    @warn "Naive numerical truncation of TTN Hamiltonian may fail for larger systems."
-    # see https://github.com/ITensor/ITensors.jl/issues/526
-    lognormT = lognorm(H)
-    H /= exp(lognormT / nv(H)) # TODO: fix broadcasting for in-place assignment
-    H = truncate(H; root_vertex, kwargs...)
-    H *= exp(lognormT / nv(H))
-  end
-  return H
-end
 
 function computeSiteSum(
   sites::IndsNetwork{V,<:Index}, ops::Sum{Scaled{C,Prod{Op}}}
diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index 81abda57..03b91b2e 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -22,7 +22,6 @@ using Test
     J2 = 2
     h = 0.5
     H = ising(c; J1=J1, J2=J2, h=h)
-
     # add combination of longer range interactions
     Hlr = copy(H)
     Hlr += 5, "Z", (1, 2), "Z", (2, 2)
@@ -54,32 +53,6 @@ using Test
       end
       @test Tttno_lr ≈ Tmpo_lr rtol = 1e-6
     end
-
-    @testset "Finite state machine" for root_vertex in leaf_vertices(is)
-      # get TTN Hamiltonian directly
-      Hfsm = ITensorNetworks.fsmTTN(H, is, root_vertex)
-      # get corresponding MPO Hamiltonian
-      Hline = MPO(relabel_sites(H, vmap), sites)
-      # compare resulting dense Hamiltonians
-      @disable_warn_order begin
-        Tttno = prod(Hline)
-        Tmpo = contract(Hfsm)
-      end
-      @test Tttno ≈ Tmpo rtol = 1e-6
-
-      # same thing for longer range interactions
-      Hfsm_lr = ITensorNetworks.fsmTTN(Hlr, is, root_vertex)
-      Hline_lr = MPO(relabel_sites(Hlr, vmap), sites)
-      @disable_warn_order begin
-        Tttno_lr = prod(Hline_lr)
-        Tmpo_lr = contract(Hfsm_lr)
-      end
-      @test Tttno_lr ≈ Tmpo_lr rtol = 1e-6
-
-      # check optional numerical truncation for finite state machine construction
-      Hfsm_trunc = ITensorNetworks.fsmTTN(H, is, root_vertex; trunc=true, cutoff=1e-10)
-      @test collect(edge_data(linkdims(Hfsm_trunc))) == [4, 3, 4, 3, 3]
-    end
   end
 
   @testset "Multiple onsite terms (regression test for issue #62)" begin
@@ -118,7 +91,6 @@ using Test
     J2 = 2
     h = 0.5
     H = heisenberg(c; J1=J1, J2=J2, h=h)
-
     # add combination of longer range interactions
     Hlr = copy(H)
     Hlr += 5, "Z", (1, 2), "Z", (2, 2)
@@ -152,4 +124,63 @@ using Test
       @test Tttno_lr ≈ Tmpo_lr rtol = 1e-6
     end
   end
+  #=
+    @testset "OpSum to TTN QN missing" begin
+      # small comb tree
+      tooth_lengths = fill(2, 3)
+      c = named_comb_tree(tooth_lengths)
+      is = siteinds("S=1/2", c; conserve_qns=true)
+      is_missing_site=copy(is)
+      is_missing_site[(2,1)]=Vector{Index}[]
+      is_noqns = copy(is)
+      for v in vertices(is)
+        is_noqns[v] = removeqns(is_noqns[v])
+      end
+
+      # linearized version
+      linear_order = [4, 1, 2, 5, 3, 6]
+      vmap = Dictionary(vertices(is)[linear_order], 1:length(linear_order))
+      sites = only.(collect(vertex_data(is)))[linear_order]
+
+      # test with next-to-nearest-neighbor Ising Hamiltonian
+      J1 = -1
+      J2 = 2
+      h = 0.5
+      H = heisenberg(c; J1=J1, J2=J2, h=h)
+      Hvac = heisenberg(is_missing_site; J1=J1, J2=J2, h=h)
+
+      # add combination of longer range interactions
+      Hlr = copy(H)
+      Hlr += 5, "Z", (1, 2), "Z", (2, 2)
+      Hlr += -4, "Z", (1, 1), "Z", (2, 2)
+      Hlr += 2.0, "Z", (2, 2), "Z", (3, 2)
+      Hlr += -1.0, "Z", (1, 2), "Z", (3, 1)
+
+      # root_vertex = (1, 2)
+      # println(leaf_vertices(is))
+
+      @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)
+        # get corresponding MPO Hamiltonian
+        Hline = MPO(relabel_sites(H, vmap), sites)
+        # compare resulting sparse Hamiltonians
+
+        @disable_warn_order begin
+          Tmpo = prod(Hline)
+          Tttno = contract(Hsvd)
+        end
+        @test Tttno ≈ Tmpo rtol = 1e-6
+
+        # this breaks for longer range interactions ###not anymore
+        Hsvd_lr = TTN(Hlr, is; root_vertex=root_vertex, method=:svd, cutoff=1e-10)
+        Hline_lr = MPO(relabel_sites(Hlr, vmap), sites)
+        @disable_warn_order begin
+          Tttno_lr = prod(Hline_lr)
+          Tmpo_lr = contract(Hsvd_lr)
+        end
+        @test Tttno_lr ≈ Tmpo_lr rtol = 1e-6
+      end
+    end
+    =#
 end

From fe75daca062b271984fa29ccd6111b6b2c1baf25 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 17:14:06 -0500
Subject: [PATCH 22/32] Remove deprecated code from opsum_to_ttn.jl. Also adapt
 TTN constructor signature to no parametric types.

---
 src/treetensornetworks/opsum_to_ttn.jl | 227 ++-----------------------
 1 file changed, 12 insertions(+), 215 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index cce18543..61dcc70c 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -34,212 +34,20 @@ end
 # 
 
 """
-    svdTTN(os::OpSum{C}, sites::IndsNetwork{V<:Index}, root_vertex::V, kwargs...) where {C,V}
+    ttn_svd(os::OpSum, sites::IndsNetwork, root_vertex, kwargs...)
 
 Construct a dense TreeTensorNetwork from a symbolic OpSum representation of a
 Hamiltonian, compressing shared interaction channels.
 """
-function svdTTN(os::OpSum, sites::IndsNetwork, root_vertex::VT; kwargs...)::TTN where {VT}
-  # Function barrier to improve type stability
-  coefficient_type = ITensors.determineValType(terms(os))
-  return svdTTN(coefficient_type, os, sites, root_vertex; kwargs...)
-end
-
-function svdTTN(
-  coefficient_type::Type{<:Number},
-  os::OpSum,
-  sites::IndsNetwork,
-  root_vertex::VT;
-  mindim::Int=1,
-  maxdim::Int=typemax(Int),
-  cutoff=eps(coefficient_type) * 10,
-)::TTN where {VT}
-  edgetype_sites = edgetype(sites)
-  coefficient_type = ITensors.determineValType(ITensors.terms(os))
-
-  # traverse tree outwards from root vertex
-  vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
-  es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
-  # 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 => Matrix{coefficient_type}(undef, 1, 1) for e in es)                                  # link isometries for SVD compression of TTN
-  inmaps = Dict(e => Dict{Vector{Op},Int}() for e in es)                                    # map from term in Hamiltonian to incoming channel index for every edge
-  outmaps = Dict(e => Dict{Vector{Op},Int}() for e in es)                                   # map from term in Hamiltonian to outgoing channel index for every edge
-  inbond_coefs = Dict(e => ITensors.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
-
-  # temporary symbolic representation of TTN Hamiltonian
-  tempTTN = Dict(v => ArrElem{Scaled{coefficient_type,Prod{Op}},degrees[v]}[] for v in vs)
-
-  # build compressed finite state machine representation
-  for v in vs
-    # for every vertex, find all edges that contain this vertex
-    edges = filter(e -> dst(e) == v || src(e) == v, es)
-    # use the corresponding ordering as index order for tensor elements at this site
-    dim_in = findfirst(e -> dst(e) == v, edges)
-    edge_in = (isnothing(dim_in) ? [] : edges[dim_in])
-    dims_out = findall(e -> src(e) == v, edges)
-    edges_out = edges[dims_out]
-    # sanity check, leaves only have single incoming or outgoing edge
-    @assert !isempty(dims_out) || !isnothing(dim_in)
-    (isempty(dims_out) || isnothing(dim_in)) && @assert is_leaf(sites, v)
-
-    for term in os
-      # loop over OpSum and pick out terms that act on current vertex
-      crosses_vertex(term, sites, v) || continue
-
-      # filter out factors that come in from the direction of the incoming edge
-      incoming = filter(
-        t -> edge_in ∈ edge_path(sites, ITensors.site(t), v), ITensors.terms(term)
-      )
-      # also store all non-incoming factors in standard order, used for channel merging
-      not_incoming = filter(
-        t -> edge_in ∉ edge_path(sites, ITensors.site(t), v), ITensors.terms(term)
-      )
-      # filter out factor that acts on current vertex
-      onsite = filter(t -> (ITensors.site(t) == v), ITensors.terms(term))
-      # for every outgoing edge, filter out factors that go out along that edge
-      outgoing = Dict(
-        e => filter(t -> e ∈ edge_path(sites, v, ITensors.site(t)), ITensors.terms(term))
-        for e in edges_out
-      )
-
-      # translate into tensor entry
-      T_inds = MVector{degrees[v]}(fill(-1, degrees[v]))
-      bond_row = -1
-      bond_col = -1
-      if !isempty(incoming)
-        bond_row = ITensors.posInLink!(inmaps[edge_in], incoming)
-        bond_col = ITensors.posInLink!(outmaps[edge_in], not_incoming) # get incoming channel
-        bond_coef = convert(coefficient_type, ITensors.coefficient(term))
-        push!(inbond_coefs[edge_in], ITensors.MatElem(bond_row, bond_col, bond_coef))
-        T_inds[dim_in] = bond_col
-      end
-      for dout in dims_out
-        T_inds[dout] = ITensors.posInLink!(outmaps[edges[dout]], outgoing[edges[dout]]) # add outgoing channel
-      end
-      # if term starts at this site, add its coefficient as a site factor
-      site_coef = one(C)
-      if (isnothing(dim_in) || T_inds[dim_in] == -1) &&
-        ITensors.argument(term) ∉ site_coef_done
-        site_coef = ITensors.coefficient(term)
-        # FIXME: maybe put a try block here for a more helpful error message?
-        site_coef = convert(coefficient_type, site_coef) #required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real
-        push!(site_coef_done, ITensors.argument(term))
-      end
-      # add onsite identity for interactions passing through vertex
-      if isempty(onsite)
-        if !ITensors.using_auto_fermion() && isfermionic(incoming, sites)
-          error("No verified fermion support for automatic TTN constructor!")
-        else
-          push!(onsite, Op("Id", v))
-        end
-      end
-      # save indices and value of symbolic tensor entry
-      el = ArrElem(T_inds, site_coef * Prod(onsite))
-      push!(tempTTN[v], el)
-    end
-    ITensors.remove_dups!(tempTTN[v])
-    # manual truncation: isometry on incoming edge
-    if !isnothing(dim_in) && !isempty(inbond_coefs[edges[dim_in]])
-      M = ITensors.toMatrix(inbond_coefs[edges[dim_in]])
-      U, S, V = svd(M)
-      P = S .^ 2
-      truncate!(P; maxdim, cutoff, mindim)
-      tdim = length(P)
-      nc = size(M, 2)
-      Vs[edges[dim_in]] = Matrix{coefficient_type}(V[1:nc, 1:tdim])
-    end
-  end
-
-  # compress this tempTTN representation into dense form
-
-  link_space = dictionary([
-    e => Index((isempty(outmaps[e]) ? 0 : size(Vs[e], 2)) + 2, edge_tag(e)) for e in es
-  ])
-
-  H = TTN(sites)
-
-  for v in vs
-
-    # redo the whole thing like before
-    edges = filter(e -> dst(e) == v || src(e) == v, es)
-    dim_in = findfirst(e -> dst(e) == v, edges)
-    dims_out = findall(e -> src(e) == v, edges)
-
-    # slice isometries at this vertex
-    Vv = [Vs[e] for e in edges]
-
-    linkinds = [link_space[e] for e in edges]
-    linkdims = dim.(linkinds)
-
-    H[v] = ITensor()
-
-    for el in tempTTN[v]
-      T_inds = el.idxs
-      t = el.val
-      (abs(coefficient(t)) > eps()) || continue
-      T = zeros(coefficient_type, linkdims...)
-      ct = convert(coefficient_type, coefficient(t))
-      terminal_dims = findall(d -> T_inds[d] == -1, 1:degrees[v])   # directions in which term starts or ends
-      normal_dims = findall(d -> T_inds[d] ≠ -1, 1:degrees[v])      # normal dimensions, do truncation thingies
-      T_inds[terminal_dims] .= 1                                  # start in channel 1
-      for dout in filter(d -> d ∈ terminal_dims, dims_out)
-        T_inds[dout] = linkdims[dout]                             # end in channel linkdims[d] for each dimension d
-      end
-      if isempty(normal_dims)
-        T[T_inds...] += ct                                        # on-site term
-      else
-        # handle channel compression isometries
-        dim_ranges = Tuple(size(Vv[d], 2) for d in normal_dims)
-        for c in CartesianIndices(dim_ranges)
-          z = ct
-          temp_inds = copy(T_inds)
-          for (i, d) in enumerate(normal_dims)
-            V_factor = Vv[d][T_inds[d], c[i]]
-            z *= (d == dim_in ? conj(V_factor) : V_factor) # conjugate incoming isemetry factor
-            temp_inds[d] = 1 + c[i]     ##the offset by one comes from the 1,...,1 padding (starting/ending...)
-          end
-          T[temp_inds...] += z
-        end
-      end
-      T = itensor(T, linkinds)
-      H[v] += T * computeSiteProd(sites, ITensors.argument(t))
-    end
-
-    # 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
-    end
-    # ending identities are a little more involved
-    if !isnothing(dim_in)
-      idT[linkdims...] = 1.0 # place identity if all channels end
-      # 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.0
-      for dout in dims_out
-        idT_end_inds[dout] = 1.0
-        idT[idT_end_inds...] = 1.0
-        idT_end_inds[dout] = linkdims[dout] # reset
-      end
-    end
-    T = itensor(idT, linkinds)
-    H[v] += T * ITensorNetworks.computeSiteProd(sites, Prod([Op("Id", v)]))
-  end
-
-  return H
-end
-
-function qn_svdTTN(
+function ttn_svd(
   os::OpSum, sites::IndsNetwork, root_vertex; kwargs...
 )::TTN
   # Function barrier to improve type stability
   coefficient_type = ITensors.determineValType(terms(os))
-  return qn_svdTTN(coefficient_type, os, sites, root_vertex; kwargs...)
+  return ttn_svd(coefficient_type, os, sites, root_vertex; kwargs...)
 end
 
-function qn_svdTTN(
+function ttn_svd(
   coefficient_type::Type{<:Number},
   os::OpSum,
   sites::IndsNetwork,
@@ -685,12 +493,12 @@ Convert an OpSum object `os` to a TreeTensorNetwork, with indices given by `site
 """
 function TTN(
   os::OpSum,
-  sites::IndsNetwork{V,<:Index};
-  root_vertex::V=default_root_vertex(sites),
+  sites::IndsNetwork;
+  root_vertex=default_root_vertex(sites),
   splitblocks=false,
   method=:svd,
   kwargs...,
-)::TTN where {V}
+)::TTN
   length(ITensors.terms(os)) == 0 && error("OpSum has no terms")
   is_tree(sites) || error("Site index graph must be a tree.")
   is_leaf(sites, root_vertex) || error("Tree root must be a leaf vertex.")
@@ -699,26 +507,15 @@ function TTN(
   os = sorteachterm(os, sites, root_vertex)
   os = ITensors.sortmergeterms(os) # not exported
   if method == :svd
-    T = qn_svdTTN(os, sites, root_vertex; kwargs...)
+    T = ttn_svd(os, sites, root_vertex; kwargs...)
   else
     error("Currently only SVD is supported as TTN constructor backend.")
   end
-  #=
-  if hasqns(first(first(vertex_data(sites))))
-    ###added feature
-    T = qn_svdTTN(os, sites, root_vertex; kwargs...)
-    if splitblocks
-      error("splitblocks not yet implemented for AbstractTreeTensorNetwork.")
-      T = ITensors.splitblocks(linkinds, T) # TODO: make this work
-    end
-    return T
+
+  if splitblocks
+    error("splitblocks not yet implemented for AbstractTreeTensorNetwork.")
+    T = ITensors.splitblocks(linkinds, T) # TODO: make this work
   end
-  =#
-  #T = svdTTN(os, sites, root_vertex; kwargs...)
-  #if splitblocks
-  #  error("splitblocks not yet implemented for AbstractTreeTensorNetwork.")
-  #  T = ITensors.splitblocks(linkinds, T) # TODO: make this work
-  #end
   return T
 end
 

From 764f4daff693d8e8d5fb6c5e19b68c023d982109 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 17:19:56 -0500
Subject: [PATCH 23/32] Restore old version as deprecated_opsum_to_ttn.jl

---
 .../deprecated_opsum_to_ttn.jl                | 636 ++++++++++++++++++
 1 file changed, 636 insertions(+)
 create mode 100644 src/treetensornetworks/deprecated_opsum_to_ttn.jl

diff --git a/src/treetensornetworks/deprecated_opsum_to_ttn.jl b/src/treetensornetworks/deprecated_opsum_to_ttn.jl
new file mode 100644
index 00000000..c0113c18
--- /dev/null
+++ b/src/treetensornetworks/deprecated_opsum_to_ttn.jl
@@ -0,0 +1,636 @@
+### This version is deprecated, we are keeping it around for later reference.
+# convert ITensors.OpSum to TreeTensorNetwork
+
+# 
+# Utility methods
+# 
+
+# linear ordering of vertices in tree graph relative to chosen root, chosen outward from root
+function find_index_in_tree(site, g::AbstractGraph, root_vertex)
+  ordering = reverse(post_order_dfs_vertices(g, root_vertex))
+  return findfirst(x -> x == site, ordering)
+end
+function find_index_in_tree(o::Op, g::AbstractGraph, root_vertex)
+  return find_index_in_tree(ITensors.site(o), g, root_vertex)
+end
+
+# determine 'support' of product operator on tree graph
+function span(t::Scaled{C,Prod{Op}}, g::AbstractGraph) where {C}
+  spn = eltype(g)[]
+  nterms = length(t)
+  for i in 1:nterms, j in i:nterms
+    path = vertex_path(g, ITensors.site(t[i]), ITensors.site(t[j]))
+    spn = union(spn, path)
+  end
+  return spn
+end
+
+# determine whether an operator string crosses a given graph vertex
+function crosses_vertex(t::Scaled{C,Prod{Op}}, g::AbstractGraph, v) where {C}
+  return v ∈ span(t, g)
+end
+
+# 
+# Tree adaptations of functionalities in ITensors.jl/src/physics/autompo/opsum_to_mpo.jl
+# 
+
+"""
+    svdTTN(os::OpSum{C}, sites::IndsNetwork{V<:Index}, root_vertex::V, kwargs...) where {C,V}
+
+Construct a dense TreeTensorNetwork from a symbolic OpSum representation of a
+Hamiltonian, compressing shared interaction channels.
+"""
+function svdTTN(
+  os::OpSum{C}, sites::IndsNetwork{VT,<:Index}, root_vertex::VT; kwargs...
+)::TTN where {C,VT}
+  mindim::Int = get(kwargs, :mindim, 1)
+  maxdim::Int = get(kwargs, :maxdim, 10000)
+  cutoff::Float64 = get(kwargs, :cutoff, 1e-15)
+
+  ValType = ITensors.determineValType(ITensors.terms(os))
+
+  # traverse tree outwards from root vertex
+  vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
+  es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
+  # some things to keep track of
+  ranks = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
+  Vs = Dict(e => Matrix{ValType}(undef, 1, 1) for e in es)                                  # link isometries for SVD compression of TTN
+  inmaps = Dict(e => Dict{Vector{Op},Int}() for e in es)                                    # map from term in Hamiltonian to incoming channel index for every edge
+  outmaps = Dict(e => Dict{Vector{Op},Int}() for e in es)                                   # map from term in Hamiltonian to outgoing channel index for every edge
+  inbond_coefs = Dict(e => ITensors.MatElem{ValType}[] 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
+
+  # temporary symbolic representation of TTN Hamiltonian
+  tempTTN = Dict(v => ArrElem{Scaled{C,Prod{Op}},ranks[v]}[] for v in vs)
+
+  # build compressed finite state machine representation
+  for v in vs
+    # for every vertex, find all edges that contain this vertex
+    edges = filter(e -> dst(e) == v || src(e) == v, es)
+    # use the corresponding ordering as index order for tensor elements at this site
+    dim_in = findfirst(e -> dst(e) == v, edges)
+    edge_in = (isnothing(dim_in) ? [] : edges[dim_in])
+    dims_out = findall(e -> src(e) == v, edges)
+    edges_out = edges[dims_out]
+
+    # sanity check, leaves only have single incoming or outgoing edge
+    @assert !isempty(dims_out) || !isnothing(dim_in)
+    (isempty(dims_out) || isnothing(dim_in)) && @assert is_leaf(sites, v)
+
+    for term in os
+      # loop over OpSum and pick out terms that act on current vertex
+      crosses_vertex(term, sites, v) || continue
+
+      # filter out factors that come in from the direction of the incoming edge
+      incoming = filter(
+        t -> edge_in ∈ edge_path(sites, ITensors.site(t), v), ITensors.terms(term)
+      )
+      # also store all non-incoming factors in standard order, used for channel merging
+      not_incoming = filter(
+        t -> edge_in ∉ edge_path(sites, ITensors.site(t), v), ITensors.terms(term)
+      )
+      # filter out factor that acts on current vertex
+      onsite = filter(t -> (ITensors.site(t) == v), ITensors.terms(term))
+      # for every outgoing edge, filter out factors that go out along that edge
+      outgoing = Dict(
+        e => filter(t -> e ∈ edge_path(sites, v, ITensors.site(t)), ITensors.terms(term))
+        for e in edges_out
+      )
+
+      # translate into tensor entry
+      T_inds = MVector{ranks[v]}(fill(-1, ranks[v]))
+      bond_row = -1
+      bond_col = -1
+      if !isempty(incoming)
+        bond_row = ITensors.posInLink!(inmaps[edge_in], incoming)
+        bond_col = ITensors.posInLink!(outmaps[edge_in], not_incoming) # get incoming channel
+        bond_coef = convert(ValType, ITensors.coefficient(term))
+        push!(inbond_coefs[edge_in], ITensors.MatElem(bond_row, bond_col, bond_coef))
+        T_inds[dim_in] = bond_col
+      end
+      for dout in dims_out
+        T_inds[dout] = ITensors.posInLink!(outmaps[edges[dout]], outgoing[edges[dout]]) # add outgoing channel
+      end
+      # if term starts at this site, add its coefficient as a site factor
+      site_coef = one(C)
+      if (isnothing(dim_in) || T_inds[dim_in] == -1) &&
+        ITensors.argument(term) ∉ site_coef_done
+        site_coef = ITensors.coefficient(term)
+        push!(site_coef_done, ITensors.argument(term))
+      end
+      # add onsite identity for interactions passing through vertex
+      if isempty(onsite)
+        if !ITensors.using_auto_fermion() && isfermionic(incoming, sites)
+          error("No verified fermion support for automatic TTN constructor!")
+        else
+          push!(onsite, Op("Id", v))
+        end
+      end
+      # save indices and value of symbolic tensor entry
+      el = ArrElem(T_inds, site_coef * Prod(onsite))
+      push!(tempTTN[v], el)
+    end
+    ITensors.remove_dups!(tempTTN[v])
+    # manual truncation: isometry on incoming edge
+    if !isnothing(dim_in) && !isempty(inbond_coefs[edges[dim_in]])
+      M = ITensors.toMatrix(inbond_coefs[edges[dim_in]])
+      U, S, V = svd(M)
+      P = S .^ 2
+      truncate!(P; maxdim=maxdim, cutoff=cutoff, mindim=mindim)
+      tdim = length(P)
+      nc = size(M, 2)
+      Vs[edges[dim_in]] = Matrix{ValType}(V[1:nc, 1:tdim])
+    end
+  end
+
+  # compress this tempTTN representation into dense form
+
+  link_space = dictionary([
+    e => Index((isempty(outmaps[e]) ? 0 : size(Vs[e], 2)) + 2, edge_tag(e)) for e in es
+  ])
+
+  H = TTN(sites)
+
+  for v in vs
+
+    # redo the whole thing like before
+    edges = filter(e -> dst(e) == v || src(e) == v, es)
+    dim_in = findfirst(e -> dst(e) == v, edges)
+    dims_out = findall(e -> src(e) == v, edges)
+
+    # slice isometries at this vertex
+    Vv = [Vs[e] for e in edges]
+
+    linkinds = [link_space[e] for e in edges]
+    linkdims = dim.(linkinds)
+
+    H[v] = ITensor()
+
+    for el in tempTTN[v]
+      T_inds = el.idxs
+      t = el.val
+      (abs(coefficient(t)) > eps()) || continue
+      T = zeros(ValType, linkdims...)
+      ct = convert(ValType, coefficient(t))
+      terminal_dims = findall(d -> T_inds[d] == -1, 1:ranks[v])   # directions in which term starts or ends
+      normal_dims = findall(d -> T_inds[d] ≠ -1, 1:ranks[v])      # normal dimensions, do truncation thingies
+      T_inds[terminal_dims] .= 1                                  # start in channel 1
+      for dout in filter(d -> d ∈ terminal_dims, dims_out)
+        T_inds[dout] = linkdims[dout]                             # end in channel linkdims[d] for each dimension d
+      end
+      if isempty(normal_dims)
+        T[T_inds...] += ct                                        # on-site term
+      else
+        # handle channel compression isometries
+        dim_ranges = Tuple(size(Vv[d], 2) for d in normal_dims)
+        for c in CartesianIndices(dim_ranges)
+          z = ct
+          temp_inds = copy(T_inds)
+          for (i, d) in enumerate(normal_dims)
+            V_factor = Vv[d][T_inds[d], c[i]]
+            z *= (d == dim_in ? conj(V_factor) : V_factor) # conjugate incoming isemetry factor
+            temp_inds[d] = 1 + c[i]
+          end
+          T[temp_inds...] += z
+        end
+      end
+      T = itensor(T, linkinds)
+      H[v] += T * computeSiteProd(sites, ITensors.argument(t))
+    end
+
+    # add starting and ending identity operators
+    idT = zeros(ValType, linkdims...)
+    if isnothing(dim_in)
+      idT[ones(Int, ranks[v])...] = 1.0 # only one real starting identity
+    end
+    # ending identities are a little more involved
+    if !isnothing(dim_in)
+      idT[linkdims...] = 1.0 # place identity if all channels end
+      # 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.0
+      for dout in dims_out
+        idT_end_inds[dout] = 1.0
+        idT[idT_end_inds...] = 1.0
+        idT_end_inds[dout] = linkdims[dout] # reset
+      end
+    end
+    T = itensor(idT, linkinds)
+    H[v] += T * ITensorNetworks.computeSiteProd(sites, Prod([Op("Id", v)]))
+  end
+
+  return H
+end
+
+# 
+# Tree adaptations of functionalities in ITensors.jl/src/physics/autompo/opsum_to_mpo_generic.jl
+# 
+
+# TODO: fix quantum number and fermion support, definitely broken
+
+# needed an extra `only` compared to ITensors version since IndsNetwork has Vector{<:Index}
+# as vertex data
+function isfermionic(t::Vector{Op}, sites::IndsNetwork{V,<:Index}) where {V}
+  p = +1
+  for op in t
+    if has_fermion_string(ITensors.name(op), only(sites[ITensors.site(op)]))
+      p *= -1
+    end
+  end
+  return (p == -1)
+end
+
+# only(site(ops[1])) in ITensors breaks for Tuple site labels, had to drop the only
+function computeSiteProd(sites::IndsNetwork{V,<:Index}, ops::Prod{Op})::ITensor where {V}
+  v = ITensors.site(ops[1])
+  T = op(sites[v], ITensors.which_op(ops[1]); ITensors.params(ops[1])...)
+  for j in 2:length(ops)
+    (ITensors.site(ops[j]) != v) && error("Mismatch of vertex labels in computeSiteProd")
+    opj = op(sites[v], ITensors.which_op(ops[j]); ITensors.params(ops[j])...)
+    T = product(T, opj)
+  end
+  return T
+end
+
+# changed `isless_site` to use tree vertex ordering relative to root
+function sorteachterm(os::OpSum, sites::IndsNetwork{V,<:Index}, root_vertex::V) where {V}
+  os = copy(os)
+  findpos(op::Op) = find_index_in_tree(op, sites, root_vertex)
+  isless_site(o1::Op, o2::Op) = findpos(o1) < findpos(o2)
+  N = nv(sites)
+  for n in eachindex(os)
+    t = os[n]
+    Nt = length(t)
+
+    if !all(map(v -> has_vertex(sites, v), ITensors.sites(t)))
+      error(
+        "The OpSum contains a term $t that does not have support on the underlying graph."
+      )
+    end
+
+    prevsite = N + 1 #keep track of whether we are switching
+    #to a new site to make sure F string
+    #is only placed at most once for each site
+
+    # Sort operators in t by site order,
+    # and keep the permutation used, perm, for analysis below
+    perm = Vector{Int}(undef, Nt)
+    sortperm!(perm, ITensors.terms(t); alg=InsertionSort, lt=isless_site)
+
+    t = coefficient(t) * Prod(ITensors.terms(t)[perm])
+
+    # Identify fermionic operators,
+    # zeroing perm for bosonic operators,
+    # and inserting string "F" operators
+    parity = +1
+    for n in Nt:-1:1
+      currsite = ITensors.site(t[n])
+      fermionic = has_fermion_string(
+        ITensors.which_op(t[n]), only(sites[ITensors.site(t[n])])
+      )
+      if !ITensors.using_auto_fermion() && (parity == -1) && (currsite < prevsite)
+        error("No verified fermion support for automatic TTN constructor!") # no verified support, just throw error
+        # Put local piece of Jordan-Wigner string emanating
+        # from fermionic operators to the right
+        # (Remaining F operators will be put in by svdMPO)
+        terms(t)[n] = Op("$(ITensors.which_op(t[n])) * F", only(ITensors.site(t[n])))
+      end
+      prevsite = currsite
+
+      if fermionic
+        error("No verified fermion support for automatic TTN constructor!") # no verified support, just throw error
+        parity = -parity
+      else
+        # Ignore bosonic operators in perm
+        # by zeroing corresponding entries
+        perm[n] = 0
+      end
+    end
+    if parity == -1
+      error("Parity-odd fermionic terms not yet supported by AutoTTN")
+    end
+
+    # Keep only fermionic op positions (non-zero entries)
+    filter!(!iszero, perm)
+    # and account for anti-commuting, fermionic operators 
+    # during above sort; put resulting sign into coef
+    t *= ITensors.parity_sign(perm)
+    ITensors.terms(os)[n] = t
+  end
+  return os
+end
+
+"""
+    TTN(os::OpSum, sites::IndsNetwork{<:Index}; kwargs...)
+    TTN(eltype::Type{<:Number}, os::OpSum, sites::IndsNetwork{<:Index}; kwargs...)
+       
+Convert an OpSum object `os` to a TreeTensorNetwork, with indices given by `sites`.
+"""
+function TTN(
+  os::OpSum,
+  sites::IndsNetwork{V,<:Index};
+  root_vertex::V=default_root_vertex(sites),
+  splitblocks=false,
+  kwargs...,
+)::TTN where {V}
+  length(ITensors.terms(os)) == 0 && error("OpSum has no terms")
+  is_tree(sites) || error("Site index graph must be a tree.")
+  is_leaf(sites, root_vertex) || error("Tree root must be a leaf vertex.")
+
+  os = deepcopy(os)
+  os = sorteachterm(os, sites, root_vertex)
+  os = ITensors.sortmergeterms(os) # not exported
+
+  if hasqns(first(first(vertex_data(sites))))
+    if !is_path_graph(sites)
+      error(
+        "OpSum → TTN constructor for QN conserving tensor networks only works for path/linear graphs right now.",
+      )
+    end
+    # Use `ITensors.MPO` for now until general TTN constructor
+    # works for QNs.
+    # TODO: Check it is a path graph and get a linear arrangement!
+    sites_linear_vertices = [only(sites[v]) for v in vertices(sites)]
+    vertices_to_linear_vertices = Dictionary(vertices(sites), eachindex(vertices(sites)))
+    os_linear_vertices = replace_vertices(os, vertices_to_linear_vertices)
+    mpo = MPO(os_linear_vertices, sites_linear_vertices)
+    tn = TTN(Dictionary(vertices(sites), [mpo[v] for v in 1:nv(sites)]))
+    return tn
+  end
+  T = svdTTN(os, sites, root_vertex; kwargs...)
+  if splitblocks
+    error("splitblocks not yet implemented for AbstractTreeTensorNetwork.")
+    T = ITensors.splitblocks(linkinds, T) # TODO: make this work
+  end
+  return T
+end
+
+function mpo(os::OpSum, external_inds::Vector; kwargs...)
+  return TTN(os, path_indsnetwork(external_inds); kwargs...)
+end
+
+# Conversion from other formats
+function TTN(o::Op, s::IndsNetwork; kwargs...)
+  return TTN(OpSum{Float64}() + o, s; kwargs...)
+end
+
+function TTN(o::Scaled{C,Op}, s::IndsNetwork; kwargs...) where {C}
+  return TTN(OpSum{C}() + o, s; kwargs...)
+end
+
+function TTN(o::Sum{Op}, s::IndsNetwork; kwargs...)
+  return TTN(OpSum{Float64}() + o, s; kwargs...)
+end
+
+function TTN(o::Prod{Op}, s::IndsNetwork; kwargs...)
+  return TTN(OpSum{Float64}() + o, s; kwargs...)
+end
+
+function TTN(o::Scaled{C,Prod{Op}}, s::IndsNetwork; kwargs...) where {C}
+  return TTN(OpSum{C}() + o, s; kwargs...)
+end
+
+function TTN(o::Sum{Scaled{C,Op}}, s::IndsNetwork; kwargs...) where {C}
+  return TTN(OpSum{C}() + o, s; kwargs...)
+end
+
+# Catch-all for leaf eltype specification
+function TTN(eltype::Type{<:Number}, os, sites::IndsNetwork; kwargs...)
+  return NDTensors.convert_scalartype(eltype, TTN(os, sites; kwargs...))
+end
+
+# 
+# Tree adaptation of functionalities in ITensors.jl/src/physics/autompo/matelem.jl
+# 
+
+#################################
+# ArrElem (simple sparse array) #
+#################################
+
+struct ArrElem{T,N}
+  idxs::MVector{N,Int}
+  val::T
+end
+
+function Base.:(==)(a1::ArrElem{T,N}, a2::ArrElem{T,N})::Bool where {T,N}
+  return (a1.idxs == a2.idxs && a1.val == a2.val)
+end
+
+function Base.isless(a1::ArrElem{T,N}, a2::ArrElem{T,N})::Bool where {T,N}
+  for n in 1:N
+    if a1.idxs[n] != a2.idxs[n]
+      return a1.idxs[n] < a2.idxs[n]
+    end
+  end
+  return a1.val < a2.val
+end
+
+# 
+# Sparse finite state machine construction
+# 
+
+# allow sparse arrays with ITensors.Sum entries
+function Base.zero(::Type{S}) where {S<:Sum}
+  return S()
+end
+Base.zero(t::Sum) = zero(typeof(t))
+
+"""
+    finite_state_machine(os::OpSum{C}, sites::IndsNetwork{V,<:Index}, root_vertex::V) where {C,V}
+
+Finite state machine generator for ITensors.OpSum Hamiltonian defined on a tree graph. The
+site Index graph must be a tree graph, and the chosen root  must be a leaf vertex of this
+tree. Returns a DataGraph of SparseArrayKit.SparseArrays.
+"""
+function finite_state_machine(
+  os::OpSum{C}, sites::IndsNetwork{V,<:Index}, root_vertex::V
+) where {C,V}
+  os = deepcopy(os)
+  os = sorteachterm(os, sites, root_vertex)
+  os = ITensors.sortmergeterms(os)
+
+  ValType = ITensors.determineValType(ITensors.terms(os))
+
+  # sparse symbolic representation of the TTN Hamiltonian as a DataGraph of SparseArrays
+  sparseTTN = DataGraph{V,SparseArray{Sum{Scaled{ValType,Prod{Op}}}}}(
+    underlying_graph(sites)
+  )
+
+  # traverse tree outwards from root vertex
+  vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
+  es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
+  # some things to keep track of
+  ranks = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
+  linkmaps = Dict(e => Dict{Prod{Op},Int}() for e in es)                                    # map from term in Hamiltonian to edge channel index for every edge
+  site_coef_done = Prod{Op}[]                                                               # list of Hamiltonian terms for which the coefficient has been added to a site factor
+  edge_orders = DataGraph{V,Vector{edgetype(sites)}}(underlying_graph(sites))               # relate indices of sparse TTN tensor to incident graph edges for each site
+
+  for v in vs
+    # collect all nontrivial entries of the TTN tensor at vertex v
+    entries = Tuple{MVector{ranks[v],Int},Scaled{ValType,Prod{Op}}}[]
+
+    # for every vertex, find all edges that contain this vertex
+    edges = filter(e -> dst(e) == v || src(e) == v, es)
+    # use the corresponding ordering as index order for tensor elements at this site
+    edge_orders[v] = edges
+    dim_in = findfirst(e -> dst(e) == v, edges)
+    edge_in = (isnothing(dim_in) ? [] : edges[dim_in])
+    dims_out = findall(e -> src(e) == v, edges)
+    edges_out = edges[dims_out]
+
+    # sanity check, leaves only have single incoming or outgoing edge
+    @assert !isempty(dims_out) || !isnothing(dim_in)
+    (isempty(dims_out) || isnothing(dim_in)) && @assert is_leaf(sites, v)
+
+    for term in os
+      # loop over OpSum and pick out terms that act on current vertex
+      crosses_vertex(term, sites, v) || continue
+
+      # filter out factors that come in from the direction of the incoming edge
+      incoming = filter(
+        t -> edge_in ∈ edge_path(sites, ITensors.site(t), v), ITensors.terms(term)
+      )
+      # filter out factor that acts on current vertex
+      onsite = filter(t -> (ITensors.site(t) == v), ITensors.terms(term))
+      # for every outgoing edge, filter out factors that go out along that edge
+      outgoing = Dict(
+        e => filter(t -> e ∈ edge_path(sites, v, ITensors.site(t)), ITensors.terms(term))
+        for e in edges_out
+      )
+
+      # translate into sparse tensor entry
+      T_inds = MVector{ranks[v]}(fill(-1, ranks[v]))
+      if !isnothing(dim_in) && !isempty(incoming)
+        T_inds[dim_in] = ITensors.posInLink!(linkmaps[edge_in], ITensors.argument(term)) # get incoming channel
+      end
+      for dout in dims_out
+        if !isempty(outgoing[edges[dout]])
+          T_inds[dout] = ITensors.posInLink!(linkmaps[edges[dout]], ITensors.argument(term)) # add outgoing channel
+        end
+      end
+      # if term starts at this site, add its coefficient as a site factor
+      site_coef = one(C)
+      if (isnothing(dim_in) || T_inds[dim_in] == -1) &&
+        ITensors.argument(term) ∉ site_coef_done
+        site_coef = ITensors.coefficient(term)
+        push!(site_coef_done, ITensors.argument(term))
+      end
+      # add onsite identity for interactions passing through vertex
+      if isempty(onsite)
+        if !ITensors.using_auto_fermion() && isfermionic(incoming, sites)
+          error("No verified fermion support for automatic TTN constructor!") # no verified support, just throw error
+        else
+          push!(onsite, Op("Id", v))
+        end
+      end
+      # save indices and value of sparse tensor entry
+      el = (T_inds, site_coef * Prod(onsite))
+      push!(entries, el)
+    end
+
+    # handle start and end of operator terms and convert to sparse array
+    linkdims = Tuple([
+      (isempty(linkmaps[e]) ? 0 : maximum(values(linkmaps[e]))) + 2 for e in edges
+    ])
+    T = SparseArray{Sum{Scaled{ValType,Prod{Op}}},ranks[v]}(undef, linkdims)
+    for (T_inds, t) in entries
+      if !isnothing(dim_in)
+        if T_inds[dim_in] == -1
+          T_inds[dim_in] = 1 # always start in first channel
+        else
+          T_inds[dim_in] += 1 # shift regular channel
+        end
+      end
+      if !isempty(dims_out)
+        end_dims = filter(d -> T_inds[d] == -1, dims_out)
+        normal_dims = filter(d -> T_inds[d] != -1, dims_out)
+        T_inds[end_dims] .= linkdims[end_dims] # always end in last channel
+        T_inds[normal_dims] .+= 1 # shift regular channels
+      end
+      T[T_inds...] += t
+    end
+    # add starting and ending identity operators
+    if isnothing(dim_in)
+      T[ones(Int, ranks[v])...] += one(ValType) * Prod([Op("Id", v)]) # only one real starting identity
+    end
+    # ending identities are a little more involved
+    if !isnothing(dim_in)
+      T[linkdims...] += one(ValType) * Prod([Op("Id", v)]) # place identity if all channels end
+      # place identity from start of incoming channel to start of each single outgoing channel
+      idT_end_inds = [linkdims...]
+      idT_end_inds[dim_in] = 1
+      for dout in dims_out
+        idT_end_inds[dout] = 1
+        T[idT_end_inds...] += one(ValType) * Prod([Op("Id", v)])
+        idT_end_inds[dout] = linkdims[dout] # reset
+      end
+    end
+    sparseTTN[v] = T
+  end
+  return sparseTTN, edge_orders
+end
+
+"""
+    fsmTTN(os::OpSum{C}, sites::IndsNetwork{V,<:Index}, root_vertex::V, kwargs...) where {C,V}
+
+Construct a dense TreeTensorNetwork from sparse finite state machine
+represenatation, without compression.
+"""
+function fsmTTN(
+  os::OpSum{C}, sites::IndsNetwork{V,<:Index}, root_vertex::V; trunc=false, kwargs...
+)::TTN where {C,V}
+  ValType = ITensors.determineValType(ITensors.terms(os))
+  # start from finite state machine
+  fsm, edge_orders = finite_state_machine(os, sites, root_vertex)
+  # some trickery to get link dimension for every edge
+  link_space = Dict{edgetype(sites),Index}()
+  function get_linkind!(link_space, e)
+    if !haskey(link_space, e)
+      d = findfirst(x -> (x == e || x == reverse(e)), edge_orders[src(e)])
+      link_space[e] = Index(size(fsm[src(e)], d), edge_tag(e))
+    end
+    return link_space[e]
+  end
+  # compress finite state machine into dense form
+  H = TTN(sites)
+  for v in vertices(sites)
+    linkinds = [get_linkind!(link_space, e) for e in edge_orders[v]]
+    linkdims = dim.(linkinds)
+    H[v] = ITensor()
+    for (T_ind, t) in nonzero_pairs(fsm[v])
+      any(map(x -> abs(coefficient(x)) > eps(), t)) || continue
+      T = zeros(ValType, linkdims...)
+      T[T_ind] += one(ValType)
+      T = itensor(T, linkinds)
+      H[v] += T * computeSiteSum(sites, t)
+    end
+  end
+  # add option for numerical truncation, but throw warning as this can fail sometimes
+  if trunc
+    @warn "Naive numerical truncation of TTN Hamiltonian may fail for larger systems."
+    # see https://github.com/ITensor/ITensors.jl/issues/526
+    lognormT = lognorm(H)
+    H /= exp(lognormT / nv(H)) # TODO: fix broadcasting for in-place assignment
+    H = truncate(H; root_vertex, kwargs...)
+    H *= exp(lognormT / nv(H))
+  end
+  return H
+end
+
+function computeSiteSum(
+  sites::IndsNetwork{V,<:Index}, ops::Sum{Scaled{C,Prod{Op}}}
+)::ITensor where {V,C}
+  ValType = ITensors.determineValType(ITensors.terms(ops))
+  v = ITensors.site(ITensors.argument(ops[1])[1])
+  T =
+    convert(ValType, coefficient(ops[1])) *
+    computeSiteProd(sites, ITensors.argument(ops[1]))
+  for j in 2:length(ops)
+    (ITensors.site(ITensors.argument(ops[j])[1]) != v) &&
+      error("Mismatch of vertex labels in computeSiteSum")
+    T +=
+      convert(ValType, coefficient(ops[j])) *
+      computeSiteProd(sites, ITensors.argument(ops[j]))
+  end
+  return T
+end
\ No newline at end of file

From 8ae0c5b69749af4aafaa1420ca3c6d602cfecb38 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Tue, 9 Jan 2024 17:43:57 -0500
Subject: [PATCH 24/32] Cosmetic edits, and remove an unused functions.

---
 src/treetensornetworks/opsum_to_ttn.jl | 19 ++++++++++---------
 test/test_opsum_to_ttn.jl              |  4 ++--
 2 files changed, 12 insertions(+), 11 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 61dcc70c..5ae6e5c5 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -41,7 +41,7 @@ Hamiltonian, compressing shared interaction channels.
 """
 function ttn_svd(
   os::OpSum, sites::IndsNetwork, root_vertex; kwargs...
-)::TTN
+)
   # Function barrier to improve type stability
   coefficient_type = ITensors.determineValType(terms(os))
   return ttn_svd(coefficient_type, os, sites, root_vertex; kwargs...)
@@ -55,7 +55,7 @@ function ttn_svd(
   mindim::Int=1,
   maxdim::Int=typemax(Int),
   cutoff=eps(coefficient_type) * 10,
-)::TTN
+)
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
   #FIXME: how are sites handled where some sites have QNs and some don't?
@@ -278,7 +278,7 @@ function ttn_svd(
       block_helper_inds[terminal_dims] .= 1
       for dout in filter(d -> d ∈ terminal_dims, dims_out)
         T_inds[dout] = sublinkdims[dout]                          # end in channel linkdims[d] for each dimension d
-        @assert T_inds[dout] == 1
+        @assert isone(T_inds[dout])
         block_helper_inds[dout] = nblocks(linkinds[dout])
       end
 
@@ -287,7 +287,7 @@ function ttn_svd(
         block_helper_inds[d] = qnblock(linkinds[d], T_qns[d])
       end
 
-      @assert all(block_helper_inds .!= -1) # check that all block indices are set
+      @assert all(≠(-1), block_helper_inds)# check that all block indices are set
       # make Block 
       theblock = Block(Tuple(block_helper_inds))
 
@@ -300,7 +300,7 @@ function ttn_svd(
 
       if isempty(normal_dims)
         M = get!(blocks, (theblock, terms(t)), zero_arr())
-        @assert length(M) == 1
+        @assert isone(length(M))
         M[] += ct
       else
         M = get!(blocks, (theblock, terms(t)), zero_arr())
@@ -329,7 +329,7 @@ function ttn_svd(
       ###FIXME: make this work if there's no physical degree of freedom at site
       Op = computeSiteProd(sites, Prod(q_op))  ###FIXME: is this implemented?
       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
-        (nnzblocks(Op) == 0) && continue  ###FIXME: this one may be breaking for no physical indices on site
+        iszero(nnzblocks(Op)) && continue  ###FIXME: this one may be breaking for no physical indices on site
       end
       sq = flux(Op)
       if !isnothing(sq)
@@ -496,7 +496,7 @@ function TTN(
   sites::IndsNetwork;
   root_vertex=default_root_vertex(sites),
   splitblocks=false,
-  method=:svd,
+  algorithm="svd",
   kwargs...,
 )::TTN
   length(ITensors.terms(os)) == 0 && error("OpSum has no terms")
@@ -506,7 +506,7 @@ function TTN(
   os = deepcopy(os)
   os = sorteachterm(os, sites, root_vertex)
   os = ITensors.sortmergeterms(os) # not exported
-  if method == :svd
+  if algorithm == "svd"
     T = ttn_svd(os, sites, root_vertex; kwargs...)
   else
     error("Currently only SVD is supported as TTN constructor backend.")
@@ -619,7 +619,7 @@ end
 Base.zero(t::Sum) = zero(typeof(t))
 
 
-
+#=
 function computeSiteSum(
   sites::IndsNetwork{V,<:Index}, ops::Sum{Scaled{C,Prod{Op}}}
 )::ITensor where {V,C}
@@ -637,3 +637,4 @@ function computeSiteSum(
   end
   return T
 end
+=#
diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index 03b91b2e..02c90587 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -45,7 +45,7 @@ using Test
       @test Tttno ≈ Tmpo rtol = 1e-6
 
       # this breaks for longer range interactions
-      Hsvd_lr = TTN(Hlr, is; root_vertex=root_vertex, method=:svd, cutoff=1e-10)
+      Hsvd_lr = TTN(Hlr, is; root_vertex=root_vertex,algorithm="svd", cutoff=1e-10)
       Hline_lr = MPO(relabel_sites(Hlr, vmap), sites)
       @disable_warn_order begin
         Tttno_lr = prod(Hline_lr)
@@ -115,7 +115,7 @@ using Test
       @test Tttno ≈ Tmpo rtol = 1e-6
 
       # this breaks for longer range interactions ###not anymore
-      Hsvd_lr = TTN(Hlr, is; root_vertex=root_vertex, method=:svd, cutoff=1e-10)
+      Hsvd_lr = TTN(Hlr, is; root_vertex=root_vertex, algorithm="svd", cutoff=1e-10)
       Hline_lr = MPO(relabel_sites(Hlr, vmap), sites)
       @disable_warn_order begin
         Tttno_lr = prod(Hline_lr)

From 3033544ab1436ab7223cf99ed69b165a3862b391 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Wed, 10 Jan 2024 10:40:53 -0500
Subject: [PATCH 25/32] Fix cutoff kwarg, remove commented out function.

---
 src/treetensornetworks/opsum_to_ttn.jl | 25 ++-----------------------
 1 file changed, 2 insertions(+), 23 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 5ae6e5c5..72ea5400 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -54,7 +54,7 @@ function ttn_svd(
   root_vertex;
   mindim::Int=1,
   maxdim::Int=typemax(Int),
-  cutoff=eps(coefficient_type) * 10,
+  cutoff=eps(real(coefficient_type)) * 10,
 )
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
@@ -263,7 +263,7 @@ function ttn_svd(
     blocks = Dict{Tuple{Block{degrees[v]},Vector{Op}},Array{coefficient_type,degrees[v]}}()
     for el in tempTTN[v]
       t = el.val
-      (abs(coefficient(t)) > eps()) || continue
+      (abs(coefficient(t)) > eps(real(coefficient_type))) || continue
       block_helper_inds = fill(-1, degrees[v]) # we manipulate T_inds later, and loose track of ending/starting information, so keep track of it here
       T_inds = el.idxs
       T_qns = el.qn_idxs
@@ -617,24 +617,3 @@ function Base.zero(::Type{S}) where {S<:Sum}
   return S()
 end
 Base.zero(t::Sum) = zero(typeof(t))
-
-
-#=
-function computeSiteSum(
-  sites::IndsNetwork{V,<:Index}, ops::Sum{Scaled{C,Prod{Op}}}
-)::ITensor where {V,C}
-  coefficient_type = ITensors.determineValType(ITensors.terms(ops))
-  v = ITensors.site(ITensors.argument(ops[1])[1])
-  T =
-    convert(coefficient_type, coefficient(ops[1])) *
-    computeSiteProd(sites, ITensors.argument(ops[1]))
-  for j in 2:length(ops)
-    (ITensors.site(ITensors.argument(ops[j])[1]) != v) &&
-      error("Mismatch of vertex labels in computeSiteSum")
-    T +=
-      convert(coefficient_type, coefficient(ops[j])) *
-      computeSiteProd(sites, ITensors.argument(ops[j]))
-  end
-  return T
-end
-=#

From 64e1c1bf7cd358a6ba1188114ee0c7b6aab21017 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Wed, 10 Jan 2024 11:05:53 -0500
Subject: [PATCH 26/32] Format and outdated comments.

---
 src/treetensornetworks/opsum_to_ttn.jl | 23 +++++------------------
 test/test_opsum_to_ttn.jl              |  2 +-
 2 files changed, 6 insertions(+), 19 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 72ea5400..428127e4 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -39,9 +39,7 @@ end
 Construct a dense TreeTensorNetwork from a symbolic OpSum representation of a
 Hamiltonian, compressing shared interaction channels.
 """
-function ttn_svd(
-  os::OpSum, sites::IndsNetwork, root_vertex; kwargs...
-)
+function ttn_svd(os::OpSum, sites::IndsNetwork, root_vertex; kwargs...)
   # Function barrier to improve type stability
   coefficient_type = ITensors.determineValType(terms(os))
   return ttn_svd(coefficient_type, os, sites, root_vertex; kwargs...)
@@ -210,35 +208,24 @@ function ttn_svd(
 
   # compress this tempTTN representation into dense form
 
-  #link_space = dictionary([
-  #  e => Index((isempty(outmaps[e]) ? 0 : size(Vs[e], 2)) + 2, edge_tag(e)) for e in es
-  #])
-  ###FIXME: that's an ugly constructor... ###QNIndex(undef) not defined
   link_space = Dict{edgetype_sites,Index}()
   linkdir = ITensors.using_auto_fermion() ? ITensors.In : ITensors.Out ###FIXME: ??
-  for e in es ###this one might also be tricky, in case there's no Vq defined on the edge
-    ###also tricky w.r.t. the link direction
-    ### ---> should we just infer whether the edge is incoming or outgoing and choose linkdir accordingly?
-    ### ===> this is basically decided by daggering in the MPO code (we should be able to stick to this pattern)
+  for e in es
     qi = Vector{Pair{QN,Int}}()
     push!(qi, QN() => 1)
     for (q, Vq) in Vs[e]
       cols = size(Vq, 2)
       if ITensors.using_auto_fermion() # <fermions>
-        push!(qi, (-q) => cols) ###the minus sign comes from the opposite linkdir ?
+        push!(qi, (-q) => cols)
       else
         push!(qi, q => cols)
       end
     end
     push!(qi, Hflux => 1)
-    link_space[e] = Index(qi...; tags=edge_tag(e), dir=linkdir) ##FIXME: maybe a single linkdir is not the right thing to do here
+    link_space[e] = Index(qi...; tags=edge_tag(e), dir=linkdir)
   end
 
   H = TTN(sites)
-  ##Copy from ITensors qn_svdMPO
-  # Find location where block of Index i
-  # matches QN q, but *not* 1 or dim(i)
-  # which are special ending/starting states
   function qnblock(i::Index, q::QN)
     for b in 2:(nblocks(i) - 1)
       flux(i, Block(b)) == q && return b
@@ -327,7 +314,7 @@ function ttn_svd(
 
     for ((b, q_op), m) in blocks
       ###FIXME: make this work if there's no physical degree of freedom at site
-      Op = computeSiteProd(sites, Prod(q_op))  ###FIXME: is this implemented?
+      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
         iszero(nnzblocks(Op)) && continue  ###FIXME: this one may be breaking for no physical indices on site
       end
diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index 02c90587..8669f6b2 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -45,7 +45,7 @@ using Test
       @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=root_vertex, algorithm="svd", cutoff=1e-10)
       Hline_lr = MPO(relabel_sites(Hlr, vmap), sites)
       @disable_warn_order begin
         Tttno_lr = prod(Hline_lr)

From 32b483daa616c2e00452357b4a45372d3e3066e5 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Wed, 10 Jan 2024 11:17:56 -0500
Subject: [PATCH 27/32] Format.

---
 src/treetensornetworks/deprecated_opsum_to_ttn.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/treetensornetworks/deprecated_opsum_to_ttn.jl b/src/treetensornetworks/deprecated_opsum_to_ttn.jl
index c0113c18..30e59e54 100644
--- a/src/treetensornetworks/deprecated_opsum_to_ttn.jl
+++ b/src/treetensornetworks/deprecated_opsum_to_ttn.jl
@@ -633,4 +633,4 @@ function computeSiteSum(
       computeSiteProd(sites, ITensors.argument(ops[j]))
   end
   return T
-end
\ No newline at end of file
+end

From 39c8c98951b215ce5805f3dc716cf822df36f52a Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Wed, 10 Jan 2024 15:07:56 -0500
Subject: [PATCH 28/32] Add compatibility with TTNs with internal vertices that
 have no siteinds.

---
 src/treetensornetworks/opsum_to_ttn.jl |  25 ++++-
 test/test_opsum_to_ttn.jl              | 124 +++++++++++++------------
 2 files changed, 87 insertions(+), 62 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 428127e4..14c9e985 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -48,7 +48,7 @@ end
 function ttn_svd(
   coefficient_type::Type{<:Number},
   os::OpSum,
-  sites::IndsNetwork,
+  sites0::IndsNetwork,
   root_vertex;
   mindim::Int=1,
   maxdim::Int=typemax(Int),
@@ -57,6 +57,7 @@ function ttn_svd(
   # check for qns on the site indices
   #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
   #FIXME: how are sites handled where some sites have QNs and some don't?
+  sites = deepcopy(sites0)
   edgetype_sites = edgetype(sites)
   vertextype_sites = vertextype(sites)
   thishasqns = any(v -> hasqns(sites[v]), vertices(sites))
@@ -93,6 +94,14 @@ function ttn_svd(
   end
 
   Hflux = -calc_qn(terms(first(terms(os))))
+  # insert dummy indices on internal vertices, these will not show up in the final tensor
+  is_internal = Dict{vertextype_sites,Bool}()
+  for v in vs
+    is_internal[v] = isempty(sites[v])
+    if isempty(sites[v])
+      push!(sites[v], Index(Hflux => 1)) ###FIXME: using Hflux may be a problem with AutoFermion, the dummy index has to be bosonic (but the QN of H should be too?) ...
+    end
+  end
   inbond_coefs = Dict(
     e => Dict{QN,Vector{ITensors.MatElem{coefficient_type}}}() for e in es
   )                         # bond coefficients for incoming edge channels
@@ -225,7 +234,7 @@ function ttn_svd(
     link_space[e] = Index(qi...; tags=edge_tag(e), dir=linkdir)
   end
 
-  H = TTN(sites)
+  H = TTN(sites0)   ###initialize TTN without the dummy indices added
   function qnblock(i::Index, q::QN)
     for b in 2:(nblocks(i) - 1)
       flux(i, Block(b)) == q && return b
@@ -342,7 +351,11 @@ function ttn_svd(
       if !thishasqns
         iT = removeqns(iT)
       end
-      H[v] += (iT * Op)
+      if is_internal[v]
+        H[v] += iT
+      else
+        H[v] += (iT * Op)
+      end
     end
 
     linkdims = dim.(linkinds)
@@ -367,7 +380,11 @@ function ttn_svd(
     if !thishasqns
       T = removeqns(T)
     end
-    H[v] += T * ITensorNetworks.computeSiteProd(sites, Prod([(Op("Id", v))]))
+    if is_internal[v]
+      H[v] += T
+    else
+      H[v] += T * ITensorNetworks.computeSiteProd(sites, Prod([(Op("Id", v))]))
+    end
   end
 
   return H
diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index 8669f6b2..dfe615e9 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -24,8 +24,8 @@ using Test
     H = ising(c; J1=J1, J2=J2, h=h)
     # add combination of longer range interactions
     Hlr = copy(H)
-    Hlr += 5, "Z", (1, 2), "Z", (2, 2)
-    Hlr += -4, "Z", (1, 1), "Z", (2, 2)
+    Hlr += 5, "Z", (1, 2), "Z", (2, 2), "Z", (3, 2)
+    Hlr += -4, "Z", (1, 1), "Z", (2, 2), "Z", (3, 1)
     Hlr += 2.0, "Z", (2, 2), "Z", (3, 2)
     Hlr += -1.0, "Z", (1, 2), "Z", (3, 1)
 
@@ -93,7 +93,7 @@ using Test
     H = heisenberg(c; J1=J1, J2=J2, h=h)
     # add combination of longer range interactions
     Hlr = copy(H)
-    Hlr += 5, "Z", (1, 2), "Z", (2, 2)
+    Hlr += 5, "Z", (1, 2), "Z", (2, 2)#, "Z", (3,2)
     Hlr += -4, "Z", (1, 1), "Z", (2, 2)
     Hlr += 2.0, "Z", (2, 2), "Z", (3, 2)
     Hlr += -1.0, "Z", (1, 2), "Z", (3, 1)
@@ -124,63 +124,71 @@ using Test
       @test Tttno_lr ≈ Tmpo_lr rtol = 1e-6
     end
   end
-  #=
-    @testset "OpSum to TTN QN missing" begin
-      # small comb tree
-      tooth_lengths = fill(2, 3)
-      c = named_comb_tree(tooth_lengths)
-      is = siteinds("S=1/2", c; conserve_qns=true)
-      is_missing_site=copy(is)
-      is_missing_site[(2,1)]=Vector{Index}[]
-      is_noqns = copy(is)
-      for v in vertices(is)
-        is_noqns[v] = removeqns(is_noqns[v])
+
+  @testset "OpSum to TTN QN missing" begin
+    # small comb tree
+    tooth_lengths = fill(2, 3)
+    c = named_comb_tree(tooth_lengths)
+    c2 = copy(c)
+    import ITensorNetworks: add_vertex!, rem_vertex!, add_edge!, rem_edge!
+    add_vertex!(c2, (-1, 1))
+    add_edge!(c2, (-1, 1) => (2, 2))
+    add_edge!(c2, (-1, 1) => (3, 1))
+    add_edge!(c2, (-1, 1) => (2, 1))
+    rem_edge!(c2, (2, 1) => (2, 2))
+    rem_edge!(c2, (2, 1) => (3, 1))
+    @show c2
+
+    is = siteinds("S=1/2", c; conserve_qns=true)
+    is_missing_site = siteinds("S=1/2", c2; conserve_qns=true)
+    is_missing_site[(-1, 1)] = Vector{Index}[]
+
+    # linearized version
+    linear_order = [4, 1, 2, 5, 3, 6]
+    vmap = Dictionary(vertices(is)[linear_order], 1:length(linear_order))
+    @show filter(d -> !isempty(d), vertex_data(is_missing_site))
+    sites = only.(filter(d -> !isempty(d), collect(vertex_data(is_missing_site))))[linear_order]
+
+    J1 = -1
+    J2 = 2
+    h = 0.5
+    H = heisenberg(c; J1=J1, J2=J2, h=h)
+    #Hvac = heisenberg(is_missing_site; J1=J1, J2=J2, h=h)
+
+    # add combination of longer range interactions
+    Hlr = copy(H)
+    Hlr += 5, "Z", (1, 2), "Z", (2, 2)
+    Hlr += -4, "Z", (1, 1), "Z", (2, 2)
+    Hlr += 2.0, "Z", (2, 2), "Z", (3, 2)
+    Hlr += -1.0, "Z", (1, 2), "Z", (3, 1)
+
+    # root_vertex = (1, 2)
+    # println(leaf_vertices(is))
+
+    @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)
+      # get corresponding MPO Hamiltonian
+      Hline = MPO(relabel_sites(H, vmap), sites)
+      # compare resulting sparse Hamiltonians
+
+      @disable_warn_order begin
+        Tmpo = prod(Hline)
+        Tttno = contract(Hsvd)
       end
 
-      # linearized version
-      linear_order = [4, 1, 2, 5, 3, 6]
-      vmap = Dictionary(vertices(is)[linear_order], 1:length(linear_order))
-      sites = only.(collect(vertex_data(is)))[linear_order]
-
-      # test with next-to-nearest-neighbor Ising Hamiltonian
-      J1 = -1
-      J2 = 2
-      h = 0.5
-      H = heisenberg(c; J1=J1, J2=J2, h=h)
-      Hvac = heisenberg(is_missing_site; J1=J1, J2=J2, h=h)
-
-      # add combination of longer range interactions
-      Hlr = copy(H)
-      Hlr += 5, "Z", (1, 2), "Z", (2, 2)
-      Hlr += -4, "Z", (1, 1), "Z", (2, 2)
-      Hlr += 2.0, "Z", (2, 2), "Z", (3, 2)
-      Hlr += -1.0, "Z", (1, 2), "Z", (3, 1)
-
-      # root_vertex = (1, 2)
-      # println(leaf_vertices(is))
-
-      @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)
-        # get corresponding MPO Hamiltonian
-        Hline = MPO(relabel_sites(H, vmap), sites)
-        # compare resulting sparse Hamiltonians
-
-        @disable_warn_order begin
-          Tmpo = prod(Hline)
-          Tttno = contract(Hsvd)
-        end
-        @test Tttno ≈ Tmpo rtol = 1e-6
-
-        # this breaks for longer range interactions ###not anymore
-        Hsvd_lr = TTN(Hlr, is; root_vertex=root_vertex, method=:svd, cutoff=1e-10)
-        Hline_lr = MPO(relabel_sites(Hlr, vmap), sites)
-        @disable_warn_order begin
-          Tttno_lr = prod(Hline_lr)
-          Tmpo_lr = contract(Hsvd_lr)
-        end
-        @test Tttno_lr ≈ Tmpo_lr rtol = 1e-6
+      @test Tttno ≈ Tmpo rtol = 1e-6
+
+      # this breaks for longer range interactions ###not anymore
+      Hsvd_lr = TTN(
+        Hlr, is_missing_site; root_vertex=root_vertex, algorithm="svd", cutoff=1e-10
+      )
+      Hline_lr = MPO(relabel_sites(Hlr, vmap), sites)
+      @disable_warn_order begin
+        Tttno_lr = prod(Hline_lr)
+        Tmpo_lr = contract(Hsvd_lr)
       end
+      @test Tttno_lr ≈ Tmpo_lr rtol = 1e-6
     end
-    =#
+  end
 end

From 0a786040f4e53a2fa8593a72e74c3dd2629b5b89 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Wed, 10 Jan 2024 15:10:28 -0500
Subject: [PATCH 29/32] Remove unnecessary comments.

---
 src/treetensornetworks/opsum_to_ttn.jl | 21 ++++++---------------
 1 file changed, 6 insertions(+), 15 deletions(-)

diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl
index 14c9e985..09703c2c 100644
--- a/src/treetensornetworks/opsum_to_ttn.jl
+++ b/src/treetensornetworks/opsum_to_ttn.jl
@@ -54,23 +54,17 @@ function ttn_svd(
   maxdim::Int=typemax(Int),
   cutoff=eps(real(coefficient_type)) * 10,
 )
-  # check for qns on the site indices
-  #FIXME: this check for whether or not any of the siteindices has QNs is somewhat ugly
-  #FIXME: how are sites handled where some sites have QNs and some don't?
   sites = deepcopy(sites0)
   edgetype_sites = edgetype(sites)
   vertextype_sites = vertextype(sites)
   thishasqns = any(v -> hasqns(sites[v]), vertices(sites))
 
   # traverse tree outwards from root vertex
-  ### FIXME: generalize to internal vertices
   vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
   es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
   # 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(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   
-  #outmaps = Dict(e => Dict{Pair{Vector{Op},QN},Int}() for e in es)                                   
   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
 
@@ -175,19 +169,17 @@ function ttn_svd(
           outmaps, edges[dout] => -outgoing_qns[edges[dout]], Dict{Vector{Op},Int}()
         )
         T_inds[dout] = ITensors.posInLink!(coutmap, outgoing[edges[dout]]) # add outgoing channel
-        T_qns[dout] = -outgoing_qns[edges[dout]]  ###minus sign to account for outgoing arrow direction
+        T_qns[dout] = -outgoing_qns[edges[dout]]  # minus sign to account for outgoing arrow direction
       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)
-        # FIXME: maybe put a try block here for a more helpful error message?
-        site_coef = convert(coefficient_type, site_coef) #required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real
+        site_coef = convert(coefficient_type, site_coef) # required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real
         push!(site_coef_done, ITensors.argument(term))
       end
       # add onsite identity for interactions passing through vertex
-      ### FIXME: generalize to internal vertices, i.e. do not add Id if this is an internal site
       if isempty(onsite)
         if !ITensors.using_auto_fermion() && isfermionic(incoming, sites)
           error("No verified fermion support for automatic TTN constructor!")
@@ -209,8 +201,8 @@ function ttn_svd(
         P = S .^ 2
         truncate!(P; maxdim, cutoff, mindim)
         tdim = length(P)
-        nc = size(M, 2) ###CHECK: verify whether this is still correct, diverges from opsum_to_mpo_qn
-        Vs[edges[dim_in]][q] = Matrix{coefficient_type}(V[1:nc, 1:tdim]) ###CHECK: verify whether this is still correct, diverges from opsum_to_mpo_qn
+        nc = size(M, 2)
+        Vs[edges[dim_in]][q] = Matrix{coefficient_type}(V[1:nc, 1:tdim])
       end
     end
   end
@@ -234,7 +226,7 @@ function ttn_svd(
     link_space[e] = Index(qi...; tags=edge_tag(e), dir=linkdir)
   end
 
-  H = TTN(sites0)   ###initialize TTN without the dummy indices added
+  H = TTN(sites0)   # initialize TTN without the dummy indices added
   function qnblock(i::Index, q::QN)
     for b in 2:(nblocks(i) - 1)
       flux(i, Block(b)) == q && return b
@@ -322,10 +314,9 @@ function ttn_svd(
     end
 
     for ((b, q_op), m) in blocks
-      ###FIXME: make this work if there's no physical degree of freedom at site
       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
-        iszero(nnzblocks(Op)) && continue  ###FIXME: this one may be breaking for no physical indices on site
+        iszero(nnzblocks(Op)) && continue
       end
       sq = flux(Op)
       if !isnothing(sq)

From a0e9fb4f45cb89d4f77fa5ccc880ba523d5221b9 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Wed, 10 Jan 2024 15:31:23 -0500
Subject: [PATCH 30/32] Remove output from tests.

---
 test/test_opsum_to_ttn.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index dfe615e9..7d986ece 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -137,7 +137,7 @@ using Test
     add_edge!(c2, (-1, 1) => (2, 1))
     rem_edge!(c2, (2, 1) => (2, 2))
     rem_edge!(c2, (2, 1) => (3, 1))
-    @show c2
+    #@show c2
 
     is = siteinds("S=1/2", c; conserve_qns=true)
     is_missing_site = siteinds("S=1/2", c2; conserve_qns=true)
@@ -146,7 +146,6 @@ using Test
     # linearized version
     linear_order = [4, 1, 2, 5, 3, 6]
     vmap = Dictionary(vertices(is)[linear_order], 1:length(linear_order))
-    @show filter(d -> !isempty(d), vertex_data(is_missing_site))
     sites = only.(filter(d -> !isempty(d), collect(vertex_data(is_missing_site))))[linear_order]
 
     J1 = -1

From 149372fff1399b185ae148dfc5ab204f95644b28 Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Wed, 10 Jan 2024 17:17:21 -0500
Subject: [PATCH 31/32] Use ITensorNetworks.model(...) instead of model(...) in
 test.

---
 test/test_opsum_to_ttn.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/test_opsum_to_ttn.jl b/test/test_opsum_to_ttn.jl
index 7d986ece..33027066 100644
--- a/test/test_opsum_to_ttn.jl
+++ b/test/test_opsum_to_ttn.jl
@@ -21,7 +21,7 @@ using Test
     J1 = -1
     J2 = 2
     h = 0.5
-    H = ising(c; J1=J1, J2=J2, h=h)
+    H = ITensorNetworks.ising(c; J1=J1, J2=J2, h=h)
     # add combination of longer range interactions
     Hlr = copy(H)
     Hlr += 5, "Z", (1, 2), "Z", (2, 2), "Z", (3, 2)
@@ -90,7 +90,7 @@ using Test
     J1 = -1
     J2 = 2
     h = 0.5
-    H = heisenberg(c; J1=J1, J2=J2, h=h)
+    H = ITensorNetworks.heisenberg(c; J1=J1, J2=J2, h=h)
     # add combination of longer range interactions
     Hlr = copy(H)
     Hlr += 5, "Z", (1, 2), "Z", (2, 2)#, "Z", (3,2)
@@ -151,7 +151,7 @@ using Test
     J1 = -1
     J2 = 2
     h = 0.5
-    H = heisenberg(c; J1=J1, J2=J2, h=h)
+    H = ITensorNetworks.heisenberg(c; J1=J1, J2=J2, h=h)
     #Hvac = heisenberg(is_missing_site; J1=J1, J2=J2, h=h)
 
     # add combination of longer range interactions

From 66d9103c3803baa45e411742918b2b78600121eb Mon Sep 17 00:00:00 2001
From: Benedikt Kloss <bene.kloss@gmail.com>
Date: Thu, 11 Jan 2024 11:25:31 -0500
Subject: [PATCH 32/32] Remove @test_broken related to lack of QN
 implementation for TTN constructor.

---
 test/test_treetensornetworks/test_solvers/test_dmrg_x.jl | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl
index 8d80e162..65871b56 100644
--- a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl
+++ b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl
@@ -49,12 +49,6 @@ end
   # Random fields h ∈ [-W, W]
   h = W * (2 * rand(nv(c)) .- 1)
 
-  if conserve_qns
-    # QN conservation for non-path graphs is currently broken
-    @test_broken TTN(ITensorNetworks.heisenberg(c; h), s)
-    continue
-  end
-
   H = TTN(ITensorNetworks.heisenberg(c; h), s)
 
   # TODO: Use `TTN(s; states=v -> rand(["↑", "↓"]))` or